3 * Exponential function, minus 1
4 * 128-bit long double precision
10 * long double x, y, expm1l();
18 * Returns e (2.71828...) raised to the x power, minus one.
20 * Range reduction is accomplished by separating the argument
21 * into an integer k and fraction f such that
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].
33 * arithmetic domain # trials peak rms
34 * IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35
38 /* Copyright 2001 by Stephen L. Moshier */
42 #include "math_private.h"
44 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
45 -.5 ln 2 < x < .5 ln 2
46 Theoretical peak relative error = 8.1e-36 */
49 P0
= 2.943520915569954073888921213330863757240E8L
,
50 P1
= -5.722847283900608941516165725053359168840E7L
,
51 P2
= 8.944630806357575461578107295909719817253E6L
,
52 P3
= -7.212432713558031519943281748462837065308E5L
,
53 P4
= 4.578962475841642634225390068461943438441E4L
,
54 P5
= -1.716772506388927649032068540558788106762E3L
,
55 P6
= 4.401308817383362136048032038528753151144E1L
,
56 P7
= -4.888737542888633647784737721812546636240E-1L,
57 Q0
= 1.766112549341972444333352727998584753865E9L
,
58 Q1
= -7.848989743695296475743081255027098295771E8L
,
59 Q2
= 1.615869009634292424463780387327037251069E8L
,
60 Q3
= -2.019684072836541751428967854947019415698E7L
,
61 Q4
= 1.682912729190313538934190635536631941751E6L
,
62 Q5
= -9.615511549171441430850103489315371768998E4L
,
63 Q6
= 3.697714952261803935521187272204485251835E3L
,
64 Q7
= -8.802340681794263968892934703309274564037E1L
,
65 /* Q8 = 1.000000000000000000000000000000000000000E0 */
68 C1
= 6.93145751953125E-1L,
69 C2
= 1.428606820309417232121458176568075500134E-6L,
70 /* ln (2^16384 * (1 - 2^-113)) */
71 maxlog
= 1.1356523406294143949491931077970764891253E4L
,
73 minarg
= -7.9018778583833765273564461846232128760607E1L
, big
= 2e4932L
;
77 __expm1l (long double x
)
79 long double px
, qx
, xx
;
81 ieee854_long_double_shape_type u
;
84 /* Detect infinity and NaN. */
87 sign
= ix
& 0x80000000;
92 if (((ix
& 0xffff) | u
.parts32
.w1
| u
.parts32
.w2
| u
.parts32
.w3
) == 0)
99 /* NaN. No invalid exception. */
103 /* expm1(+- 0) = +- 0. */
104 if ((ix
== 0) && (u
.parts32
.w1
| u
.parts32
.w2
| u
.parts32
.w3
) == 0)
113 return (4.0/big
- 1.0L);
115 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
116 xx
= C1
+ C2
; /* ln 2. */
117 px
= __floorl (0.5 + x
/ xx
);
119 /* remainder times ln 2 */
123 /* Approximate exp(remainder ln 2). */
126 + P5
) * x
+ P4
) * x
+ P3
) * x
+ P2
) * x
+ P1
) * x
+ P0
) * x
;
130 + Q6
) * x
+ Q5
) * x
+ Q4
) * x
+ Q3
) * x
+ Q2
) * x
+ Q1
) * x
+ Q0
;
133 qx
= x
+ (0.5 * xx
+ xx
* px
/ qx
);
135 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
137 We have qx = exp(remainder ln 2) - 1, so
138 exp(x) - 1 = 2^k (qx + 1) - 1
139 = 2^k qx + 2^k - 1. */
141 px
= ldexpl (1.0L, k
);
142 x
= px
* qx
+ (px
- 1.0);
146 weak_alias (__expm1l
, expm1l
)
147 #ifdef NO_LONG_DOUBLE
148 strong_alias (__expm1
, __expm1l
) weak_alias (__expm1
, expm1l
)