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
63 static 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^16384 * (1 - 2^-113)) */
86 maxlog = 1.1356523406294143949491931077970764891253E4Q,
87 /* ln 2^-114 */
88 minarg = -7.9018778583833765273564461846232128760607E1Q;
89
90
91 __float128
expm1q(__float128 x)92 expm1q (__float128 x)
93 {
94 __float128 px, qx, xx;
95 int32_t ix, sign;
96 ieee854_float128 u;
97 int k;
98
99 /* Detect infinity and NaN. */
100 u.value = x;
101 ix = u.words32.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 expq (x);
108 }
109 if (ix >= 0x7fff0000)
110 {
111 /* Infinity. */
112 if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
113 {
114 if (sign)
115 return -1.0Q;
116 else
117 return x;
118 }
119 /* NaN. No invalid exception. */
120 return x;
121 }
122
123 /* expm1(+- 0) = +- 0. */
124 if ((ix == 0) && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
125 return x;
126
127 /* Overflow. */
128 if (x > maxlog)
129 {
130 errno = ERANGE;
131 return (HUGE_VALQ * HUGE_VALQ);
132 }
133
134 /* Minimum value. */
135 if (x < minarg)
136 return (4.0/HUGE_VALQ - 1.0Q);
137
138 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
139 xx = C1 + C2; /* ln 2. */
140 px = floorq (0.5 + x / xx);
141 k = px;
142 /* remainder times ln 2 */
143 x -= px * C1;
144 x -= px * C2;
145
146 /* Approximate exp(remainder ln 2). */
147 px = (((((((P7 * x
148 + P6) * x
149 + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
150
151 qx = (((((((x
152 + Q7) * x
153 + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
154
155 xx = x * x;
156 qx = x + (0.5 * xx + xx * px / qx);
157
158 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
159
160 We have qx = exp(remainder ln 2) - 1, so
161 exp(x) - 1 = 2^k (qx + 1) - 1
162 = 2^k qx + 2^k - 1. */
163
164 px = ldexpq (1.0Q, k);
165 x = px * qx + (px - 1.0);
166 return x;
167 }
168