Define IPTOS_CLASS_* macros according to RFC 2474.
[glibc.git] / sysdeps / ieee754 / ldbl-128ibm / s_expm1l.c
blobf631eddf59059b947515b27b4bda198c9b50ed2c
1 /* expm1l.c
3 * Exponential function, minus 1
4 * 128-bit long double precision
8 * SYNOPSIS:
10 * long double x, y, expm1l();
12 * y = expm1l( x );
16 * DESCRIPTION:
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
23 * x k f
24 * e = 2 e.
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].
30 * ACCURACY:
32 * Relative error:
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
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.
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.
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 */
54 #include <errno.h>
55 #include "math.h"
56 #include "math_private.h"
57 #include <math_ldbl_opt.h>
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 */
63 static const long double
64 P0 = 2.943520915569954073888921213330863757240E8L,
65 P1 = -5.722847283900608941516165725053359168840E7L,
66 P2 = 8.944630806357575461578107295909719817253E6L,
67 P3 = -7.212432713558031519943281748462837065308E5L,
68 P4 = 4.578962475841642634225390068461943438441E4L,
69 P5 = -1.716772506388927649032068540558788106762E3L,
70 P6 = 4.401308817383362136048032038528753151144E1L,
71 P7 = -4.888737542888633647784737721812546636240E-1L,
72 Q0 = 1.766112549341972444333352727998584753865E9L,
73 Q1 = -7.848989743695296475743081255027098295771E8L,
74 Q2 = 1.615869009634292424463780387327037251069E8L,
75 Q3 = -2.019684072836541751428967854947019415698E7L,
76 Q4 = 1.682912729190313538934190635536631941751E6L,
77 Q5 = -9.615511549171441430850103489315371768998E4L,
78 Q6 = 3.697714952261803935521187272204485251835E3L,
79 Q7 = -8.802340681794263968892934703309274564037E1L,
80 /* Q8 = 1.000000000000000000000000000000000000000E0 */
81 /* C1 + C2 = ln 2 */
83 C1 = 6.93145751953125E-1L,
84 C2 = 1.428606820309417232121458176568075500134E-6L,
85 /* ln (2^16384 * (1 - 2^-113)) */
86 maxlog = 1.1356523406294143949491931077970764891253E4L,
87 /* ln 2^-114 */
88 minarg = -7.9018778583833765273564461846232128760607E1L, big = 2e307L;
91 long double
92 __expm1l (long double x)
94 long double px, qx, xx;
95 int32_t ix, sign;
96 ieee854_long_double_shape_type u;
97 int k;
99 /* Detect infinity and NaN. */
100 u.value = x;
101 ix = u.parts32.w0;
102 sign = ix & 0x80000000;
103 ix &= 0x7fffffff;
104 if (ix >= 0x7ff00000)
106 /* Infinity. */
107 if (((ix & 0xfffff) | u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
109 if (sign)
110 return -1.0L;
111 else
112 return x;
114 /* NaN. No invalid exception. */
115 return x;
118 /* expm1(+- 0) = +- 0. */
119 if ((ix == 0) && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
120 return x;
122 /* Overflow. */
123 if (x > maxlog)
125 __set_errno (ERANGE);
126 return (big * big);
129 /* Minimum value. */
130 if (x < minarg)
131 return (4.0/big - 1.0L);
133 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
134 xx = C1 + C2; /* ln 2. */
135 px = __floorl (0.5 + x / xx);
136 k = px;
137 /* remainder times ln 2 */
138 x -= px * C1;
139 x -= px * C2;
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;
146 qx = (((((((x
147 + Q7) * x
148 + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
150 xx = x * x;
151 qx = x + (0.5 * xx + xx * px / qx);
153 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
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. */
159 px = __ldexpl (1.0L, k);
160 x = px * qx + (px - 1.0);
161 return x;
163 libm_hidden_def (__expm1l)
164 long_double_symbol (libm, __expm1l, expm1l);