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
55
56 #include <errno.h>
57 #include <float.h>
58 #include <math.h>
59 #include <math_private.h>
60 #include <math-underflow.h>
61 #include <libm-alias-ldouble.h>
62
63 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
64 -.5 ln 2 < x < .5 ln 2
65 Theoretical peak relative error = 8.1e-36 */
66
67 static const _Float128
68 P0 = L(2.943520915569954073888921213330863757240E8),
69 P1 = L(-5.722847283900608941516165725053359168840E7),
70 P2 = L(8.944630806357575461578107295909719817253E6),
71 P3 = L(-7.212432713558031519943281748462837065308E5),
72 P4 = L(4.578962475841642634225390068461943438441E4),
73 P5 = L(-1.716772506388927649032068540558788106762E3),
74 P6 = L(4.401308817383362136048032038528753151144E1),
75 P7 = L(-4.888737542888633647784737721812546636240E-1),
76 Q0 = L(1.766112549341972444333352727998584753865E9),
77 Q1 = L(-7.848989743695296475743081255027098295771E8),
78 Q2 = L(1.615869009634292424463780387327037251069E8),
79 Q3 = L(-2.019684072836541751428967854947019415698E7),
80 Q4 = L(1.682912729190313538934190635536631941751E6),
81 Q5 = L(-9.615511549171441430850103489315371768998E4),
82 Q6 = L(3.697714952261803935521187272204485251835E3),
83 Q7 = L(-8.802340681794263968892934703309274564037E1),
84 /* Q8 = 1.000000000000000000000000000000000000000E0 */
85 /* C1 + C2 = ln 2 */
86
87 C1 = L(6.93145751953125E-1),
88 C2 = L(1.428606820309417232121458176568075500134E-6),
89 /* ln 2^-114 */
90 minarg = L(-7.9018778583833765273564461846232128760607E1), big = L(1e4932);
91
92
93 _Float128
94 __expm1l (_Float128 x)
95 {
96 _Float128 px, qx, xx;
97 int32_t ix, sign;
98 ieee854_long_double_shape_type u;
99 int k;
100
101 /* Detect infinity and NaN. */
102 u.value = x;
103 ix = u.parts32.w0;
104 sign = ix & 0x80000000;
105 ix &= 0x7fffffff;
106 if (!sign && ix >= 0x40060000)
107 {
108 /* If num is positive and exp >= 6 use plain exp. */
109 return __expl (x);
110 }
111 if (ix >= 0x7fff0000)
112 {
113 /* Infinity (which must be negative infinity). */
114 if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
115 return -1;
116 /* NaN. Invalid exception if signaling. */
117 return x + x;
118 }
119
120 /* expm1(+- 0) = +- 0. */
121 if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
122 return x;
123
124 /* Minimum value. */
125 if (x < minarg)
126 return (4.0/big - 1);
127
128 /* Avoid internal underflow when result does not underflow, while
129 ensuring underflow (without returning a zero of the wrong sign)
130 when the result does underflow. */
131 if (fabsl (x) < L(0x1p-113))
132 {
133 math_check_force_underflow (x);
134 return x;
135 }
136
137 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
138 xx = C1 + C2; /* ln 2. */
139 px = floorl (0.5 + x / xx);
140 k = px;
141 /* remainder times ln 2 */
142 x -= px * C1;
143 x -= px * C2;
144
145 /* Approximate exp(remainder ln 2). */
146 px = (((((((P7 * x
147 + P6) * x
148 + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
149
150 qx = (((((((x
151 + Q7) * x
152 + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
153
154 xx = x * x;
155 qx = x + (0.5 * xx + xx * px / qx);
156
157 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
158
159 We have qx = exp(remainder ln 2) - 1, so
160 exp(x) - 1 = 2^k (qx + 1) - 1
161 = 2^k qx + 2^k - 1. */
162
163 px = __ldexpl (1, k);
164 x = px * qx + (px - 1.0);
165 return x;
166 }
167 libm_hidden_def (__expm1l)
168 libm_alias_ldouble (__expm1, expm1)