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