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