| s_expl.c (241051) | s_expl.c (241516) |
|---|---|
| 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 241051 2012-09-29 16:40:12Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 241516 2012-10-13 19:53:11Z 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 * --- 20 unchanged lines hidden (view full) --- 59static volatile const long double tiny = 0x1p-10000L; 60 61static const union IEEEl2bits 62/* log(2**16384 - 0.5) rounded towards zero: */ 63o_threshold = LD80C(0xb17217f7d1cf79ab, 13, 11356.5234062941439488L), 64/* log(2**(-16381-64-1)) rounded towards zero: */ 65u_threshold = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L); 66 | 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 * --- 20 unchanged lines hidden (view full) --- 59static volatile const long double tiny = 0x1p-10000L; 60 61static const union IEEEl2bits 62/* log(2**16384 - 0.5) rounded towards zero: */ 63o_threshold = LD80C(0xb17217f7d1cf79ab, 13, 11356.5234062941439488L), 64/* log(2**(-16381-64-1)) rounded towards zero: */ 65u_threshold = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L); 66 |
| 67static const double __aligned(64) | 67static const double |
| 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 */ 74L1 = 5.4152123484527692e-3, /* 0x162e42ff000000.0p-60 */ 75L2 = -3.2819649005320973e-13, /* -0x1718432a1b0e26.0p-94 */ --- 5 unchanged lines hidden (view full) --- 81P2 = 0.5, 82P3 = 1.6666666666666119e-1, /* 0x15555555555490.0p-55 */ 83P4 = 4.1666666666665887e-2, /* 0x155555555554e5.0p-57 */ 84P5 = 8.3333354987869413e-3, /* 0x1111115b789919.0p-59 */ 85P6 = 1.3888891738560272e-3; /* 0x16c16c651633ae.0p-62 */ 86 87/* 88 * 2^(i/INTERVALS) for i in [0,INTERVALS] is represented by two values where | 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 */ 74L1 = 5.4152123484527692e-3, /* 0x162e42ff000000.0p-60 */ 75L2 = -3.2819649005320973e-13, /* -0x1718432a1b0e26.0p-94 */ --- 5 unchanged lines hidden (view full) --- 81P2 = 0.5, 82P3 = 1.6666666666666119e-1, /* 0x15555555555490.0p-55 */ 83P4 = 4.1666666666665887e-2, /* 0x155555555554e5.0p-57 */ 84P5 = 8.3333354987869413e-3, /* 0x1111115b789919.0p-59 */ 85P6 = 1.3888891738560272e-3; /* 0x16c16c651633ae.0p-62 */ 86 87/* 88 * 2^(i/INTERVALS) for i in [0,INTERVALS] is represented by two values where |
| 89 * the first 53 bits of the significand is stored in hi and the next 53 90 * bits are in lo. Tang's paper states that the trailing 6 bits of hi should | 89 * the first 53 bits of the significand are stored in hi and the next 53 90 * bits are in lo. Tang's paper states that the trailing 6 bits of hi must |
| 91 * be zero for his algorithm in both single and double precision, because 92 * the table is re-used in the implementation of expm1() where a floating | 91 * be zero for his algorithm in both single and double precision, because 92 * the table is re-used in the implementation of expm1() where a floating |
| 93 * point addition involving hi must be exact. The conversion of a 53-bit 94 * double into a 64-bit long double gives 11 trailing bit, which are zero. | 93 * point addition involving hi must be exact. Here hi is double, so 94 * converting it to long double gives 11 trailing zero bits. |
| 95 */ 96static const struct { 97 double hi; 98 double lo; | 95 */ 96static const struct { 97 double hi; 98 double lo; |
| 99} s[INTERVALS] __aligned(16) = { | 99/* XXX should rename 's'. */ 100} s[INTERVALS] = { |
| 100 0x1p+0, 0x0p+0, 101 0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54, 102 0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53, 103 0x1.04315e86e7f84p+0, 0x1.7ae71f3441b49p-53, 104 0x1.059b0d3158574p+0, 0x1.d73e2a475b465p-55, 105 0x1.0706b29ddf6ddp+0, 0x1.8db880753b0f6p-53, 106 0x1.0874518759bc8p+0, 0x1.186be4bb284ffp-57, 107 0x1.09e3ecac6f383p+0, 0x1.1487818316136p-54, --- 196 unchanged lines hidden --- | 101 0x1p+0, 0x0p+0, 102 0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54, 103 0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53, 104 0x1.04315e86e7f84p+0, 0x1.7ae71f3441b49p-53, 105 0x1.059b0d3158574p+0, 0x1.d73e2a475b465p-55, 106 0x1.0706b29ddf6ddp+0, 0x1.8db880753b0f6p-53, 107 0x1.0874518759bc8p+0, 0x1.186be4bb284ffp-57, 108 0x1.09e3ecac6f383p+0, 0x1.1487818316136p-54, --- 196 unchanged lines hidden --- |