]>
Commit | Line | Data |
---|---|---|
f964490f | 1 | /* Quad-precision floating point e^x. |
bcf01e6d | 2 | Copyright (C) 1999,2004,2006, 2008, 2011 Free Software Foundation, Inc. |
f964490f RM |
3 | This file is part of the GNU C Library. |
4 | Contributed by Jakub Jelinek <jj@ultra.linux.cz> | |
5 | Partly based on double-precision code | |
6 | by Geoffrey Keating <geoffk@ozemail.com.au> | |
7 | ||
8 | The GNU C Library is free software; you can redistribute it and/or | |
9 | modify it under the terms of the GNU Lesser General Public | |
10 | License as published by the Free Software Foundation; either | |
11 | version 2.1 of the License, or (at your option) any later version. | |
12 | ||
13 | The GNU C Library is distributed in the hope that it will be useful, | |
14 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
15 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
16 | Lesser General Public License for more details. | |
17 | ||
18 | You should have received a copy of the GNU Lesser General Public | |
59ba27a6 PE |
19 | License along with the GNU C Library; if not, see |
20 | <http://www.gnu.org/licenses/>. */ | |
f964490f RM |
21 | |
22 | /* The basic design here is from | |
23 | Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with | |
24 | Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991, | |
25 | pp. 410-423. | |
26 | ||
27 | We work with number pairs where the first number is the high part and | |
28 | the second one is the low part. Arithmetic with the high part numbers must | |
29 | be exact, without any roundoff errors. | |
30 | ||
31 | The input value, X, is written as | |
32 | X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x | |
33 | - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl | |
34 | ||
35 | where: | |
36 | - n is an integer, 16384 >= n >= -16495; | |
37 | - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205 | |
38 | - t1 is an integer, 89 >= t1 >= -89 | |
39 | - t2 is an integer, 65 >= t2 >= -65 | |
40 | - |arg1[t1]-t1/256.0| < 2^-53 | |
41 | - |arg2[t2]-t2/32768.0| < 2^-53 | |
42 | - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53 | |
43 | ||
44 | Then e^x is approximated as | |
45 | ||
46 | e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1) | |
47 | + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1) | |
48 | * p (x + xl + n * ln(2)_1)) | |
49 | where: | |
50 | - p(x) is a polynomial approximating e(x)-1 | |
51 | - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table | |
52 | - e^(arg2[t2]_0 + arg2[t2]_1) likewise | |
53 | - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1. | |
54 | ||
55 | If it happens that n_1 == 0 (this is the usual case), that multiplication | |
56 | is omitted. | |
57 | */ | |
58 | ||
59 | #ifndef _GNU_SOURCE | |
60 | #define _GNU_SOURCE | |
61 | #endif | |
62 | #include <float.h> | |
63 | #include <ieee754.h> | |
64 | #include <math.h> | |
65 | #include <fenv.h> | |
66 | #include <inttypes.h> | |
67 | #include <math_private.h> | |
68 | #include <sysdeps/ieee754/ldbl-128/t_expl.h> | |
69 | ||
70 | static const long double C[] = { | |
71 | /* Smallest integer x for which e^x overflows. */ | |
72 | #define himark C[0] | |
73 | 709.08956571282405153382846025171462914L, | |
74 | ||
75 | /* Largest integer x for which e^x underflows. */ | |
76 | #define lomark C[1] | |
77 | -709.08956571282405153382846025171462914L, | |
78 | ||
79 | /* 3x2^96 */ | |
80 | #define THREEp96 C[2] | |
81 | 59421121885698253195157962752.0L, | |
82 | ||
83 | /* 3x2^103 */ | |
84 | #define THREEp103 C[3] | |
85 | 30423614405477505635920876929024.0L, | |
86 | ||
87 | /* 3x2^111 */ | |
88 | #define THREEp111 C[4] | |
89 | 7788445287802241442795744493830144.0L, | |
90 | ||
91 | /* 1/ln(2) */ | |
92 | #define M_1_LN2 C[5] | |
93 | 1.44269504088896340735992468100189204L, | |
94 | ||
95 | /* first 93 bits of ln(2) */ | |
96 | #define M_LN2_0 C[6] | |
97 | 0.693147180559945309417232121457981864L, | |
98 | ||
99 | /* ln2_0 - ln(2) */ | |
100 | #define M_LN2_1 C[7] | |
101 | -1.94704509238074995158795957333327386E-31L, | |
102 | ||
103 | /* very small number */ | |
104 | #define TINY C[8] | |
105 | 1.0e-308L, | |
106 | ||
107 | /* 2^16383 */ | |
108 | #define TWO1023 C[9] | |
109 | 8.988465674311579538646525953945123668E+307L, | |
110 | ||
111 | /* 256 */ | |
112 | #define TWO8 C[10] | |
113 | 256.0L, | |
114 | ||
115 | /* 32768 */ | |
116 | #define TWO15 C[11] | |
117 | 32768.0L, | |
118 | ||
119 | /* Chebyshev polynom coeficients for (exp(x)-1)/x */ | |
120 | #define P1 C[12] | |
121 | #define P2 C[13] | |
122 | #define P3 C[14] | |
123 | #define P4 C[15] | |
124 | #define P5 C[16] | |
125 | #define P6 C[17] | |
126 | 0.5L, | |
127 | 1.66666666666666666666666666666666683E-01L, | |
128 | 4.16666666666666666666654902320001674E-02L, | |
129 | 8.33333333333333333333314659767198461E-03L, | |
130 | 1.38888888889899438565058018857254025E-03L, | |
131 | 1.98412698413981650382436541785404286E-04L, | |
132 | }; | |
133 | ||
134 | long double | |
135 | __ieee754_expl (long double x) | |
136 | { | |
137 | /* Check for usual case. */ | |
138 | if (isless (x, himark) && isgreater (x, lomark)) | |
139 | { | |
140 | int tval1, tval2, unsafe, n_i, exponent2; | |
141 | long double x22, n, result, xl; | |
142 | union ibm_extended_long_double ex2_u, scale_u; | |
143 | fenv_t oldenv; | |
144 | ||
145 | feholdexcept (&oldenv); | |
146 | #ifdef FE_TONEAREST | |
147 | fesetround (FE_TONEAREST); | |
148 | #endif | |
149 | ||
246ec411 | 150 | n = __roundl (x*M_1_LN2); |
f964490f RM |
151 | x = x-n*M_LN2_0; |
152 | xl = n*M_LN2_1; | |
153 | ||
246ec411 | 154 | tval1 = __roundl (x*TWO8); |
f964490f RM |
155 | x -= __expl_table[T_EXPL_ARG1+2*tval1]; |
156 | xl -= __expl_table[T_EXPL_ARG1+2*tval1+1]; | |
157 | ||
246ec411 | 158 | tval2 = __roundl (x*TWO15); |
f964490f RM |
159 | x -= __expl_table[T_EXPL_ARG2+2*tval2]; |
160 | xl -= __expl_table[T_EXPL_ARG2+2*tval2+1]; | |
161 | ||
162 | x = x + xl; | |
163 | ||
164 | /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */ | |
165 | ex2_u.d = __expl_table[T_EXPL_RES1 + tval1] | |
166 | * __expl_table[T_EXPL_RES2 + tval2]; | |
167 | n_i = (int)n; | |
168 | /* 'unsafe' is 1 iff n_1 != 0. */ | |
169 | unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1; | |
170 | ex2_u.ieee.exponent += n_i >> unsafe; | |
171 | /* Fortunately, there are no subnormal lowpart doubles in | |
172 | __expl_table, only normal values and zeros. | |
173 | But after scaling it can be subnormal. */ | |
174 | exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe); | |
175 | if (ex2_u.ieee.exponent2 == 0) | |
176 | /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */; | |
177 | else if (exponent2 > 0) | |
178 | ex2_u.ieee.exponent2 = exponent2; | |
179 | else if (exponent2 <= -54) | |
180 | { | |
181 | ex2_u.ieee.exponent2 = 0; | |
182 | ex2_u.ieee.mantissa2 = 0; | |
183 | ex2_u.ieee.mantissa3 = 0; | |
184 | } | |
185 | else | |
186 | { | |
187 | static const double | |
188 | two54 = 1.80143985094819840000e+16, /* 4350000000000000 */ | |
189 | twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */ | |
190 | ex2_u.dd[1] *= two54; | |
191 | ex2_u.ieee.exponent2 += n_i >> unsafe; | |
192 | ex2_u.dd[1] *= twom54; | |
193 | } | |
194 | ||
195 | /* Compute scale = 2^n_1. */ | |
196 | scale_u.d = 1.0L; | |
197 | scale_u.ieee.exponent += n_i - (n_i >> unsafe); | |
198 | ||
199 | /* Approximate e^x2 - 1, using a seventh-degree polynomial, | |
200 | with maximum error in [-2^-16-2^-53,2^-16+2^-53] | |
201 | less than 4.8e-39. */ | |
202 | x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6))))); | |
203 | ||
204 | /* Return result. */ | |
205 | fesetenv (&oldenv); | |
206 | ||
207 | result = x22 * ex2_u.d + ex2_u.d; | |
208 | ||
209 | /* Now we can test whether the result is ultimate or if we are unsure. | |
210 | In the later case we should probably call a mpn based routine to give | |
211 | the ultimate result. | |
212 | Empirically, this routine is already ultimate in about 99.9986% of | |
213 | cases, the test below for the round to nearest case will be false | |
214 | in ~ 99.9963% of cases. | |
215 | Without proc2 routine maximum error which has been seen is | |
216 | 0.5000262 ulp. | |
217 | ||
218 | union ieee854_long_double ex3_u; | |
219 | ||
220 | #ifdef FE_TONEAREST | |
221 | fesetround (FE_TONEAREST); | |
222 | #endif | |
223 | ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d; | |
224 | ex2_u.d = result; | |
225 | ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS | |
226 | - ex2_u.ieee.exponent; | |
227 | n_i = abs (ex3_u.d); | |
228 | n_i = (n_i + 1) / 2; | |
229 | fesetenv (&oldenv); | |
230 | #ifdef FE_TONEAREST | |
231 | if (fegetround () == FE_TONEAREST) | |
232 | n_i -= 0x4000; | |
233 | #endif | |
234 | if (!n_i) { | |
235 | return __ieee754_expl_proc2 (origx); | |
236 | } | |
237 | */ | |
238 | if (!unsafe) | |
239 | return result; | |
240 | else | |
241 | return result * scale_u.d; | |
242 | } | |
243 | /* Exceptional cases: */ | |
244 | else if (isless (x, himark)) | |
245 | { | |
246 | if (__isinfl (x)) | |
247 | /* e^-inf == 0, with no error. */ | |
248 | return 0; | |
249 | else | |
250 | /* Underflow */ | |
251 | return TINY * TINY; | |
252 | } | |
253 | else | |
254 | /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ | |
255 | return TWO1023*x; | |
256 | } | |
bcf01e6d | 257 | strong_alias (__ieee754_expl, __expl_finite) |