| s_expl.c (251327) | s_expl.c (251328) |
|---|---|
| 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/ld80/s_expl.c 251327 2013-06-03 17:51:08Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 251328 2013-06-03 18:07:04Z kargl $"); |
| 31 32/** 33 * Compute the exponential of x for Intel 80-bit format. This is based on: 34 * 35 * PTP Tang, "Table-driven implementation of the exponential function 36 * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 37 * 144-157 (1989). 38 * --- 17 unchanged lines hidden (view full) --- 56static const long double 57huge = 0x1p10000L, 58twom10000 = 0x1p-10000L; 59/* XXX Prevent gcc from erroneously constant folding this: */ 60static volatile const long double tiny = 0x1p-10000L; 61 62static const union IEEEl2bits 63/* log(2**16384 - 0.5) rounded towards zero: */ | 31 32/** 33 * Compute the exponential of x for Intel 80-bit format. This is based on: 34 * 35 * PTP Tang, "Table-driven implementation of the exponential function 36 * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 37 * 144-157 (1989). 38 * --- 17 unchanged lines hidden (view full) --- 56static const long double 57huge = 0x1p10000L, 58twom10000 = 0x1p-10000L; 59/* XXX Prevent gcc from erroneously constant folding this: */ 60static volatile const long double tiny = 0x1p-10000L; 61 62static const union IEEEl2bits 63/* log(2**16384 - 0.5) rounded towards zero: */ |
| 64o_threshold = LD80C(0xb17217f7d1cf79ab, 13, 11356.5234062941439488L), | 64/* log(2**16384 - 0.5 + 1) rounded towards zero for expm1l() is the same: */ 65o_thresholdu = LD80C(0xb17217f7d1cf79ab, 13, 11356.5234062941439488L), 66#define o_threshold (o_thresholdu.e) |
| 65/* log(2**(-16381-64-1)) rounded towards zero: */ | 67/* log(2**(-16381-64-1)) rounded towards zero: */ |
| 66u_threshold = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L); | 68u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L); 69#define u_threshold (u_thresholdu.e) |
| 67 68static const double 69/* 70 * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication). L1 must 71 * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest 72 * bits zero so that multiplication of it by n is exact. 73 */ 74INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */ --- 166 unchanged lines hidden (view full) --- 241 hx = u.xbits.expsign; 242 ix = hx & 0x7fff; 243 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 244 if (ix == BIAS + LDBL_MAX_EXP) { 245 if (hx & 0x8000 && u.xbits.man == 1ULL << 63) 246 return (0.0L); /* x is -Inf */ 247 return (x + x); /* x is +Inf, NaN or unsupported */ 248 } | 70 71static const double 72/* 73 * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication). L1 must 74 * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest 75 * bits zero so that multiplication of it by n is exact. 76 */ 77INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */ --- 166 unchanged lines hidden (view full) --- 244 hx = u.xbits.expsign; 245 ix = hx & 0x7fff; 246 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 247 if (ix == BIAS + LDBL_MAX_EXP) { 248 if (hx & 0x8000 && u.xbits.man == 1ULL << 63) 249 return (0.0L); /* x is -Inf */ 250 return (x + x); /* x is +Inf, NaN or unsupported */ 251 } |
| 249 if (x > o_threshold.e) | 252 if (x > o_threshold) |
| 250 return (huge * huge); | 253 return (huge * huge); |
| 251 if (x < u_threshold.e) | 254 if (x < u_threshold) |
| 252 return (tiny * tiny); 253 } else if (ix < BIAS - 66) { /* |x| < 0x1p-66 */ 254 /* includes pseudo-denormals */ 255 if (huge + x > 1.0L) /* trigger inexact iff x != 0 */ 256 return (1.0L + x); 257 } 258 259 ENTERI(); --- 46 unchanged lines hidden --- | 255 return (tiny * tiny); 256 } else if (ix < BIAS - 66) { /* |x| < 0x1p-66 */ 257 /* includes pseudo-denormals */ 258 if (huge + x > 1.0L) /* trigger inexact iff x != 0 */ 259 return (1.0L + x); 260 } 261 262 ENTERI(); --- 46 unchanged lines hidden --- |