1 /* Single-precision floating point e^x.
2 Copyright (C) 1997-2014 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <http://www.gnu.org/licenses/>. */
22 The input value, x, is written as
24 x = n * ln(2) + t/512 + delta[t] + x;
27 - n is an integer, 127 >= n >= -150;
28 - t is an integer, 177 >= t >= -177
29 - delta is based on a table entry, delta[t] < 2^-28
30 - x is whatever is left, |x| < 2^-10
32 Then e^x is approximated as
34 e^x = 2^n ( e^(t/512 + delta[t])
35 + ( e^(t/512 + delta[t])
36 * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
39 - p(x) is a polynomial approximating e(x)-1;
40 - e^(t/512 + delta[t]) is obtained from a table.
42 The table used is the same one as for the double precision version;
43 since we have the table, we might as well use it.
45 It turns out to be faster to do calculations in double precision than
46 to perform an 'accurate table method' expf, because of the range reduction
47 overhead (compare exp2f).
54 #include <math_private.h>
56 extern const float __exp_deltatable
[178];
57 extern const double __exp_atable
[355] /* __attribute__((mode(DF))) */;
59 static const float TWOM100
= 7.88860905e-31;
60 static const float TWO127
= 1.7014118346e+38;
63 __ieee754_expf (float x
)
65 static const float himark
= 88.72283935546875;
66 static const float lomark
= -103.972084045410;
67 /* Check for usual case. */
68 if (isless (x
, himark
) && isgreater (x
, lomark
))
70 static const float THREEp42
= 13194139533312.0;
71 static const float THREEp22
= 12582912.0;
74 static const float M_1_LN2
= 1.44269502163f
;
77 static const double M_LN2
= .6931471805599452862;
80 double x22
, t
, result
, dx
;
82 union ieee754_double ex2_u
;
85 SET_RESTORE_ROUND_NOEXF (FE_TONEAREST
);
88 n
= x
* M_1_LN2
+ THREEp22
;
92 /* Calculate t/512. */
97 /* Compute tval = t. */
98 tval
= (int) (t
* 512.0);
101 delta
= - __exp_deltatable
[tval
];
103 delta
= __exp_deltatable
[-tval
];
105 /* Compute ex2 = 2^n e^(t/512+delta[t]). */
106 ex2_u
.d
= __exp_atable
[tval
+177];
107 ex2_u
.ieee
.exponent
+= (int) n
;
109 /* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
110 with maximum error in [-2^-10-2^-28,2^-10+2^-28]
112 x22
= (0.5000000496709180453 * dx
+ 1.0000001192102037084) * dx
+ delta
;
116 result
= x22
* ex2_u
.d
+ ex2_u
.d
;
117 return (float) result
;
119 /* Exceptional cases: */
120 else if (isless (x
, himark
))
123 /* e^-inf == 0, with no error. */
127 return TWOM100
* TWOM100
;
130 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */
133 strong_alias (__ieee754_expf
, __expf_finite
)