1/* expm1l.c
2 *
3 * Exponential function, minus 1
4 * 128-bit __float128 precision
5 *
6 *
7 *
8 * SYNOPSIS:
9 *
10 * __float128 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, write to the Free Software
52 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
53
54
55
56#include <errno.h>
57#include "quadmath-imp.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
63static const __float128
64 P0 = 2.943520915569954073888921213330863757240E8Q,
65 P1 = -5.722847283900608941516165725053359168840E7Q,
66 P2 = 8.944630806357575461578107295909719817253E6Q,
67 P3 = -7.212432713558031519943281748462837065308E5Q,
68 P4 = 4.578962475841642634225390068461943438441E4Q,
69 P5 = -1.716772506388927649032068540558788106762E3Q,
70 P6 = 4.401308817383362136048032038528753151144E1Q,
71 P7 = -4.888737542888633647784737721812546636240E-1Q,
72 Q0 = 1.766112549341972444333352727998584753865E9Q,
73 Q1 = -7.848989743695296475743081255027098295771E8Q,
74 Q2 = 1.615869009634292424463780387327037251069E8Q,
75 Q3 = -2.019684072836541751428967854947019415698E7Q,
76 Q4 = 1.682912729190313538934190635536631941751E6Q,
77 Q5 = -9.615511549171441430850103489315371768998E4Q,
78 Q6 = 3.697714952261803935521187272204485251835E3Q,
79 Q7 = -8.802340681794263968892934703309274564037E1Q,
80 /* Q8 = 1.000000000000000000000000000000000000000E0 */
81/* C1 + C2 = ln 2 */
82
83 C1 = 6.93145751953125E-1Q,
84 C2 = 1.428606820309417232121458176568075500134E-6Q,
85/* ln 2^-114 */
86 minarg = -7.9018778583833765273564461846232128760607E1Q;
87
88
89__float128
90expm1q (__float128 x)
91{
92 __float128 px, qx, xx;
93 int32_t ix, sign;
94 ieee854_float128 u;
95 int k;
96
97 /* Detect infinity and NaN. */
98 u.value = x;
99 ix = u.words32.w0;
100 sign = ix & 0x80000000;
101 ix &= 0x7fffffff;
102 if (!sign && ix >= 0x40060000)
103 {
104 /* If num is positive and exp >= 6 use plain exp. */
105 return expq (x);
106 }
107 if (ix >= 0x7fff0000)
108 {
109 /* Infinity (which must be negative infinity). */
110 if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
111 return -1.0Q;
112 /* NaN. Invalid exception if signaling. */
113 return x + x;
114 }
115
116 /* expm1(+- 0) = +- 0. */
117 if ((ix == 0) && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
118 return x;
119
120 /* Minimum value. */
121 if (x < minarg)
122 return (4.0/HUGE_VALQ - 1.0Q);
123
124 /* Avoid internal underflow when result does not underflow, while
125 ensuring underflow (without returning a zero of the wrong sign)
126 when the result does underflow. */
127 if (fabsq (x) < 0x1p-113Q)
128 {
129 math_check_force_underflow (x);
130 return x;
131 }
132
133 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
134 xx = C1 + C2; /* ln 2. */
135 px = floorq (0.5 + x / xx);
136 k = px;
137 /* remainder times ln 2 */
138 x -= px * C1;
139 x -= px * C2;
140
141 /* Approximate exp(remainder ln 2). */
142 px = (((((((P7 * x
143 + P6) * x
144 + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
145
146 qx = (((((((x
147 + Q7) * x
148 + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
149
150 xx = x * x;
151 qx = x + (0.5 * xx + xx * px / qx);
152
153 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
154
155 We have qx = exp(remainder ln 2) - 1, so
156 exp(x) - 1 = 2^k (qx + 1) - 1
157 = 2^k qx + 2^k - 1. */
158
159 px = ldexpq (1.0Q, k);
160 x = px * qx + (px - 1.0);
161 return x;
162}
163