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