1 /* Quad-precision floating point e^x.
2 Copyright (C) 1999,2004,2006, 2008 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>
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.
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.
18 You should have received a copy of the GNU Lesser General Public
19 License along with the GNU C Library; if not, write to the Free
20 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
23 /* The basic design here is from
24 Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with
25 Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991,
28 We work with number pairs where the first number is the high part and
29 the second one is the low part. Arithmetic with the high part numbers must
30 be exact, without any roundoff errors.
32 The input value, X, is written as
33 X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x
34 - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl
37 - n is an integer, 16384 >= n >= -16495;
38 - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205
39 - t1 is an integer, 89 >= t1 >= -89
40 - t2 is an integer, 65 >= t2 >= -65
41 - |arg1[t1]-t1/256.0| < 2^-53
42 - |arg2[t2]-t2/32768.0| < 2^-53
43 - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53
45 Then e^x is approximated as
47 e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
48 + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
49 * p (x + xl + n * ln(2)_1))
51 - p(x) is a polynomial approximating e(x)-1
52 - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table
53 - e^(arg2[t2]_0 + arg2[t2]_1) likewise
54 - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1.
56 If it happens that n_1 == 0 (this is the usual case), that multiplication
68 #include <math_private.h>
69 #include <sysdeps/ieee754/ldbl-128/t_expl.h>
71 static const long double C
[] = {
72 /* Smallest integer x for which e^x overflows. */
74 709.08956571282405153382846025171462914L,
76 /* Largest integer x for which e^x underflows. */
78 -709.08956571282405153382846025171462914L,
82 59421121885698253195157962752.0L,
85 #define THREEp103 C[3]
86 30423614405477505635920876929024.0L,
89 #define THREEp111 C[4]
90 7788445287802241442795744493830144.0L,
94 1.44269504088896340735992468100189204L,
96 /* first 93 bits of ln(2) */
98 0.693147180559945309417232121457981864L,
102 -1.94704509238074995158795957333327386E-31L,
104 /* very small number */
110 8.988465674311579538646525953945123668E+307L,
120 /* Chebyshev polynom coeficients for (exp(x)-1)/x */
128 1.66666666666666666666666666666666683E-01L,
129 4.16666666666666666666654902320001674E-02L,
130 8.33333333333333333333314659767198461E-03L,
131 1.38888888889899438565058018857254025E-03L,
132 1.98412698413981650382436541785404286E-04L,
136 __ieee754_expl (long double x
)
138 /* Check for usual case. */
139 if (isless (x
, himark
) && isgreater (x
, lomark
))
141 int tval1
, tval2
, unsafe
, n_i
, exponent2
;
142 long double x22
, n
, result
, xl
;
143 union ibm_extended_long_double ex2_u
, scale_u
;
146 feholdexcept (&oldenv
);
148 fesetround (FE_TONEAREST
);
151 n
= __roundl (x
*M_1_LN2
);
155 tval1
= __roundl (x
*TWO8
);
156 x
-= __expl_table
[T_EXPL_ARG1
+2*tval1
];
157 xl
-= __expl_table
[T_EXPL_ARG1
+2*tval1
+1];
159 tval2
= __roundl (x
*TWO15
);
160 x
-= __expl_table
[T_EXPL_ARG2
+2*tval2
];
161 xl
-= __expl_table
[T_EXPL_ARG2
+2*tval2
+1];
165 /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */
166 ex2_u
.d
= __expl_table
[T_EXPL_RES1
+ tval1
]
167 * __expl_table
[T_EXPL_RES2
+ tval2
];
169 /* 'unsafe' is 1 iff n_1 != 0. */
170 unsafe
= fabsl(n_i
) >= -LDBL_MIN_EXP
- 1;
171 ex2_u
.ieee
.exponent
+= n_i
>> unsafe
;
172 /* Fortunately, there are no subnormal lowpart doubles in
173 __expl_table, only normal values and zeros.
174 But after scaling it can be subnormal. */
175 exponent2
= ex2_u
.ieee
.exponent2
+ (n_i
>> unsafe
);
176 if (ex2_u
.ieee
.exponent2
== 0)
177 /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */;
178 else if (exponent2
> 0)
179 ex2_u
.ieee
.exponent2
= exponent2
;
180 else if (exponent2
<= -54)
182 ex2_u
.ieee
.exponent2
= 0;
183 ex2_u
.ieee
.mantissa2
= 0;
184 ex2_u
.ieee
.mantissa3
= 0;
189 two54
= 1.80143985094819840000e+16, /* 4350000000000000 */
190 twom54
= 5.55111512312578270212e-17; /* 3C90000000000000 */
191 ex2_u
.dd
[1] *= two54
;
192 ex2_u
.ieee
.exponent2
+= n_i
>> unsafe
;
193 ex2_u
.dd
[1] *= twom54
;
196 /* Compute scale = 2^n_1. */
198 scale_u
.ieee
.exponent
+= n_i
- (n_i
>> unsafe
);
200 /* Approximate e^x2 - 1, using a seventh-degree polynomial,
201 with maximum error in [-2^-16-2^-53,2^-16+2^-53]
202 less than 4.8e-39. */
203 x22
= x
+ x
*x
*(P1
+x
*(P2
+x
*(P3
+x
*(P4
+x
*(P5
+x
*P6
)))));
208 result
= x22
* ex2_u
.d
+ ex2_u
.d
;
210 /* Now we can test whether the result is ultimate or if we are unsure.
211 In the later case we should probably call a mpn based routine to give
213 Empirically, this routine is already ultimate in about 99.9986% of
214 cases, the test below for the round to nearest case will be false
215 in ~ 99.9963% of cases.
216 Without proc2 routine maximum error which has been seen is
219 union ieee854_long_double ex3_u;
222 fesetround (FE_TONEAREST);
224 ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
226 ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
227 - ex2_u.ieee.exponent;
232 if (fegetround () == FE_TONEAREST)
236 return __ieee754_expl_proc2 (origx);
242 return result
* scale_u
.d
;
244 /* Exceptional cases: */
245 else if (isless (x
, himark
))
248 /* e^-inf == 0, with no error. */
255 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */