1 | /* expm1l.c |
2 | * |
3 | * Exponential function, minus 1 |
4 | * 128-bit long double precision |
5 | * |
6 | * |
7 | * |
8 | * SYNOPSIS: |
9 | * |
10 | * long double x, y, expm1l(); |
11 | * |
12 | * y = expm1l( x ); |
13 | * |
14 | * |
15 | * |
16 | * DESCRIPTION: |
17 | * |
18 | * Returns e (2.71828...) raised to the x power, minus one. |
19 | * |
20 | * Range reduction is accomplished by separating the argument |
21 | * into an integer k and fraction f such that |
22 | * |
23 | * x k f |
24 | * e = 2 e. |
25 | * |
26 | * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1 |
27 | * in the basic range [-0.5 ln 2, 0.5 ln 2]. |
28 | * |
29 | * |
30 | * ACCURACY: |
31 | * |
32 | * Relative error: |
33 | * arithmetic domain # trials peak rms |
34 | * IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35 |
35 | * |
36 | */ |
37 | |
38 | /* Copyright 2001 by Stephen L. Moshier |
39 | |
40 | This library is free software; you can redistribute it and/or |
41 | modify it under the terms of the GNU Lesser General Public |
42 | License as published by the Free Software Foundation; either |
43 | version 2.1 of the License, or (at your option) any later version. |
44 | |
45 | This library is distributed in the hope that it will be useful, |
46 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
47 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
48 | Lesser General Public License for more details. |
49 | |
50 | You should have received a copy of the GNU Lesser General Public |
51 | License along with this library; if not, see |
52 | <https://www.gnu.org/licenses/>. */ |
53 | |
54 | #include <errno.h> |
55 | #include <math.h> |
56 | #include <math_private.h> |
57 | #include <math_ldbl_opt.h> |
58 | |
59 | /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) |
60 | -.5 ln 2 < x < .5 ln 2 |
61 | Theoretical peak relative error = 8.1e-36 */ |
62 | |
63 | static const long double |
64 | P0 = 2.943520915569954073888921213330863757240E8L, |
65 | P1 = -5.722847283900608941516165725053359168840E7L, |
66 | P2 = 8.944630806357575461578107295909719817253E6L, |
67 | P3 = -7.212432713558031519943281748462837065308E5L, |
68 | P4 = 4.578962475841642634225390068461943438441E4L, |
69 | P5 = -1.716772506388927649032068540558788106762E3L, |
70 | P6 = 4.401308817383362136048032038528753151144E1L, |
71 | P7 = -4.888737542888633647784737721812546636240E-1L, |
72 | Q0 = 1.766112549341972444333352727998584753865E9L, |
73 | Q1 = -7.848989743695296475743081255027098295771E8L, |
74 | Q2 = 1.615869009634292424463780387327037251069E8L, |
75 | Q3 = -2.019684072836541751428967854947019415698E7L, |
76 | Q4 = 1.682912729190313538934190635536631941751E6L, |
77 | Q5 = -9.615511549171441430850103489315371768998E4L, |
78 | Q6 = 3.697714952261803935521187272204485251835E3L, |
79 | Q7 = -8.802340681794263968892934703309274564037E1L, |
80 | /* Q8 = 1.000000000000000000000000000000000000000E0 */ |
81 | /* C1 + C2 = ln 2 */ |
82 | |
83 | C1 = 6.93145751953125E-1L, |
84 | C2 = 1.428606820309417232121458176568075500134E-6L, |
85 | /* ln 2^-114 */ |
86 | minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e290L; |
87 | |
88 | |
89 | long double |
90 | __expm1l (long double x) |
91 | { |
92 | long double px, qx, xx; |
93 | int32_t ix, lx, sign; |
94 | int k; |
95 | double xhi; |
96 | |
97 | /* Detect infinity and NaN. */ |
98 | xhi = ldbl_high (x); |
99 | EXTRACT_WORDS (ix, lx, xhi); |
100 | sign = ix & 0x80000000; |
101 | ix &= 0x7fffffff; |
102 | if (!sign && ix >= 0x40600000) |
103 | return __expl (x: x); |
104 | if (ix >= 0x7ff00000) |
105 | { |
106 | /* Infinity (which must be negative infinity). */ |
107 | if (((ix - 0x7ff00000) | lx) == 0) |
108 | return -1.0L; |
109 | /* NaN. Invalid exception if signaling. */ |
110 | return x + x; |
111 | } |
112 | |
113 | /* expm1(+- 0) = +- 0. */ |
114 | if ((ix | lx) == 0) |
115 | return x; |
116 | |
117 | /* Minimum value. */ |
118 | if (x < minarg) |
119 | return (4.0/big - 1.0L); |
120 | |
121 | /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ |
122 | xx = C1 + C2; /* ln 2. */ |
123 | px = floorl (0.5 + x / xx); |
124 | k = px; |
125 | /* remainder times ln 2 */ |
126 | x -= px * C1; |
127 | x -= px * C2; |
128 | |
129 | /* Approximate exp(remainder ln 2). */ |
130 | px = (((((((P7 * x |
131 | + P6) * x |
132 | + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x; |
133 | |
134 | qx = (((((((x |
135 | + Q7) * x |
136 | + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0; |
137 | |
138 | xx = x * x; |
139 | qx = x + (0.5 * xx + xx * px / qx); |
140 | |
141 | /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). |
142 | |
143 | We have qx = exp(remainder ln 2) - 1, so |
144 | exp(x) - 1 = 2^k (qx + 1) - 1 |
145 | = 2^k qx + 2^k - 1. */ |
146 | |
147 | px = __ldexpl (x: 1.0L, exponent: k); |
148 | x = px * qx + (px - 1.0); |
149 | return x; |
150 | } |
151 | libm_hidden_def (__expm1l) |
152 | long_double_symbol (libm, __expm1l, expm1l); |
153 | |