ldbl-128 fmal compile fix
[glibc.git] / sysdeps / ieee754 / ldbl-128 / s_log1pl.c
blob5e9cf490550428d8112ccaf7fc1736b61985368a
1 /* log1pl.c
3 * Relative error logarithm
4 * Natural logarithm of 1+x, 128-bit long double precision
8 * SYNOPSIS:
10 * long double x, y, log1pl();
12 * y = log1pl( x );
16 * DESCRIPTION:
18 * Returns the base e (2.718...) logarithm of 1+x.
20 * The argument 1+x is separated into its exponent and fractional
21 * parts. If the exponent is between -1 and +1, the logarithm
22 * of the fraction is approximated by
24 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
26 * Otherwise, setting z = 2(w-1)/(w+1),
28 * log(w) = z + z^3 P(z)/Q(z).
32 * ACCURACY:
34 * Relative error:
35 * arithmetic domain # trials peak rms
36 * IEEE -1, 8 100000 1.9e-34 4.3e-35
39 /* Copyright 2001 by Stephen L. Moshier
41 This library is free software; you can redistribute it and/or
42 modify it under the terms of the GNU Lesser General Public
43 License as published by the Free Software Foundation; either
44 version 2.1 of the License, or (at your option) any later version.
46 This library is distributed in the hope that it will be useful,
47 but WITHOUT ANY WARRANTY; without even the implied warranty of
48 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
49 Lesser General Public License for more details.
51 You should have received a copy of the GNU Lesser General Public
52 License along with this library; if not, write to the Free Software
53 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
56 #include "math.h"
57 #include "math_private.h"
59 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
60 * 1/sqrt(2) <= 1+x < sqrt(2)
61 * Theoretical peak relative error = 5.3e-37,
62 * relative peak error spread = 2.3e-14
64 static const long double
65 P12 = 1.538612243596254322971797716843006400388E-6L,
66 P11 = 4.998469661968096229986658302195402690910E-1L,
67 P10 = 2.321125933898420063925789532045674660756E1L,
68 P9 = 4.114517881637811823002128927449878962058E2L,
69 P8 = 3.824952356185897735160588078446136783779E3L,
70 P7 = 2.128857716871515081352991964243375186031E4L,
71 P6 = 7.594356839258970405033155585486712125861E4L,
72 P5 = 1.797628303815655343403735250238293741397E5L,
73 P4 = 2.854829159639697837788887080758954924001E5L,
74 P3 = 3.007007295140399532324943111654767187848E5L,
75 P2 = 2.014652742082537582487669938141683759923E5L,
76 P1 = 7.771154681358524243729929227226708890930E4L,
77 P0 = 1.313572404063446165910279910527789794488E4L,
78 /* Q12 = 1.000000000000000000000000000000000000000E0L, */
79 Q11 = 4.839208193348159620282142911143429644326E1L,
80 Q10 = 9.104928120962988414618126155557301584078E2L,
81 Q9 = 9.147150349299596453976674231612674085381E3L,
82 Q8 = 5.605842085972455027590989944010492125825E4L,
83 Q7 = 2.248234257620569139969141618556349415120E5L,
84 Q6 = 6.132189329546557743179177159925690841200E5L,
85 Q5 = 1.158019977462989115839826904108208787040E6L,
86 Q4 = 1.514882452993549494932585972882995548426E6L,
87 Q3 = 1.347518538384329112529391120390701166528E6L,
88 Q2 = 7.777690340007566932935753241556479363645E5L,
89 Q1 = 2.626900195321832660448791748036714883242E5L,
90 Q0 = 3.940717212190338497730839731583397586124E4L;
92 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
93 * where z = 2(x-1)/(x+1)
94 * 1/sqrt(2) <= x < sqrt(2)
95 * Theoretical peak relative error = 1.1e-35,
96 * relative peak error spread 1.1e-9
98 static const long double
99 R5 = -8.828896441624934385266096344596648080902E-1L,
100 R4 = 8.057002716646055371965756206836056074715E1L,
101 R3 = -2.024301798136027039250415126250455056397E3L,
102 R2 = 2.048819892795278657810231591630928516206E4L,
103 R1 = -8.977257995689735303686582344659576526998E4L,
104 R0 = 1.418134209872192732479751274970992665513E5L,
105 /* S6 = 1.000000000000000000000000000000000000000E0L, */
106 S5 = -1.186359407982897997337150403816839480438E2L,
107 S4 = 3.998526750980007367835804959888064681098E3L,
108 S3 = -5.748542087379434595104154610899551484314E4L,
109 S2 = 4.001557694070773974936904547424676279307E5L,
110 S1 = -1.332535117259762928288745111081235577029E6L,
111 S0 = 1.701761051846631278975701529965589676574E6L;
113 /* C1 + C2 = ln 2 */
114 static const long double C1 = 6.93145751953125E-1L;
115 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
117 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
118 /* ln (2^16384 * (1 - 2^-113)) */
119 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
120 static const long double zero = 0.0L;
122 long double
123 __log1pl (long double xm1)
125 long double x, y, z, r, s;
126 ieee854_long_double_shape_type u;
127 int32_t hx;
128 int e;
130 /* Test for NaN or infinity input. */
131 u.value = xm1;
132 hx = u.parts32.w0;
133 if (hx >= 0x7fff0000)
134 return xm1;
136 /* log1p(+- 0) = +- 0. */
137 if (((hx & 0x7fffffff) == 0)
138 && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
139 return xm1;
141 x = xm1 + 1.0L;
143 /* log1p(-1) = -inf */
144 if (x <= 0.0L)
146 if (x == 0.0L)
147 return (-1.0L / (x - x));
148 else
149 return (zero / (x - x));
152 /* Separate mantissa from exponent. */
154 /* Use frexp used so that denormal numbers will be handled properly. */
155 x = __frexpl (x, &e);
157 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
158 where z = 2(x-1)/x+1). */
159 if ((e > 2) || (e < -2))
161 if (x < sqrth)
162 { /* 2( 2x-1 )/( 2x+1 ) */
163 e -= 1;
164 z = x - 0.5L;
165 y = 0.5L * z + 0.5L;
167 else
168 { /* 2 (x-1)/(x+1) */
169 z = x - 0.5L;
170 z -= 0.5L;
171 y = 0.5L * x + 0.5L;
173 x = z / y;
174 z = x * x;
175 r = ((((R5 * z
176 + R4) * z
177 + R3) * z
178 + R2) * z
179 + R1) * z
180 + R0;
181 s = (((((z
182 + S5) * z
183 + S4) * z
184 + S3) * z
185 + S2) * z
186 + S1) * z
187 + S0;
188 z = x * (z * r / s);
189 z = z + e * C2;
190 z = z + x;
191 z = z + e * C1;
192 return (z);
196 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
198 if (x < sqrth)
200 e -= 1;
201 if (e != 0)
202 x = 2.0L * x - 1.0L; /* 2x - 1 */
203 else
204 x = xm1;
206 else
208 if (e != 0)
209 x = x - 1.0L;
210 else
211 x = xm1;
213 z = x * x;
214 r = (((((((((((P12 * x
215 + P11) * x
216 + P10) * x
217 + P9) * x
218 + P8) * x
219 + P7) * x
220 + P6) * x
221 + P5) * x
222 + P4) * x
223 + P3) * x
224 + P2) * x
225 + P1) * x
226 + P0;
227 s = (((((((((((x
228 + Q11) * x
229 + Q10) * x
230 + Q9) * x
231 + Q8) * x
232 + Q7) * x
233 + Q6) * x
234 + Q5) * x
235 + Q4) * x
236 + Q3) * x
237 + Q2) * x
238 + Q1) * x
239 + Q0;
240 y = x * (z * r / s);
241 y = y + e * C2;
242 z = y - 0.5L * z;
243 z = z + x;
244 z = z + e * C1;
245 return (z);
248 weak_alias (__log1pl, log1pl)