1 /* log2l.c
2 * Base 2 logarithm, 128-bit long double precision
3 *
4 *
5 *
6 * SYNOPSIS:
7 *
8 * long double x, y, log2l();
9 *
10 * y = log2l( x );
11 *
12 *
13 *
14 * DESCRIPTION:
15 *
16 * Returns the base 2 logarithm of x.
17 *
18 * The argument is separated into its exponent and fractional
19 * parts. If the exponent is between -1 and +1, the (natural)
20 * logarithm of the fraction is approximated by
21 *
22 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
23 *
24 * Otherwise, setting z = 2(x-1)/x+1),
25 *
26 * log(x) = z + z^3 P(z)/Q(z).
27 *
28 *
29 *
30 * ACCURACY:
31 *
32 * Relative error:
33 * arithmetic domain # trials peak rms
34 * IEEE 0.5, 2.0 100,000 2.6e-34 4.9e-35
35 * IEEE exp(+-10000) 100,000 9.6e-35 4.0e-35
36 *
37 * In the tests over the interval exp(+-10000), the logarithms
38 * of the random arguments were uniformly distributed over
39 * [-10000, +10000].
40 *
41 */
42
43 /*
44 Cephes Math Library Release 2.2: January, 1991
45 Copyright 1984, 1991 by Stephen L. Moshier
46 Adapted for glibc November, 2001
47
48 This library is free software; you can redistribute it and/or
49 modify it under the terms of the GNU Lesser General Public
50 License as published by the Free Software Foundation; either
51 version 2.1 of the License, or (at your option) any later version.
52
53 This library is distributed in the hope that it will be useful,
54 but WITHOUT ANY WARRANTY; without even the implied warranty of
55 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
56 Lesser General Public License for more details.
57
58 You should have received a copy of the GNU Lesser General Public
59 License along with this library; if not, see <http://www.gnu.org/licenses/>.
60 */
61
62 #include "quadmath-imp.h"
63
64 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
65 * 1/sqrt(2) <= x < sqrt(2)
66 * Theoretical peak relative error = 5.3e-37,
67 * relative peak error spread = 2.3e-14
68 */
69 static const __float128 P[13] =
70 {
71 1.313572404063446165910279910527789794488E4Q,
72 7.771154681358524243729929227226708890930E4Q,
73 2.014652742082537582487669938141683759923E5Q,
74 3.007007295140399532324943111654767187848E5Q,
75 2.854829159639697837788887080758954924001E5Q,
76 1.797628303815655343403735250238293741397E5Q,
77 7.594356839258970405033155585486712125861E4Q,
78 2.128857716871515081352991964243375186031E4Q,
79 3.824952356185897735160588078446136783779E3Q,
80 4.114517881637811823002128927449878962058E2Q,
81 2.321125933898420063925789532045674660756E1Q,
82 4.998469661968096229986658302195402690910E-1Q,
83 1.538612243596254322971797716843006400388E-6Q
84 };
85 static const __float128 Q[12] =
86 {
87 3.940717212190338497730839731583397586124E4Q,
88 2.626900195321832660448791748036714883242E5Q,
89 7.777690340007566932935753241556479363645E5Q,
90 1.347518538384329112529391120390701166528E6Q,
91 1.514882452993549494932585972882995548426E6Q,
92 1.158019977462989115839826904108208787040E6Q,
93 6.132189329546557743179177159925690841200E5Q,
94 2.248234257620569139969141618556349415120E5Q,
95 5.605842085972455027590989944010492125825E4Q,
96 9.147150349299596453976674231612674085381E3Q,
97 9.104928120962988414618126155557301584078E2Q,
98 4.839208193348159620282142911143429644326E1Q
99 /* 1.000000000000000000000000000000000000000E0L, */
100 };
101
102 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
103 * where z = 2(x-1)/(x+1)
104 * 1/sqrt(2) <= x < sqrt(2)
105 * Theoretical peak relative error = 1.1e-35,
106 * relative peak error spread 1.1e-9
107 */
108 static const __float128 R[6] =
109 {
110 1.418134209872192732479751274970992665513E5Q,
111 -8.977257995689735303686582344659576526998E4Q,
112 2.048819892795278657810231591630928516206E4Q,
113 -2.024301798136027039250415126250455056397E3Q,
114 8.057002716646055371965756206836056074715E1Q,
115 -8.828896441624934385266096344596648080902E-1Q
116 };
117 static const __float128 S[6] =
118 {
119 1.701761051846631278975701529965589676574E6Q,
120 -1.332535117259762928288745111081235577029E6Q,
121 4.001557694070773974936904547424676279307E5Q,
122 -5.748542087379434595104154610899551484314E4Q,
123 3.998526750980007367835804959888064681098E3Q,
124 -1.186359407982897997337150403816839480438E2Q
125 /* 1.000000000000000000000000000000000000000E0L, */
126 };
127
128 static const __float128
129 /* log2(e) - 1 */
130 LOG2EA = 4.4269504088896340735992468100189213742664595E-1Q,
131 /* sqrt(2)/2 */
132 SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
133
134
135 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
136
137 static __float128
138 neval (__float128 x, const __float128 *p, int n)
139 {
140 __float128 y;
141
142 p += n;
143 y = *p--;
144 do
145 {
146 y = y * x + *p--;
147 }
148 while (--n > 0);
149 return y;
150 }
151
152
153 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
154
155 static __float128
156 deval (__float128 x, const __float128 *p, int n)
157 {
158 __float128 y;
159
160 p += n;
161 y = x + *p--;
162 do
163 {
164 y = y * x + *p--;
165 }
166 while (--n > 0);
167 return y;
168 }
169
170
171
172 __float128
173 log2q (__float128 x)
174 {
175 __float128 z;
176 __float128 y;
177 int e;
178 int64_t hx, lx;
179
180 /* Test for domain */
181 GET_FLT128_WORDS64 (hx, lx, x);
182 if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
183 return (-1 / fabsq (x)); /* log2l(+-0)=-inf */
184 if (hx < 0)
185 return (x - x) / (x - x);
186 if (hx >= 0x7fff000000000000LL)
187 return (x + x);
188
189 if (x == 1)
190 return 0;
191
192 /* separate mantissa from exponent */
193
194 /* Note, frexp is used so that denormal numbers
195 * will be handled properly.
196 */
197 x = frexpq (x, &e);
198
199
200 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
201 * where z = 2(x-1)/x+1)
202 */
203 if ((e > 2) || (e < -2))
204 {
205 if (x < SQRTH)
206 { /* 2( 2x-1 )/( 2x+1 ) */
207 e -= 1;
208 z = x - 0.5Q;
209 y = 0.5Q * z + 0.5Q;
210 }
211 else
212 { /* 2 (x-1)/(x+1) */
213 z = x - 0.5Q;
214 z -= 0.5Q;
215 y = 0.5Q * x + 0.5Q;
216 }
217 x = z / y;
218 z = x * x;
219 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
220 goto done;
221 }
222
223
224 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
225
226 if (x < SQRTH)
227 {
228 e -= 1;
229 x = 2.0 * x - 1; /* 2x - 1 */
230 }
231 else
232 {
233 x = x - 1;
234 }
235 z = x * x;
236 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
237 y = y - 0.5 * z;
238
239 done:
240
241 /* Multiply log of fraction by log2(e)
242 * and base 2 exponent by 1
243 */
244 z = y * LOG2EA;
245 z += x * LOG2EA;
246 z += y;
247 z += x;
248 z += e;
249 return (z);
250 }