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 <https://www.gnu.org/licenses/>.
60 */
61
62 #include <math.h>
63 #include <math_private.h>
64 #include <libm-alias-finite.h>
65
66 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
67 * 1/sqrt(2) <= x < sqrt(2)
68 * Theoretical peak relative error = 5.3e-37,
69 * relative peak error spread = 2.3e-14
70 */
71 static const _Float128 P[13] =
72 {
73 L(1.313572404063446165910279910527789794488E4),
74 L(7.771154681358524243729929227226708890930E4),
75 L(2.014652742082537582487669938141683759923E5),
76 L(3.007007295140399532324943111654767187848E5),
77 L(2.854829159639697837788887080758954924001E5),
78 L(1.797628303815655343403735250238293741397E5),
79 L(7.594356839258970405033155585486712125861E4),
80 L(2.128857716871515081352991964243375186031E4),
81 L(3.824952356185897735160588078446136783779E3),
82 L(4.114517881637811823002128927449878962058E2),
83 L(2.321125933898420063925789532045674660756E1),
84 L(4.998469661968096229986658302195402690910E-1),
85 L(1.538612243596254322971797716843006400388E-6)
86 };
87 static const _Float128 Q[12] =
88 {
89 L(3.940717212190338497730839731583397586124E4),
90 L(2.626900195321832660448791748036714883242E5),
91 L(7.777690340007566932935753241556479363645E5),
92 L(1.347518538384329112529391120390701166528E6),
93 L(1.514882452993549494932585972882995548426E6),
94 L(1.158019977462989115839826904108208787040E6),
95 L(6.132189329546557743179177159925690841200E5),
96 L(2.248234257620569139969141618556349415120E5),
97 L(5.605842085972455027590989944010492125825E4),
98 L(9.147150349299596453976674231612674085381E3),
99 L(9.104928120962988414618126155557301584078E2),
100 L(4.839208193348159620282142911143429644326E1)
101 /* 1.000000000000000000000000000000000000000E0L, */
102 };
103
104 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
105 * where z = 2(x-1)/(x+1)
106 * 1/sqrt(2) <= x < sqrt(2)
107 * Theoretical peak relative error = 1.1e-35,
108 * relative peak error spread 1.1e-9
109 */
110 static const _Float128 R[6] =
111 {
112 L(1.418134209872192732479751274970992665513E5),
113 L(-8.977257995689735303686582344659576526998E4),
114 L(2.048819892795278657810231591630928516206E4),
115 L(-2.024301798136027039250415126250455056397E3),
116 L(8.057002716646055371965756206836056074715E1),
117 L(-8.828896441624934385266096344596648080902E-1)
118 };
119 static const _Float128 S[6] =
120 {
121 L(1.701761051846631278975701529965589676574E6),
122 L(-1.332535117259762928288745111081235577029E6),
123 L(4.001557694070773974936904547424676279307E5),
124 L(-5.748542087379434595104154610899551484314E4),
125 L(3.998526750980007367835804959888064681098E3),
126 L(-1.186359407982897997337150403816839480438E2)
127 /* 1.000000000000000000000000000000000000000E0L, */
128 };
129
130 static const _Float128
131 /* log2(e) - 1 */
132 LOG2EA = L(4.4269504088896340735992468100189213742664595E-1),
133 /* sqrt(2)/2 */
134 SQRTH = L(7.071067811865475244008443621048490392848359E-1);
135
136
137 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
138
139 static _Float128
140 neval (_Float128 x, const _Float128 *p, int n)
141 {
142 _Float128 y;
143
144 p += n;
145 y = *p--;
146 do
147 {
148 y = y * x + *p--;
149 }
150 while (--n > 0);
151 return y;
152 }
153
154
155 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
156
157 static _Float128
158 deval (_Float128 x, const _Float128 *p, int n)
159 {
160 _Float128 y;
161
162 p += n;
163 y = x + *p--;
164 do
165 {
166 y = y * x + *p--;
167 }
168 while (--n > 0);
169 return y;
170 }
171
172
173
174 _Float128
175 __ieee754_log2l (_Float128 x)
176 {
177 _Float128 z;
178 _Float128 y;
179 int e;
180 int64_t hx, lx;
181
182 /* Test for domain */
183 GET_LDOUBLE_WORDS64 (hx, lx, x);
184 if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
185 return (-1 / fabsl (x)); /* log2l(+-0)=-inf */
186 if (hx < 0)
187 return (x - x) / (x - x);
188 if (hx >= 0x7fff000000000000LL)
189 return (x + x);
190
191 if (x == 1)
192 return 0;
193
194 /* separate mantissa from exponent */
195
196 /* Note, frexp is used so that denormal numbers
197 * will be handled properly.
198 */
199 x = __frexpl (x, &e);
200
201
202 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
203 * where z = 2(x-1)/x+1)
204 */
205 if ((e > 2) || (e < -2))
206 {
207 if (x < SQRTH)
208 { /* 2( 2x-1 )/( 2x+1 ) */
209 e -= 1;
210 z = x - L(0.5);
211 y = L(0.5) * z + L(0.5);
212 }
213 else
214 { /* 2 (x-1)/(x+1) */
215 z = x - L(0.5);
216 z -= L(0.5);
217 y = L(0.5) * x + L(0.5);
218 }
219 x = z / y;
220 z = x * x;
221 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
222 goto done;
223 }
224
225
226 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
227
228 if (x < SQRTH)
229 {
230 e -= 1;
231 x = 2.0 * x - 1; /* 2x - 1 */
232 }
233 else
234 {
235 x = x - 1;
236 }
237 z = x * x;
238 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
239 y = y - 0.5 * z;
240
241 done:
242
243 /* Multiply log of fraction by log2(e)
244 * and base 2 exponent by 1
245 */
246 z = y * LOG2EA;
247 z += x * LOG2EA;
248 z += y;
249 z += x;
250 z += e;
251 return (z);
252 }
253 libm_alias_finite (__ieee754_log2l, __log2l)