]>
Commit | Line | Data |
---|---|---|
dc30f461 | 1 | /* Single-precision floating point e^x. |
f7a9f785 | 2 | Copyright (C) 1997-2016 Free Software Foundation, Inc. |
dc30f461 UD |
3 | This file is part of the GNU C Library. |
4 | Contributed by Geoffrey Keating <geoffk@ozemail.com.au> | |
f7eac6eb | 5 | |
dc30f461 | 6 | The GNU C Library is free software; you can redistribute it and/or |
41bdb6e2 AJ |
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. | |
f7eac6eb | 10 | |
dc30f461 UD |
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 | |
41bdb6e2 | 14 | Lesser General Public License for more details. |
f7eac6eb | 15 | |
41bdb6e2 | 16 | You should have received a copy of the GNU Lesser General Public |
59ba27a6 PE |
17 | License along with the GNU C Library; if not, see |
18 | <http://www.gnu.org/licenses/>. */ | |
dc30f461 UD |
19 | |
20 | /* How this works: | |
21 | ||
22 | The input value, x, is written as | |
23 | ||
24 | x = n * ln(2) + t/512 + delta[t] + x; | |
25 | ||
26 | where: | |
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 | |
31 | ||
32 | Then e^x is approximated as | |
33 | ||
34 | e^x = 2^n ( e^(t/512 + delta[t]) | |
bcf01e6d UD |
35 | + ( e^(t/512 + delta[t]) |
36 | * ( p(x + delta[t] + n * ln(2)) - delta ) ) ) | |
dc30f461 UD |
37 | |
38 | where | |
39 | - p(x) is a polynomial approximating e(x)-1; | |
40 | - e^(t/512 + delta[t]) is obtained from a table. | |
41 | ||
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. | |
44 | ||
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). | |
48 | */ | |
dc30f461 UD |
49 | #include <float.h> |
50 | #include <ieee754.h> | |
51 | #include <math.h> | |
52 | #include <fenv.h> | |
53 | #include <inttypes.h> | |
54 | #include <math_private.h> | |
55 | ||
56 | extern const float __exp_deltatable[178]; | |
57 | extern const double __exp_atable[355] /* __attribute__((mode(DF))) */; | |
58 | ||
d9a8d0ab UD |
59 | static const float TWOM100 = 7.88860905e-31; |
60 | static const float TWO127 = 1.7014118346e+38; | |
dc30f461 UD |
61 | |
62 | float | |
63 | __ieee754_expf (float x) | |
f7eac6eb | 64 | { |
dc30f461 UD |
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)) | |
69 | { | |
66000494 UD |
70 | static const float THREEp42 = 13194139533312.0; |
71 | static const float THREEp22 = 12582912.0; | |
dc30f461 UD |
72 | /* 1/ln(2). */ |
73 | #undef M_1_LN2 | |
74 | static const float M_1_LN2 = 1.44269502163f; | |
75 | /* ln(2) */ | |
76 | #undef M_LN2 | |
77 | static const double M_LN2 = .6931471805599452862; | |
78 | ||
79 | int tval; | |
80 | double x22, t, result, dx; | |
81 | float n, delta; | |
82 | union ieee754_double ex2_u; | |
dc30f461 | 83 | |
eb92c487 RH |
84 | { |
85 | SET_RESTORE_ROUND_NOEXF (FE_TONEAREST); | |
f7eac6eb | 86 | |
eb92c487 RH |
87 | /* Calculate n. */ |
88 | n = x * M_1_LN2 + THREEp22; | |
89 | n -= THREEp22; | |
90 | dx = x - n*M_LN2; | |
dc30f461 | 91 | |
eb92c487 RH |
92 | /* Calculate t/512. */ |
93 | t = dx + THREEp42; | |
94 | t -= THREEp42; | |
95 | dx -= t; | |
dc30f461 | 96 | |
eb92c487 RH |
97 | /* Compute tval = t. */ |
98 | tval = (int) (t * 512.0); | |
66000494 | 99 | |
eb92c487 RH |
100 | if (t >= 0) |
101 | delta = - __exp_deltatable[tval]; | |
102 | else | |
103 | delta = __exp_deltatable[-tval]; | |
dc30f461 | 104 | |
eb92c487 RH |
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; | |
dc30f461 | 108 | |
eb92c487 RH |
109 | /* Approximate e^(dx+delta) - 1, using a second-degree polynomial, |
110 | with maximum error in [-2^-10-2^-28,2^-10+2^-28] | |
111 | less than 5e-11. */ | |
112 | x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta; | |
113 | } | |
dc30f461 UD |
114 | |
115 | /* Return result. */ | |
dc30f461 UD |
116 | result = x22 * ex2_u.d + ex2_u.d; |
117 | return (float) result; | |
118 | } | |
119 | /* Exceptional cases: */ | |
120 | else if (isless (x, himark)) | |
121 | { | |
d81f90cc | 122 | if (isinf (x)) |
dc30f461 UD |
123 | /* e^-inf == 0, with no error. */ |
124 | return 0; | |
125 | else | |
126 | /* Underflow */ | |
127 | return TWOM100 * TWOM100; | |
128 | } | |
129 | else | |
130 | /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ | |
131 | return TWO127*x; | |
f7eac6eb | 132 | } |
bcf01e6d | 133 | strong_alias (__ieee754_expf, __expf_finite) |