| s_expl.c (240861) | s_expl.c (240864) |
|---|---|
| 1/*- 2 * Copyright (c) 2009-2012 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-2012 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 240861 2012-09-23 17:36:01Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 240864 2012-09-23 18:06:27Z kargl $"); |
| 31 | 31 |
| 32/* | 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 * 39 * where the 32 table entries have been expanded to INTERVALS (see below). 40 */ --- 24 unchanged lines hidden (view full) --- 65u_threshold = LD80C(0xb21dfe7f09e2baa9, 13, 1, -11399.4985314888605581L); 66 67static const double __aligned(64) 68/* 69 * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication). L1 must 70 * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest 71 * bits zero so that multiplication of it by n is exact. 72 */ | 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 * 39 * where the 32 table entries have been expanded to INTERVALS (see below). 40 */ --- 24 unchanged lines hidden (view full) --- 65u_threshold = LD80C(0xb21dfe7f09e2baa9, 13, 1, -11399.4985314888605581L); 66 67static const double __aligned(64) 68/* 69 * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication). L1 must 70 * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest 71 * bits zero so that multiplication of it by n is exact. 72 */ |
| 73INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */ |
|
| 73L1 = 5.4152123484527692e-3, /* 0x162e42ff000000.0p-60 */ 74L2 = -3.2819649005320973e-13, /* -0x1718432a1b0e26.0p-94 */ | 74L1 = 5.4152123484527692e-3, /* 0x162e42ff000000.0p-60 */ 75L2 = -3.2819649005320973e-13, /* -0x1718432a1b0e26.0p-94 */ |
| 75INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */ | |
| 76/* 77 * Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]: 78 * |exp(x) - p(x)| < 2**-77.2 79 * (0.002708 is ln2/(2*INTERVALS) rounded up a little). 80 */ 81P2 = 0.5, 82P3 = 1.6666666666666119e-1, /* 0x15555555555490.0p-55 */ 83P4 = 4.1666666666665887e-2, /* 0x155555555554e5.0p-57 */ --- 142 unchanged lines hidden (view full) --- 226 0x1.fa7c1819e90d8p+0, 0x1.74853f3a5931ep-55, 227 0x1.fd3c22b8f71f1p+0, 0x1.2eb74966579e7p-57 228}; 229 230long double 231expl(long double x) 232{ 233 union IEEEl2bits u, v; | 76/* 77 * Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]: 78 * |exp(x) - p(x)| < 2**-77.2 79 * (0.002708 is ln2/(2*INTERVALS) rounded up a little). 80 */ 81P2 = 0.5, 82P3 = 1.6666666666666119e-1, /* 0x15555555555490.0p-55 */ 83P4 = 4.1666666666665887e-2, /* 0x155555555554e5.0p-57 */ --- 142 unchanged lines hidden (view full) --- 226 0x1.fa7c1819e90d8p+0, 0x1.74853f3a5931ep-55, 227 0x1.fd3c22b8f71f1p+0, 0x1.2eb74966579e7p-57 228}; 229 230long double 231expl(long double x) 232{ 233 union IEEEl2bits u, v; |
| 234 long double fn, r, r1, r2, q, t, t23, t45, twopk, twopkp10000, z; | 234 long double fn, q, r, r1, r2, t, t23, t45, twopk, twopkp10000, z; |
| 235 int k, n, n2; 236 uint16_t hx, ix; 237 238 /* Filter out exceptional cases. */ 239 u.e = x; 240 hx = u.xbits.expsign; 241 ix = hx & 0x7fff; 242 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ --- 20 unchanged lines hidden (view full) --- 263 r = x - fn * L1 - fn * L2; /* r = r1 + r2 done independently. */ 264#if defined(HAVE_EFFICIENT_IRINTL) 265 n = irintl(fn); 266#elif defined(HAVE_EFFICIENT_IRINT) 267 n = irint(fn); 268#else 269 n = (int)fn; 270#endif | 235 int k, n, n2; 236 uint16_t hx, ix; 237 238 /* Filter out exceptional cases. */ 239 u.e = x; 240 hx = u.xbits.expsign; 241 ix = hx & 0x7fff; 242 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ --- 20 unchanged lines hidden (view full) --- 263 r = x - fn * L1 - fn * L2; /* r = r1 + r2 done independently. */ 264#if defined(HAVE_EFFICIENT_IRINTL) 265 n = irintl(fn); 266#elif defined(HAVE_EFFICIENT_IRINT) 267 n = irint(fn); 268#else 269 n = (int)fn; 270#endif |
| 271 n2 = (unsigned)n % INTERVALS; /* Tang's j. */ | 271 n2 = (unsigned)n % INTERVALS; |
| 272 k = (n - n2) / INTERVALS; 273 r1 = x - fn * L1; 274 r2 = -fn * L2; 275 276 /* Prepare scale factors. */ 277 v.xbits.man = 1ULL << 63; 278 if (k >= LDBL_MIN_EXP) { 279 v.xbits.expsign = BIAS + k; --- 24 unchanged lines hidden --- | 272 k = (n - n2) / INTERVALS; 273 r1 = x - fn * L1; 274 r2 = -fn * L2; 275 276 /* Prepare scale factors. */ 277 v.xbits.man = 1ULL << 63; 278 if (k >= LDBL_MIN_EXP) { 279 v.xbits.expsign = BIAS + k; --- 24 unchanged lines hidden --- |