1 /* logll.c
2 *
3 * Natural logarithm for 128-bit long double precision.
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * long double x, y, logl();
10 *
11 * y = logl( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns the base e (2.718...) logarithm of x.
18 *
19 * The argument is separated into its exponent and fractional
20 * parts. Use of a lookup table increases the speed of the routine.
21 * The program uses logarithms tabulated at intervals of 1/128 to
22 * cover the domain from approximately 0.7 to 1.4.
23 *
24 * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
25 * log(1+x) = x - 0.5 x^2 + x^3 P(x) .
26 *
27 *
28 *
29 * ACCURACY:
30 *
31 * Relative error:
32 * arithmetic domain # trials peak rms
33 * IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35
34 * IEEE 0.125, 8 100000 1.2e-34 4.1e-35
35 *
36 *
37 * WARNING:
38 *
39 * This program uses integer operations on bit fields of floating-point
40 * numbers. It does not work with data structures other than the
41 * structure assumed.
42 *
43 */
44
45 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
46
47 This library is free software; you can redistribute it and/or
48 modify it under the terms of the GNU Lesser General Public
49 License as published by the Free Software Foundation; either
50 version 2.1 of the License, or (at your option) any later version.
51
52 This library is distributed in the hope that it will be useful,
53 but WITHOUT ANY WARRANTY; without even the implied warranty of
54 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
55 Lesser General Public License for more details.
56
57 You should have received a copy of the GNU Lesser General Public
58 License along with this library; if not, see
59 <https://www.gnu.org/licenses/>. */
60
61 #include <math.h>
62 #include <math_private.h>
63 #include <libm-alias-finite.h>
64
65 /* log(1+x) = x - .5 x^2 + x^3 l(x)
66 -.0078125 <= x <= +.0078125
67 peak relative error 1.2e-37 */
68 static const long double
69 l3 = 3.333333333333333333333333333333336096926E-1L,
70 l4 = -2.499999999999999999999999999486853077002E-1L,
71 l5 = 1.999999999999999999999999998515277861905E-1L,
72 l6 = -1.666666666666666666666798448356171665678E-1L,
73 l7 = 1.428571428571428571428808945895490721564E-1L,
74 l8 = -1.249999999999999987884655626377588149000E-1L,
75 l9 = 1.111111111111111093947834982832456459186E-1L,
76 l10 = -1.000000000000532974938900317952530453248E-1L,
77 l11 = 9.090909090915566247008015301349979892689E-2L,
78 l12 = -8.333333211818065121250921925397567745734E-2L,
79 l13 = 7.692307559897661630807048686258659316091E-2L,
80 l14 = -7.144242754190814657241902218399056829264E-2L,
81 l15 = 6.668057591071739754844678883223432347481E-2L;
82
83 /* Lookup table of ln(t) - (t-1)
84 t = 0.5 + (k+26)/128)
85 k = 0, ..., 91 */
86 static const long double logtbl[92] = {
87 -5.5345593589352099112142921677820359632418E-2L,
88 -5.2108257402767124761784665198737642086148E-2L,
89 -4.8991686870576856279407775480686721935120E-2L,
90 -4.5993270766361228596215288742353061431071E-2L,
91 -4.3110481649613269682442058976885699556950E-2L,
92 -4.0340872319076331310838085093194799765520E-2L,
93 -3.7682072451780927439219005993827431503510E-2L,
94 -3.5131785416234343803903228503274262719586E-2L,
95 -3.2687785249045246292687241862699949178831E-2L,
96 -3.0347913785027239068190798397055267411813E-2L,
97 -2.8110077931525797884641940838507561326298E-2L,
98 -2.5972247078357715036426583294246819637618E-2L,
99 -2.3932450635346084858612873953407168217307E-2L,
100 -2.1988775689981395152022535153795155900240E-2L,
101 -2.0139364778244501615441044267387667496733E-2L,
102 -1.8382413762093794819267536615342902718324E-2L,
103 -1.6716169807550022358923589720001638093023E-2L,
104 -1.5138929457710992616226033183958974965355E-2L,
105 -1.3649036795397472900424896523305726435029E-2L,
106 -1.2244881690473465543308397998034325468152E-2L,
107 -1.0924898127200937840689817557742469105693E-2L,
108 -9.6875626072830301572839422532631079809328E-3L,
109 -8.5313926245226231463436209313499745894157E-3L,
110 -7.4549452072765973384933565912143044991706E-3L,
111 -6.4568155251217050991200599386801665681310E-3L,
112 -5.5356355563671005131126851708522185605193E-3L,
113 -4.6900728132525199028885749289712348829878E-3L,
114 -3.9188291218610470766469347968659624282519E-3L,
115 -3.2206394539524058873423550293617843896540E-3L,
116 -2.5942708080877805657374888909297113032132E-3L,
117 -2.0385211375711716729239156839929281289086E-3L,
118 -1.5522183228760777967376942769773768850872E-3L,
119 -1.1342191863606077520036253234446621373191E-3L,
120 -7.8340854719967065861624024730268350459991E-4L,
121 -4.9869831458030115699628274852562992756174E-4L,
122 -2.7902661731604211834685052867305795169688E-4L,
123 -1.2335696813916860754951146082826952093496E-4L,
124 -3.0677461025892873184042490943581654591817E-5L,
125 #define ZERO logtbl[38]
126 0.0000000000000000000000000000000000000000E0L,
127 -3.0359557945051052537099938863236321874198E-5L,
128 -1.2081346403474584914595395755316412213151E-4L,
129 -2.7044071846562177120083903771008342059094E-4L,
130 -4.7834133324631162897179240322783590830326E-4L,
131 -7.4363569786340080624467487620270965403695E-4L,
132 -1.0654639687057968333207323853366578860679E-3L,
133 -1.4429854811877171341298062134712230604279E-3L,
134 -1.8753781835651574193938679595797367137975E-3L,
135 -2.3618380914922506054347222273705859653658E-3L,
136 -2.9015787624124743013946600163375853631299E-3L,
137 -3.4938307889254087318399313316921940859043E-3L,
138 -4.1378413103128673800485306215154712148146E-3L,
139 -4.8328735414488877044289435125365629849599E-3L,
140 -5.5782063183564351739381962360253116934243E-3L,
141 -6.3731336597098858051938306767880719015261E-3L,
142 -7.2169643436165454612058905294782949315193E-3L,
143 -8.1090214990427641365934846191367315083867E-3L,
144 -9.0486422112807274112838713105168375482480E-3L,
145 -1.0035177140880864314674126398350812606841E-2L,
146 -1.1067990155502102718064936259435676477423E-2L,
147 -1.2146457974158024928196575103115488672416E-2L,
148 -1.3269969823361415906628825374158424754308E-2L,
149 -1.4437927104692837124388550722759686270765E-2L,
150 -1.5649743073340777659901053944852735064621E-2L,
151 -1.6904842527181702880599758489058031645317E-2L,
152 -1.8202661505988007336096407340750378994209E-2L,
153 -1.9542647000370545390701192438691126552961E-2L,
154 -2.0924256670080119637427928803038530924742E-2L,
155 -2.2346958571309108496179613803760727786257E-2L,
156 -2.3810230892650362330447187267648486279460E-2L,
157 -2.5313561699385640380910474255652501521033E-2L,
158 -2.6856448685790244233704909690165496625399E-2L,
159 -2.8438398935154170008519274953860128449036E-2L,
160 -3.0058928687233090922411781058956589863039E-2L,
161 -3.1717563112854831855692484086486099896614E-2L,
162 -3.3413836095418743219397234253475252001090E-2L,
163 -3.5147290019036555862676702093393332533702E-2L,
164 -3.6917475563073933027920505457688955423688E-2L,
165 -3.8723951502862058660874073462456610731178E-2L,
166 -4.0566284516358241168330505467000838017425E-2L,
167 -4.2444048996543693813649967076598766917965E-2L,
168 -4.4356826869355401653098777649745233339196E-2L,
169 -4.6304207416957323121106944474331029996141E-2L,
170 -4.8285787106164123613318093945035804818364E-2L,
171 -5.0301169421838218987124461766244507342648E-2L,
172 -5.2349964705088137924875459464622098310997E-2L,
173 -5.4431789996103111613753440311680967840214E-2L,
174 -5.6546268881465384189752786409400404404794E-2L,
175 -5.8693031345788023909329239565012647817664E-2L,
176 -6.0871713627532018185577188079210189048340E-2L,
177 -6.3081958078862169742820420185833800925568E-2L,
178 -6.5323413029406789694910800219643791556918E-2L,
179 -6.7595732653791419081537811574227049288168E-2L
180 };
181
182 /* ln(2) = ln2a + ln2b with extended precision. */
183 static const long double
184 ln2a = 6.93145751953125e-1L,
185 ln2b = 1.4286068203094172321214581765680755001344E-6L;
186
187 static const long double
188 ldbl_epsilon = 0x1p-106L;
189
190 long double
191 __ieee754_logl(long double x)
192 {
193 long double z, y, w, t;
194 unsigned int m;
195 int k, e;
196 double xhi;
197 uint32_t hx, lx;
198
199 xhi = ldbl_high (x);
200 EXTRACT_WORDS (hx, lx, xhi);
201 m = hx;
202
203 /* Check for IEEE special cases. */
204 k = m & 0x7fffffff;
205 /* log(0) = -infinity. */
206 if ((k | lx) == 0)
207 {
208 return -0.5L / ZERO;
209 }
210 /* log ( x < 0 ) = NaN */
211 if (m & 0x80000000)
212 {
213 return (x - x) / ZERO;
214 }
215 /* log (infinity or NaN) */
216 if (k >= 0x7ff00000)
217 {
218 return x + x;
219 }
220
221 /* On this interval the table is not used due to cancellation error. */
222 if ((x <= 1.0078125L) && (x >= 0.9921875L))
223 {
224 if (x == 1.0L)
225 return 0.0L;
226 z = x - 1.0L;
227 k = 64;
228 t = 1.0L;
229 e = 0;
230 }
231 else
232 {
233 /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */
234 unsigned int w0;
235 e = (int) (m >> 20) - (int) 0x3fe;
236 if (e == -1022)
237 {
238 x *= 0x1p106L;
239 xhi = ldbl_high (x);
240 EXTRACT_WORDS (hx, lx, xhi);
241 m = hx;
242 e = (int) (m >> 20) - (int) 0x3fe - 106;
243 }
244 m &= 0xfffff;
245 w0 = m | 0x3fe00000;
246 m |= 0x100000;
247 /* Find lookup table index k from high order bits of the significand. */
248 if (m < 0x168000)
249 {
250 k = (m - 0xff000) >> 13;
251 /* t is the argument 0.5 + (k+26)/128
252 of the nearest item to u in the lookup table. */
253 INSERT_WORDS (xhi, 0x3ff00000 + (k << 13), 0);
254 t = xhi;
255 w0 += 0x100000;
256 e -= 1;
257 k += 64;
258 }
259 else
260 {
261 k = (m - 0xfe000) >> 14;
262 INSERT_WORDS (xhi, 0x3fe00000 + (k << 14), 0);
263 t = xhi;
264 }
265 x = __scalbnl (x, ((int) ((w0 - hx) * 2)) >> 21);
266 /* log(u) = log( t u/t ) = log(t) + log(u/t)
267 log(t) is tabulated in the lookup table.
268 Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t.
269 cf. Cody & Waite. */
270 z = (x - t) / t;
271 }
272 /* Series expansion of log(1+z). */
273 w = z * z;
274 /* Avoid spurious underflows. */
275 if (__glibc_unlikely (fabsl (z) <= ldbl_epsilon))
276 y = 0.0L;
277 else
278 {
279 y = ((((((((((((l15 * z
280 + l14) * z
281 + l13) * z
282 + l12) * z
283 + l11) * z
284 + l10) * z
285 + l9) * z
286 + l8) * z
287 + l7) * z
288 + l6) * z
289 + l5) * z
290 + l4) * z
291 + l3) * z * w;
292 y -= 0.5 * w;
293 }
294 y += e * ln2b; /* Base 2 exponent offset times ln(2). */
295 y += z;
296 y += logtbl[k-26]; /* log(t) - (t-1) */
297 y += (t - 1.0L);
298 y += e * ln2a;
299 return y;
300 }
301 libm_alias_finite (__ieee754_logl, __logl)