1/* log10q.c
2 *
3 * Common logarithm, 128-bit __float128 precision
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * __float128 x, y, log10l();
10 *
11 * y = log10q( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns the base 10 logarithm of x.
18 *
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
22 *
23 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
24 *
25 * Otherwise, setting z = 2(x-1)/x+1),
26 *
27 * log(x) = z + z^3 P(z)/Q(z).
28 *
29 *
30 *
31 * ACCURACY:
32 *
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
37 *
38 * In the tests over the interval exp(+-10000), the logarithms
39 * of the random arguments were uniformly distributed over
40 * [-10000, +10000].
41 *
42 */
43
44/*
45 Cephes Math Library Release 2.2: January, 1991
46 Copyright 1984, 1991 by Stephen L. Moshier
47 Adapted for glibc November, 2001
48
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.
53
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.
58
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
62
63 */
64
65#include "quadmath-imp.h"
66
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
71 */
72static const __float128 P[13] =
73{
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
87};
88static const __float128 Q[12] =
89{
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, */
103};
104
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
110 */
111static const __float128 R[6] =
112{
113 1.418134209872192732479751274970992665513E5Q,
114 -8.977257995689735303686582344659576526998E4Q,
115 2.048819892795278657810231591630928516206E4Q,
116 -2.024301798136027039250415126250455056397E3Q,
117 8.057002716646055371965756206836056074715E1Q,
118 -8.828896441624934385266096344596648080902E-1Q
119};
120static const __float128 S[6] =
121{
122 1.701761051846631278975701529965589676574E6Q,
123 -1.332535117259762928288745111081235577029E6Q,
124 4.001557694070773974936904547424676279307E5Q,
125 -5.748542087379434595104154610899551484314E4Q,
126 3.998526750980007367835804959888064681098E3Q,
127 -1.186359407982897997337150403816839480438E2Q
128/* 1.000000000000000000000000000000000000000E0Q, */
129};
130
131static const __float128
132/* log10(2) */
133L102A = 0.3125Q,
134L102B = -1.14700043360188047862611052755069732318101185E-2Q,
135/* log10(e) */
136L10EA = 0.5Q,
137L10EB = -6.570551809674817234887108108339491770560299E-2Q,
138/* sqrt(2)/2 */
139SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
140
141
142
143/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
144
145static __float128
146neval (__float128 x, const __float128 *p, int n)
147{
148 __float128 y;
149
150 p += n;
151 y = *p--;
152 do
153 {
154 y = y * x + *p--;
155 }
156 while (--n > 0);
157 return y;
158}
159
160
161/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
162
163static __float128
164deval (__float128 x, const __float128 *p, int n)
165{
166 __float128 y;
167
168 p += n;
169 y = x + *p--;
170 do
171 {
172 y = y * x + *p--;
173 }
174 while (--n > 0);
175 return y;
176}
177
178
179
180__float128
181log10q (__float128 x)
182{
183 __float128 z;
184 __float128 y;
185 int e;
186 int64_t hx, lx;
187
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);
196
197 if (x == 1.0Q)
198 return 0.0Q;
199
200/* separate mantissa from exponent */
201
202/* Note, frexp is used so that denormal numbers
203 * will be handled properly.
204 */
205 x = frexpq (x, &e);
206
207
208/* logarithm using log(x) = z + z**3 P(z)/Q(z),
209 * where z = 2(x-1)/x+1)
210 */
211 if ((e > 2) || (e < -2))
212 {
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;
218 }
219 else
220 { /* 2 (x-1)/(x+1) */
221 z = x - 0.5Q;
222 z -= 0.5Q;
223 y = 0.5Q * x + 0.5Q;
224 }
225 x = z / y;
226 z = x * x;
227 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
228 goto done;
229 }
230
231
232/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
233
234 if (x < SQRTH)
235 {
236 e -= 1;
237 x = 2.0 * x - 1.0Q; /* 2x - 1 */
238 }
239 else
240 {
241 x = x - 1.0Q;
242 }
243 z = x * x;
244 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
245 y = y - 0.5 * z;
246
247done:
248
249 /* Multiply log of fraction by log10(e)
250 * and base 2 exponent by log10(2).
251 */
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);
259}
260