PR 53456 Fix incorrect ChangeLog entry.
[official-gcc.git] / libquadmath / math / log1pq.c
bloba466dc8921eefca510eefedf8996e21c393a805e
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 "quadmath-imp.h"
58 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
59 * 1/sqrt(2) <= 1+x < sqrt(2)
60 * Theoretical peak relative error = 5.3e-37,
61 * relative peak error spread = 2.3e-14
63 static const __float128
64 P12 = 1.538612243596254322971797716843006400388E-6Q,
65 P11 = 4.998469661968096229986658302195402690910E-1Q,
66 P10 = 2.321125933898420063925789532045674660756E1Q,
67 P9 = 4.114517881637811823002128927449878962058E2Q,
68 P8 = 3.824952356185897735160588078446136783779E3Q,
69 P7 = 2.128857716871515081352991964243375186031E4Q,
70 P6 = 7.594356839258970405033155585486712125861E4Q,
71 P5 = 1.797628303815655343403735250238293741397E5Q,
72 P4 = 2.854829159639697837788887080758954924001E5Q,
73 P3 = 3.007007295140399532324943111654767187848E5Q,
74 P2 = 2.014652742082537582487669938141683759923E5Q,
75 P1 = 7.771154681358524243729929227226708890930E4Q,
76 P0 = 1.313572404063446165910279910527789794488E4Q,
77 /* Q12 = 1.000000000000000000000000000000000000000E0Q, */
78 Q11 = 4.839208193348159620282142911143429644326E1Q,
79 Q10 = 9.104928120962988414618126155557301584078E2Q,
80 Q9 = 9.147150349299596453976674231612674085381E3Q,
81 Q8 = 5.605842085972455027590989944010492125825E4Q,
82 Q7 = 2.248234257620569139969141618556349415120E5Q,
83 Q6 = 6.132189329546557743179177159925690841200E5Q,
84 Q5 = 1.158019977462989115839826904108208787040E6Q,
85 Q4 = 1.514882452993549494932585972882995548426E6Q,
86 Q3 = 1.347518538384329112529391120390701166528E6Q,
87 Q2 = 7.777690340007566932935753241556479363645E5Q,
88 Q1 = 2.626900195321832660448791748036714883242E5Q,
89 Q0 = 3.940717212190338497730839731583397586124E4Q;
91 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
92 * where z = 2(x-1)/(x+1)
93 * 1/sqrt(2) <= x < sqrt(2)
94 * Theoretical peak relative error = 1.1e-35,
95 * relative peak error spread 1.1e-9
97 static const __float128
98 R5 = -8.828896441624934385266096344596648080902E-1Q,
99 R4 = 8.057002716646055371965756206836056074715E1Q,
100 R3 = -2.024301798136027039250415126250455056397E3Q,
101 R2 = 2.048819892795278657810231591630928516206E4Q,
102 R1 = -8.977257995689735303686582344659576526998E4Q,
103 R0 = 1.418134209872192732479751274970992665513E5Q,
104 /* S6 = 1.000000000000000000000000000000000000000E0Q, */
105 S5 = -1.186359407982897997337150403816839480438E2Q,
106 S4 = 3.998526750980007367835804959888064681098E3Q,
107 S3 = -5.748542087379434595104154610899551484314E4Q,
108 S2 = 4.001557694070773974936904547424676279307E5Q,
109 S1 = -1.332535117259762928288745111081235577029E6Q,
110 S0 = 1.701761051846631278975701529965589676574E6Q;
112 /* C1 + C2 = ln 2 */
113 static const __float128 C1 = 6.93145751953125E-1Q;
114 static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q;
116 static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q;
117 static const __float128 zero = 0.0Q;
120 __float128
121 log1pq (__float128 xm1)
123 __float128 x, y, z, r, s;
124 ieee854_float128 u;
125 int32_t hx;
126 int e;
128 /* Test for NaN or infinity input. */
129 u.value = xm1;
130 hx = u.words32.w0;
131 if (hx >= 0x7fff0000)
132 return xm1;
134 /* log1p(+- 0) = +- 0. */
135 if (((hx & 0x7fffffff) == 0)
136 && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
137 return xm1;
139 x = xm1 + 1.0Q;
141 /* log1p(-1) = -inf */
142 if (x <= 0.0Q)
144 if (x == 0.0Q)
145 return (-1.0Q / (x - x));
146 else
147 return (zero / (x - x));
150 /* Separate mantissa from exponent. */
152 /* Use frexp used so that denormal numbers will be handled properly. */
153 x = frexpq (x, &e);
155 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
156 where z = 2(x-1)/x+1). */
157 if ((e > 2) || (e < -2))
159 if (x < sqrth)
160 { /* 2( 2x-1 )/( 2x+1 ) */
161 e -= 1;
162 z = x - 0.5Q;
163 y = 0.5Q * z + 0.5Q;
165 else
166 { /* 2 (x-1)/(x+1) */
167 z = x - 0.5Q;
168 z -= 0.5Q;
169 y = 0.5Q * x + 0.5Q;
171 x = z / y;
172 z = x * x;
173 r = ((((R5 * z
174 + R4) * z
175 + R3) * z
176 + R2) * z
177 + R1) * z
178 + R0;
179 s = (((((z
180 + S5) * z
181 + S4) * z
182 + S3) * z
183 + S2) * z
184 + S1) * z
185 + S0;
186 z = x * (z * r / s);
187 z = z + e * C2;
188 z = z + x;
189 z = z + e * C1;
190 return (z);
194 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
196 if (x < sqrth)
198 e -= 1;
199 if (e != 0)
200 x = 2.0Q * x - 1.0Q; /* 2x - 1 */
201 else
202 x = xm1;
204 else
206 if (e != 0)
207 x = x - 1.0Q;
208 else
209 x = xm1;
211 z = x * x;
212 r = (((((((((((P12 * x
213 + P11) * x
214 + P10) * x
215 + P9) * x
216 + P8) * x
217 + P7) * x
218 + P6) * x
219 + P5) * x
220 + P4) * x
221 + P3) * x
222 + P2) * x
223 + P1) * x
224 + P0;
225 s = (((((((((((x
226 + Q11) * x
227 + Q10) * x
228 + Q9) * x
229 + Q8) * x
230 + Q7) * x
231 + Q6) * x
232 + Q5) * x
233 + Q4) * x
234 + Q3) * x
235 + Q2) * x
236 + Q1) * x
237 + Q0;
238 y = x * (z * r / s);
239 y = y + e * C2;
240 z = y - 0.5Q * z;
241 z = z + x;
242 z = z + e * C1;
243 return (z);