1/* Quad-precision floating point e^x.
2 Copyright (C) 1999-2020 Free Software Foundation, Inc.
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
19 License along with the GNU C Library; if not, see
20 <https://www.gnu.org/licenses/>. */
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-barriers.h>
68#include <math_private.h>
69#include <math-underflow.h>
70#include <stdlib.h>
71#include "t_expl.h"
72#include <libm-alias-finite.h>
73
74static const _Float128 C[] = {
75/* Smallest integer x for which e^x overflows. */
76#define himark C[0]
77 L(11356.523406294143949491931077970765),
78
79/* Largest integer x for which e^x underflows. */
80#define lomark C[1]
81L(-11433.4627433362978788372438434526231),
82
83/* 3x2^96 */
84#define THREEp96 C[2]
85 L(59421121885698253195157962752.0),
86
87/* 3x2^103 */
88#define THREEp103 C[3]
89 L(30423614405477505635920876929024.0),
90
91/* 3x2^111 */
92#define THREEp111 C[4]
93 L(7788445287802241442795744493830144.0),
94
95/* 1/ln(2) */
96#define M_1_LN2 C[5]
97 L(1.44269504088896340735992468100189204),
98
99/* first 93 bits of ln(2) */
100#define M_LN2_0 C[6]
101 L(0.693147180559945309417232121457981864),
102
103/* ln2_0 - ln(2) */
104#define M_LN2_1 C[7]
105L(-1.94704509238074995158795957333327386E-31),
106
107/* very small number */
108#define TINY C[8]
109 L(1.0e-4900),
110
111/* 2^16383 */
112#define TWO16383 C[9]
113 L(5.94865747678615882542879663314003565E+4931),
114
115/* 256 */
116#define TWO8 C[10]
117 256,
118
119/* 32768 */
120#define TWO15 C[11]
121 32768,
122
123/* Chebyshev polynom coefficients for (exp(x)-1)/x */
124#define P1 C[12]
125#define P2 C[13]
126#define P3 C[14]
127#define P4 C[15]
128#define P5 C[16]
129#define P6 C[17]
130 L(0.5),
131 L(1.66666666666666666666666666666666683E-01),
132 L(4.16666666666666666666654902320001674E-02),
133 L(8.33333333333333333333314659767198461E-03),
134 L(1.38888888889899438565058018857254025E-03),
135 L(1.98412698413981650382436541785404286E-04),
136};
137
138_Float128
139__ieee754_expl (_Float128 x)
140{
141 /* Check for usual case. */
142 if (isless (x, himark) && isgreater (x, lomark))
143 {
144 int tval1, tval2, unsafe, n_i;
145 _Float128 x22, n, t, result, xl;
146 union ieee854_long_double ex2_u, scale_u;
147 fenv_t oldenv;
148
149 feholdexcept (&oldenv);
150#ifdef FE_TONEAREST
151 fesetround (FE_TONEAREST);
152#endif
153
154 /* Calculate n. */
155 n = x * M_1_LN2 + THREEp111;
156 n -= THREEp111;
157 x = x - n * M_LN2_0;
158 xl = n * M_LN2_1;
159
160 /* Calculate t/256. */
161 t = x + THREEp103;
162 t -= THREEp103;
163
164 /* Compute tval1 = t. */
165 tval1 = (int) (t * TWO8);
166
167 x -= __expl_table[T_EXPL_ARG1+2*tval1];
168 xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
169
170 /* Calculate t/32768. */
171 t = x + THREEp96;
172 t -= THREEp96;
173
174 /* Compute tval2 = t. */
175 tval2 = (int) (t * TWO15);
176
177 x -= __expl_table[T_EXPL_ARG2+2*tval2];
178 xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
179
180 x = x + xl;
181
182 /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */
183 ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
184 * __expl_table[T_EXPL_RES2 + tval2];
185 n_i = (int)n;
186 /* 'unsafe' is 1 iff n_1 != 0. */
187 unsafe = abs(n_i) >= 15000;
188 ex2_u.ieee.exponent += n_i >> unsafe;
189
190 /* Compute scale = 2^n_1. */
191 scale_u.d = 1;
192 scale_u.ieee.exponent += n_i - (n_i >> unsafe);
193
194 /* Approximate e^x2 - 1, using a seventh-degree polynomial,
195 with maximum error in [-2^-16-2^-53,2^-16+2^-53]
196 less than 4.8e-39. */
197 x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
198 math_force_eval (x22);
199
200 /* Return result. */
201 fesetenv (&oldenv);
202
203 result = x22 * ex2_u.d + ex2_u.d;
204
205 /* Now we can test whether the result is ultimate or if we are unsure.
206 In the later case we should probably call a mpn based routine to give
207 the ultimate result.
208 Empirically, this routine is already ultimate in about 99.9986% of
209 cases, the test below for the round to nearest case will be false
210 in ~ 99.9963% of cases.
211 Without proc2 routine maximum error which has been seen is
212 0.5000262 ulp.
213
214 union ieee854_long_double ex3_u;
215
216 #ifdef FE_TONEAREST
217 fesetround (FE_TONEAREST);
218 #endif
219 ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
220 ex2_u.d = result;
221 ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
222 - ex2_u.ieee.exponent;
223 n_i = abs (ex3_u.d);
224 n_i = (n_i + 1) / 2;
225 fesetenv (&oldenv);
226 #ifdef FE_TONEAREST
227 if (fegetround () == FE_TONEAREST)
228 n_i -= 0x4000;
229 #endif
230 if (!n_i) {
231 return __ieee754_expl_proc2 (origx);
232 }
233 */
234 if (!unsafe)
235 return result;
236 else
237 {
238 result *= scale_u.d;
239 math_check_force_underflow_nonneg (result);
240 return result;
241 }
242 }
243 /* Exceptional cases: */
244 else if (isless (x, himark))
245 {
246 if (isinf (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 TWO16383*x;
256}
257libm_alias_finite (__ieee754_expl, __expl)
258