1
2 /* @(#)e_exp.c 5.1 93/09/24 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14 /* exp(x)
15 * Returns the exponential of x.
16 *
17 * Method
18 * 1. Argument reduction:
19 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
20 * Given x, find r and integer k such that
21 *
22 * x = k*ln2 + r, |r| <= 0.5*ln2.
23 *
24 * Here r will be represented as r = hi-lo for better
25 * accuracy.
26 *
27 * 2. Approximation of exp(r) by a special rational function on
28 * the interval [0,0.34658]:
29 * Write
30 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
31 * We use a special Reme algorithm on [0,0.34658] to generate
32 * a polynomial of degree 5 to approximate R. The maximum error
33 * of this polynomial approximation is bounded by 2**-59. In
34 * other words,
35 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
36 * (where z=r*r, and the values of P1 to P5 are listed below)
37 * and
38 * | 5 | -59
39 * | 2.0+P1*z+...+P5*z - R(z) | <= 2
40 * | |
41 * The computation of exp(r) thus becomes
42 * 2*r
43 * exp(r) = 1 + -------
44 * R - r
45 * r*R1(r)
46 * = 1 + r + ----------- (for better accuracy)
47 * 2 - R1(r)
48 * where
49 * 2 4 10
50 * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
51 *
52 * 3. Scale back to obtain exp(x):
53 * From step 1, we have
54 * exp(x) = 2^k * exp(r)
55 *
56 * Special cases:
57 * exp(INF) is INF, exp(NaN) is NaN;
58 * exp(-INF) is 0, and
59 * for finite argument, only exp(0)=1 is exact.
60 *
61 * Accuracy:
62 * according to an error analysis, the error is always less than
63 * 1 ulp (unit in the last place).
64 *
65 * Misc. info.
66 * For IEEE double
67 * if x > 7.09782712893383973096e+02 then exp(x) overflow
68 * if x < -7.45133219101941108420e+02 then exp(x) underflow
69 *
70 * Constants:
71 * The hexadecimal values are the intended ones for the following
72 * constants. The decimal values may be used, provided that the
73 * compiler will convert from decimal to binary accurately enough
74 * to produce the hexadecimal values shown.
75 */
76
77 #include "fdlibm.h"
78 #include "math_config.h"
79 #if __OBSOLETE_MATH_DOUBLE
80
81 #ifdef _NEED_FLOAT64
82
83 static const __float64
84 one = _F_64(1.0),
85 halF[2] = {_F_64(0.5),_F_64(-0.5),},
86 huge = _F_64(1.0e+300),
87 twom1000 = _F_64(9.33263618503218878990e-302), /* 2**-1000=0x01700000,0*/
88 o_threshold = _F_64(7.09782712893383973096e+02), /* 0x40862E42, 0xFEFA39EF */
89 u_threshold = _F_64(-7.45133219101941108420e+02), /* 0xc0874910, 0xD52D3051 */
90 ln2HI[2] ={ _F_64(6.93147180369123816490e-01), /* 0x3fe62e42, 0xfee00000 */
91 _F_64(-6.93147180369123816490e-01),},/* 0xbfe62e42, 0xfee00000 */
92 ln2LO[2] ={ _F_64(1.90821492927058770002e-10), /* 0x3dea39ef, 0x35793c76 */
93 _F_64(-1.90821492927058770002e-10),},/* 0xbdea39ef, 0x35793c76 */
94 invln2 = _F_64(1.44269504088896338700e+00), /* 0x3ff71547, 0x652b82fe */
95 P1 = _F_64(1.66666666666666019037e-01), /* 0x3FC55555, 0x5555553E */
96 P2 = _F_64(-2.77777777770155933842e-03), /* 0xBF66C16C, 0x16BEBD93 */
97 P3 = _F_64(6.61375632143793436117e-05), /* 0x3F11566A, 0xAF25DE2C */
98 P4 = _F_64(-1.65339022054652515390e-06), /* 0xBEBBBD41, 0xC5D26BF1 */
99 P5 = _F_64(4.13813679705723846039e-08); /* 0x3E663769, 0x72BEA4D0 */
100
101 __float64
exp64(__float64 x)102 exp64(__float64 x) /* default IEEE double exp */
103 {
104 __float64 y, hi, lo, c, t;
105 __int32_t k = 0, xsb;
106 __uint32_t hx;
107
108 GET_HIGH_WORD(hx, x);
109 xsb = (hx >> 31) & 1; /* sign bit of x */
110 hx &= 0x7fffffff; /* high word of |x| */
111
112 /* filter out non-finite argument */
113 if (hx >= 0x40862E42) { /* if |x|>=709.78... */
114 if (hx >= 0x7ff00000) {
115 __uint32_t lx;
116 GET_LOW_WORD(lx, x);
117 if (((hx & 0xfffff) | lx) != 0)
118 return x + x; /* NaN */
119 else
120 return (xsb == 0) ? x : _F_64(0.0); /* exp(+-inf)={inf,0} */
121 }
122 if (x > o_threshold)
123 return __math_oflow(0); /* overflow */
124 if (x < u_threshold)
125 return __math_uflow(0); /* underflow */
126 }
127
128 /* argument reduction */
129 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
130 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
131 hi = x - ln2HI[xsb];
132 lo = ln2LO[xsb];
133 k = 1 - xsb - xsb;
134 } else {
135 k = invln2 * x + halF[xsb];
136 t = k;
137 hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */
138 lo = t * ln2LO[0];
139 }
140 x = hi - lo;
141 } else if (hx < 0x3df00000) { /* when |x|<2**-32 */
142 if (huge + x > one)
143 return one + x; /* trigger inexact */
144 }
145
146 /* x is now in primary range */
147 t = x * x;
148 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
149 if (k == 0)
150 return one - ((x * c) / (c - _F_64(2.0)) - x);
151 else
152 y = one - ((lo - (x * c) / (_F_64(2.0) - c)) - hi);
153 if (k >= -1021) {
154 __uint32_t hy;
155 GET_HIGH_WORD(hy, y);
156 SET_HIGH_WORD(y, hy + (k << 20)); /* add k to y's exponent */
157 return y;
158 } else {
159 __uint32_t hy;
160 GET_HIGH_WORD(hy, y);
161 SET_HIGH_WORD(y, hy + ((k + 1000) << 20)); /* add k to y's exponent */
162 return y * twom1000;
163 }
164 }
165
166 _MATH_ALIAS_d_d(exp)
167
168 #endif /* _NEED_FLOAT64 */
169 #else
170 #include "../common/exp.c"
171 #endif /* __OBSOLETE_MATH_DOUBLE */
172
