| s_expl.c (251338) | s_expl.c (251339) |
|---|---|
| 1/*- 2 * Copyright (c) 2009-2013 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 13 unchanged lines hidden (view full) --- 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 * 26 * Optimized by Bruce D. Evans. 27 */ 28 29#include <sys/cdefs.h> | 1/*- 2 * Copyright (c) 2009-2013 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 13 unchanged lines hidden (view full) --- 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 * 26 * Optimized by Bruce D. Evans. 27 */ 28 29#include <sys/cdefs.h> |
| 30__FBSDID("$FreeBSD: head/lib/msun/ld128/s_expl.c 251338 2013-06-03 18:57:35Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld128/s_expl.c 251339 2013-06-03 19:13:44Z kargl $"); |
| 31 32/* 33 * ld128 version of s_expl.c. See ../ld80/s_expl.c for most comments. 34 */ 35 36#include <float.h> 37 38#include "fpmath.h" --- 184 unchanged lines hidden (view full) --- 223long double 224expl(long double x) 225{ 226 union IEEEl2bits u, v; 227 long double q, r, r1, t, twopk, twopkp10000; 228 double dr, fn, r2; 229 230 int k, n, n2; | 31 32/* 33 * ld128 version of s_expl.c. See ../ld80/s_expl.c for most comments. 34 */ 35 36#include <float.h> 37 38#include "fpmath.h" --- 184 unchanged lines hidden (view full) --- 223long double 224expl(long double x) 225{ 226 union IEEEl2bits u, v; 227 long double q, r, r1, t, twopk, twopkp10000; 228 double dr, fn, r2; 229 230 int k, n, n2; |
| 231 uint32_t hx, ix; | 231 uint16_t hx, ix; |
| 232 233 /* Filter out exceptional cases. */ 234 u.e = x; 235 hx = u.xbits.expsign; 236 ix = hx & 0x7fff; 237 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 238 if (ix == BIAS + LDBL_MAX_EXP) { 239 if (hx & 0x8000) /* x is -Inf or -NaN */ 240 return (-1 / x); 241 return (x + x); /* x is +Inf or NaN */ 242 } 243 if (x > o_threshold) 244 return (huge * huge); 245 if (x < u_threshold) 246 return (tiny * tiny); 247 } else if (ix < BIAS - 114) { /* |x| < 0x1p-114 */ 248 return (1 + x); /* 1 with inexact iff x != 0 */ 249 } 250 | 232 233 /* Filter out exceptional cases. */ 234 u.e = x; 235 hx = u.xbits.expsign; 236 ix = hx & 0x7fff; 237 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 238 if (ix == BIAS + LDBL_MAX_EXP) { 239 if (hx & 0x8000) /* x is -Inf or -NaN */ 240 return (-1 / x); 241 return (x + x); /* x is +Inf or NaN */ 242 } 243 if (x > o_threshold) 244 return (huge * huge); 245 if (x < u_threshold) 246 return (tiny * tiny); 247 } else if (ix < BIAS - 114) { /* |x| < 0x1p-114 */ 248 return (1 + x); /* 1 with inexact iff x != 0 */ 249 } 250 |
| 251 ENTERI(); 252 |
|
| 251 /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */ 252 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 253 /* XXX assume no extra precision for the additions, as for trig fns. */ 254 /* XXX this set of comments is now quadruplicated. */ 255 fn = (double)x * INV_L + 0x1.8p52 - 0x1.8p52; 256#if defined(HAVE_EFFICIENT_IRINT) 257 n = irint(fn); 258#else 259 n = (int)fn; 260#endif 261 n2 = (unsigned)n % INTERVALS; 262 k = n >> LOG2_INTERVALS; 263 r1 = x - fn * L1; 264 r2 = fn * -L2; | 253 /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */ 254 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 255 /* XXX assume no extra precision for the additions, as for trig fns. */ 256 /* XXX this set of comments is now quadruplicated. */ 257 fn = (double)x * INV_L + 0x1.8p52 - 0x1.8p52; 258#if defined(HAVE_EFFICIENT_IRINT) 259 n = irint(fn); 260#else 261 n = (int)fn; 262#endif 263 n2 = (unsigned)n % INTERVALS; 264 k = n >> LOG2_INTERVALS; 265 r1 = x - fn * L1; 266 r2 = fn * -L2; |
| 267 r = r1 + r2; |
|
| 265 266 /* Prepare scale factors. */ | 268 269 /* Prepare scale factors. */ |
| 267 v.xbits.manh = 0; 268 v.xbits.manl = 0; | 270 /* XXX sparc64 multiplication is so slow that scalbnl() is faster. */ 271 v.e = 1; |
| 269 if (k >= LDBL_MIN_EXP) { 270 v.xbits.expsign = BIAS + k; 271 twopk = v.e; 272 } else { 273 v.xbits.expsign = BIAS + k + 10000; 274 twopkp10000 = v.e; 275 } 276 277 /* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */ 278 dr = r; 279 q = r2 + r * r * (A2 + r * (A3 + r * (A4 + r * (A5 + r * (A6 + 280 dr * (A7 + dr * (A8 + dr * (A9 + dr * A10)))))))); 281 t = tbl[n2].lo + tbl[n2].hi; 282 t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi; 283 284 /* Scale by 2**k. */ 285 if (k >= LDBL_MIN_EXP) { 286 if (k == LDBL_MAX_EXP) | 272 if (k >= LDBL_MIN_EXP) { 273 v.xbits.expsign = BIAS + k; 274 twopk = v.e; 275 } else { 276 v.xbits.expsign = BIAS + k + 10000; 277 twopkp10000 = v.e; 278 } 279 280 /* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */ 281 dr = r; 282 q = r2 + r * r * (A2 + r * (A3 + r * (A4 + r * (A5 + r * (A6 + 283 dr * (A7 + dr * (A8 + dr * (A9 + dr * A10)))))))); 284 t = tbl[n2].lo + tbl[n2].hi; 285 t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi; 286 287 /* Scale by 2**k. */ 288 if (k >= LDBL_MIN_EXP) { 289 if (k == LDBL_MAX_EXP) |
| 287 return (t * 2.0L * 0x1p16383L); 288 return (t * twopk); | 290 RETURNI(t * 2 * 0x1p16383L); 291 RETURNI(t * twopk); |
| 289 } else { | 292 } else { |
| 290 return (t * twopkp10000 * twom10000); | 293 RETURNI(t * twopkp10000 * twom10000); |
| 291 } 292} | 294 } 295} |