Fix compilation failure with C++98 compilers
[official-gcc.git] / libquadmath / math / log10q.c
blob3afb1121267765bf6716ca206170e403775cdcc8
1 /* log10q.c
3 * Common logarithm, 128-bit __float128 precision
7 * SYNOPSIS:
9 * __float128 x, y, log10l();
11 * y = log10q( x );
15 * DESCRIPTION:
17 * Returns the base 10 logarithm of x.
19 * The argument is separated into its exponent and fractional
20 * parts. If the exponent is between -1 and +1, the logarithm
21 * of the fraction is approximated by
23 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
25 * Otherwise, setting z = 2(x-1)/x+1),
27 * log(x) = z + z^3 P(z)/Q(z).
31 * ACCURACY:
33 * Relative error:
34 * arithmetic domain # trials peak rms
35 * IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35
36 * IEEE exp(+-10000) 30000 1.0e-34 4.1e-35
38 * In the tests over the interval exp(+-10000), the logarithms
39 * of the random arguments were uniformly distributed over
40 * [-10000, +10000].
45 Cephes Math Library Release 2.2: January, 1991
46 Copyright 1984, 1991 by Stephen L. Moshier
47 Adapted for glibc November, 2001
49 This library is free software; you can redistribute it and/or
50 modify it under the terms of the GNU Lesser General Public
51 License as published by the Free Software Foundation; either
52 version 2.1 of the License, or (at your option) any later version.
54 This library is distributed in the hope that it will be useful,
55 but WITHOUT ANY WARRANTY; without even the implied warranty of
56 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
57 Lesser General Public License for more details.
59 You should have received a copy of the GNU Lesser General Public
60 License along with this library; if not, write to the Free Software
61 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
65 #include "quadmath-imp.h"
67 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
68 * 1/sqrt(2) <= x < sqrt(2)
69 * Theoretical peak relative error = 5.3e-37,
70 * relative peak error spread = 2.3e-14
72 static const __float128 P[13] =
74 1.313572404063446165910279910527789794488E4Q,
75 7.771154681358524243729929227226708890930E4Q,
76 2.014652742082537582487669938141683759923E5Q,
77 3.007007295140399532324943111654767187848E5Q,
78 2.854829159639697837788887080758954924001E5Q,
79 1.797628303815655343403735250238293741397E5Q,
80 7.594356839258970405033155585486712125861E4Q,
81 2.128857716871515081352991964243375186031E4Q,
82 3.824952356185897735160588078446136783779E3Q,
83 4.114517881637811823002128927449878962058E2Q,
84 2.321125933898420063925789532045674660756E1Q,
85 4.998469661968096229986658302195402690910E-1Q,
86 1.538612243596254322971797716843006400388E-6Q
88 static const __float128 Q[12] =
90 3.940717212190338497730839731583397586124E4Q,
91 2.626900195321832660448791748036714883242E5Q,
92 7.777690340007566932935753241556479363645E5Q,
93 1.347518538384329112529391120390701166528E6Q,
94 1.514882452993549494932585972882995548426E6Q,
95 1.158019977462989115839826904108208787040E6Q,
96 6.132189329546557743179177159925690841200E5Q,
97 2.248234257620569139969141618556349415120E5Q,
98 5.605842085972455027590989944010492125825E4Q,
99 9.147150349299596453976674231612674085381E3Q,
100 9.104928120962988414618126155557301584078E2Q,
101 4.839208193348159620282142911143429644326E1Q
102 /* 1.000000000000000000000000000000000000000E0Q, */
105 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
106 * where z = 2(x-1)/(x+1)
107 * 1/sqrt(2) <= x < sqrt(2)
108 * Theoretical peak relative error = 1.1e-35,
109 * relative peak error spread 1.1e-9
111 static const __float128 R[6] =
113 1.418134209872192732479751274970992665513E5Q,
114 -8.977257995689735303686582344659576526998E4Q,
115 2.048819892795278657810231591630928516206E4Q,
116 -2.024301798136027039250415126250455056397E3Q,
117 8.057002716646055371965756206836056074715E1Q,
118 -8.828896441624934385266096344596648080902E-1Q
120 static const __float128 S[6] =
122 1.701761051846631278975701529965589676574E6Q,
123 -1.332535117259762928288745111081235577029E6Q,
124 4.001557694070773974936904547424676279307E5Q,
125 -5.748542087379434595104154610899551484314E4Q,
126 3.998526750980007367835804959888064681098E3Q,
127 -1.186359407982897997337150403816839480438E2Q
128 /* 1.000000000000000000000000000000000000000E0Q, */
131 static const __float128
132 /* log10(2) */
133 L102A = 0.3125Q,
134 L102B = -1.14700043360188047862611052755069732318101185E-2Q,
135 /* log10(e) */
136 L10EA = 0.5Q,
137 L10EB = -6.570551809674817234887108108339491770560299E-2Q,
138 /* sqrt(2)/2 */
139 SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
143 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
145 static __float128
146 neval (__float128 x, const __float128 *p, int n)
148 __float128 y;
150 p += n;
151 y = *p--;
154 y = y * x + *p--;
156 while (--n > 0);
157 return y;
161 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
163 static __float128
164 deval (__float128 x, const __float128 *p, int n)
166 __float128 y;
168 p += n;
169 y = x + *p--;
172 y = y * x + *p--;
174 while (--n > 0);
175 return y;
180 __float128
181 log10q (__float128 x)
183 __float128 z;
184 __float128 y;
185 int e;
186 int64_t hx, lx;
188 /* Test for domain */
189 GET_FLT128_WORDS64 (hx, lx, x);
190 if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
191 return (-1.0Q / fabsq (x)); /* log10l(+-0)=-inf */
192 if (hx < 0)
193 return (x - x) / (x - x);
194 if (hx >= 0x7fff000000000000LL)
195 return (x + x);
197 if (x == 1.0Q)
198 return 0.0Q;
200 /* separate mantissa from exponent */
202 /* Note, frexp is used so that denormal numbers
203 * will be handled properly.
205 x = frexpq (x, &e);
208 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
209 * where z = 2(x-1)/x+1)
211 if ((e > 2) || (e < -2))
213 if (x < SQRTH)
214 { /* 2( 2x-1 )/( 2x+1 ) */
215 e -= 1;
216 z = x - 0.5Q;
217 y = 0.5Q * z + 0.5Q;
219 else
220 { /* 2 (x-1)/(x+1) */
221 z = x - 0.5Q;
222 z -= 0.5Q;
223 y = 0.5Q * x + 0.5Q;
225 x = z / y;
226 z = x * x;
227 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
228 goto done;
232 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
234 if (x < SQRTH)
236 e -= 1;
237 x = 2.0 * x - 1.0Q; /* 2x - 1 */
239 else
241 x = x - 1.0Q;
243 z = x * x;
244 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
245 y = y - 0.5 * z;
247 done:
249 /* Multiply log of fraction by log10(e)
250 * and base 2 exponent by log10(2).
252 z = y * L10EB;
253 z += x * L10EB;
254 z += e * L102B;
255 z += y * L10EA;
256 z += x * L10EA;
257 z += e * L102A;
258 return (z);