1 /* Copyright JS Foundation and other contributors, http://js.foundation
2 *
3 * Licensed under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License.
5 * You may obtain a copy of the License at
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 *
15 * This file is based on work under the following copyright and permission
16 * notice:
17 *
18 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
19 *
20 * Permission to use, copy, modify, and distribute this
21 * software is freely granted, provided that this notice
22 * is preserved.
23 *
24 * @(#)e_exp.c 1.6 04/04/22
25 */
26
27 #include "jerry-libm-internal.h"
28
29 /* exp(x)
30 * Returns the exponential of x.
31 *
32 * Method:
33 * 1. Argument reduction:
34 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
35 * Given x, find r and integer k such that
36 *
37 * x = k*ln2 + r, |r| <= 0.5*ln2.
38 *
39 * Here r will be represented as r = hi-lo for better
40 * accuracy.
41 *
42 * 2. Approximation of exp(r) by a special rational function on
43 * the interval [0,0.34658]:
44 * Write
45 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
46 * We use a special Remes algorithm on [0,0.34658] to generate
47 * a polynomial of degree 5 to approximate R. The maximum error
48 * of this polynomial approximation is bounded by 2**-59. In
49 * other words,
50 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
51 * (where z=r*r, and the values of P1 to P5 are listed below)
52 * and
53 * | 5 | -59
54 * | 2.0+P1*z+...+P5*z - R(z) | <= 2
55 * | |
56 * The computation of exp(r) thus becomes
57 * 2*r
58 * exp(r) = 1 + -------
59 * R - r
60 * r*R1(r)
61 * = 1 + r + ----------- (for better accuracy)
62 * 2 - R1(r)
63 * where
64 * 2 4 10
65 * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
66 *
67 * 3. Scale back to obtain exp(x):
68 * From step 1, we have
69 * exp(x) = 2^k * exp(r)
70 *
71 * Special cases:
72 * exp(INF) is INF, exp(NaN) is NaN;
73 * exp(-INF) is 0, and
74 * for finite argument, only exp(0)=1 is exact.
75 *
76 * Accuracy:
77 * according to an error analysis, the error is always less than
78 * 1 ulp (unit in the last place).
79 *
80 * Misc. info:
81 * For IEEE double
82 * if x > 7.09782712893383973096e+02 then exp(x) overflow
83 * if x < -7.45133219101941108420e+02 then exp(x) underflow
84 *
85 * Constants:
86 * The hexadecimal values are the intended ones for the following
87 * constants. The decimal values may be used, provided that the
88 * compiler will convert from decimal to binary accurately enough
89 * to produce the hexadecimal values shown.
90 */
91
92 static const double halF[2] =
93 {
94 0.5,
95 -0.5,
96 };
97 static const double ln2HI[2] =
98 {
99 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
100 -6.93147180369123816490e-01, /* 0xbfe62e42, 0xfee00000 */
101 };
102 static const double ln2LO[2] =
103 {
104 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
105 -1.90821492927058770002e-10, /* 0xbdea39ef, 0x35793c76 */
106 };
107
108 #define one 1.0
109 #define huge 1.0e+300
110 #define twom1000 9.33263618503218878990e-302 /* 2**-1000=0x01700000,0 */
111 #define o_threshold 7.09782712893383973096e+02 /* 0x40862E42, 0xFEFA39EF */
112 #define u_threshold -7.45133219101941108420e+02 /* 0xc0874910, 0xD52D3051 */
113 #define invln2 1.44269504088896338700e+00 /* 0x3ff71547, 0x652b82fe */
114 #define P1 1.66666666666666019037e-01 /* 0x3FC55555, 0x5555553E */
115 #define P2 -2.77777777770155933842e-03 /* 0xBF66C16C, 0x16BEBD93 */
116 #define P3 6.61375632143793436117e-05 /* 0x3F11566A, 0xAF25DE2C */
117 #define P4 -1.65339022054652515390e-06 /* 0xBEBBBD41, 0xC5D26BF1 */
118 #define P5 4.13813679705723846039e-08 /* 0x3E663769, 0x72BEA4D0 */
119
120 double
exp(double x)121 exp (double x) /* default IEEE double exp */
122 {
123 double hi, lo, c, t;
124 int k = 0, xsb;
125 unsigned hx;
126
127 hx = __HI (x); /* high word of x */
128 xsb = (hx >> 31) & 1; /* sign bit of x */
129 hx &= 0x7fffffff; /* high word of |x| */
130
131 /* filter out non-finite argument */
132 if (hx >= 0x40862E42) /* if |x| >= 709.78... */
133 {
134 if (hx >= 0x7ff00000)
135 {
136 if (((hx & 0xfffff) | __LO (x)) != 0) /* NaN */
137 {
138 return x + x;
139 }
140 else /* exp(+-inf) = {inf,0} */
141 {
142 return (xsb == 0) ? x : 0.0;
143 }
144 }
145 if (x > o_threshold) /* overflow */
146 {
147 return huge * huge;
148 }
149 if (x < u_threshold) /* underflow */
150 {
151 return twom1000 * twom1000;
152 }
153 }
154
155 /* argument reduction */
156 if (hx > 0x3fd62e42) /* if |x| > 0.5 ln2 */
157 {
158 if (hx < 0x3FF0A2B2) /* and |x| < 1.5 ln2 */
159 {
160 hi = x - ln2HI[xsb];
161 lo = ln2LO[xsb];
162 k = 1 - xsb - xsb;
163 }
164 else
165 {
166 k = (int) (invln2 * x + halF[xsb]);
167 t = k;
168 hi = x - t * ln2HI[0]; /* t * ln2HI is exact here */
169 lo = t * ln2LO[0];
170 }
171 x = hi - lo;
172 }
173 else if (hx < 0x3e300000) /* when |x| < 2**-28 */
174 {
175 if (huge + x > one) /* trigger inexact */
176 {
177 return one + x;
178 }
179 }
180 else
181 {
182 k = 0;
183 }
184
185 double_accessor ret;
186
187 /* x is now in primary range */
188 t = x * x;
189 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
190 if (k == 0)
191 {
192 return one - ((x * c) / (c - 2.0) - x);
193 }
194 else
195 {
196 ret.dbl = one - ((lo - (x * c) / (2.0 - c)) - hi);
197 }
198 if (k >= -1021)
199 {
200 ret.as_int.hi += (((unsigned int) k) << 20); /* add k to y's exponent */
201 return ret.dbl;
202 }
203 else
204 {
205 ret.as_int.hi += ((k + 1000) << 20); /* add k to y's exponent */
206 return ret.dbl * twom1000;
207 }
208 } /* exp */
209
210 #undef one
211 #undef huge
212 #undef twom1000
213 #undef o_threshold
214 #undef u_threshold
215 #undef invln2
216 #undef P1
217 #undef P2
218 #undef P3
219 #undef P4
220 #undef P5
221