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
__expm1l(long double x)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);
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 (1.0L, 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