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