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 _Float128 |
69 | l3 = L(3.333333333333333333333333333333336096926E-1), |
70 | l4 = L(-2.499999999999999999999999999486853077002E-1), |
71 | l5 = L(1.999999999999999999999999998515277861905E-1), |
72 | l6 = L(-1.666666666666666666666798448356171665678E-1), |
73 | l7 = L(1.428571428571428571428808945895490721564E-1), |
74 | l8 = L(-1.249999999999999987884655626377588149000E-1), |
75 | l9 = L(1.111111111111111093947834982832456459186E-1), |
76 | l10 = L(-1.000000000000532974938900317952530453248E-1), |
77 | l11 = L(9.090909090915566247008015301349979892689E-2), |
78 | l12 = L(-8.333333211818065121250921925397567745734E-2), |
79 | l13 = L(7.692307559897661630807048686258659316091E-2), |
80 | l14 = L(-7.144242754190814657241902218399056829264E-2), |
81 | l15 = L(6.668057591071739754844678883223432347481E-2); |
82 | |
83 | /* Lookup table of ln(t) - (t-1) |
84 | t = 0.5 + (k+26)/128) |
85 | k = 0, ..., 91 */ |
86 | static const _Float128 logtbl[92] = { |
87 | L(-5.5345593589352099112142921677820359632418E-2), |
88 | L(-5.2108257402767124761784665198737642086148E-2), |
89 | L(-4.8991686870576856279407775480686721935120E-2), |
90 | L(-4.5993270766361228596215288742353061431071E-2), |
91 | L(-4.3110481649613269682442058976885699556950E-2), |
92 | L(-4.0340872319076331310838085093194799765520E-2), |
93 | L(-3.7682072451780927439219005993827431503510E-2), |
94 | L(-3.5131785416234343803903228503274262719586E-2), |
95 | L(-3.2687785249045246292687241862699949178831E-2), |
96 | L(-3.0347913785027239068190798397055267411813E-2), |
97 | L(-2.8110077931525797884641940838507561326298E-2), |
98 | L(-2.5972247078357715036426583294246819637618E-2), |
99 | L(-2.3932450635346084858612873953407168217307E-2), |
100 | L(-2.1988775689981395152022535153795155900240E-2), |
101 | L(-2.0139364778244501615441044267387667496733E-2), |
102 | L(-1.8382413762093794819267536615342902718324E-2), |
103 | L(-1.6716169807550022358923589720001638093023E-2), |
104 | L(-1.5138929457710992616226033183958974965355E-2), |
105 | L(-1.3649036795397472900424896523305726435029E-2), |
106 | L(-1.2244881690473465543308397998034325468152E-2), |
107 | L(-1.0924898127200937840689817557742469105693E-2), |
108 | L(-9.6875626072830301572839422532631079809328E-3), |
109 | L(-8.5313926245226231463436209313499745894157E-3), |
110 | L(-7.4549452072765973384933565912143044991706E-3), |
111 | L(-6.4568155251217050991200599386801665681310E-3), |
112 | L(-5.5356355563671005131126851708522185605193E-3), |
113 | L(-4.6900728132525199028885749289712348829878E-3), |
114 | L(-3.9188291218610470766469347968659624282519E-3), |
115 | L(-3.2206394539524058873423550293617843896540E-3), |
116 | L(-2.5942708080877805657374888909297113032132E-3), |
117 | L(-2.0385211375711716729239156839929281289086E-3), |
118 | L(-1.5522183228760777967376942769773768850872E-3), |
119 | L(-1.1342191863606077520036253234446621373191E-3), |
120 | L(-7.8340854719967065861624024730268350459991E-4), |
121 | L(-4.9869831458030115699628274852562992756174E-4), |
122 | L(-2.7902661731604211834685052867305795169688E-4), |
123 | L(-1.2335696813916860754951146082826952093496E-4), |
124 | L(-3.0677461025892873184042490943581654591817E-5), |
125 | #define ZERO logtbl[38] |
126 | L(0.0000000000000000000000000000000000000000E0), |
127 | L(-3.0359557945051052537099938863236321874198E-5), |
128 | L(-1.2081346403474584914595395755316412213151E-4), |
129 | L(-2.7044071846562177120083903771008342059094E-4), |
130 | L(-4.7834133324631162897179240322783590830326E-4), |
131 | L(-7.4363569786340080624467487620270965403695E-4), |
132 | L(-1.0654639687057968333207323853366578860679E-3), |
133 | L(-1.4429854811877171341298062134712230604279E-3), |
134 | L(-1.8753781835651574193938679595797367137975E-3), |
135 | L(-2.3618380914922506054347222273705859653658E-3), |
136 | L(-2.9015787624124743013946600163375853631299E-3), |
137 | L(-3.4938307889254087318399313316921940859043E-3), |
138 | L(-4.1378413103128673800485306215154712148146E-3), |
139 | L(-4.8328735414488877044289435125365629849599E-3), |
140 | L(-5.5782063183564351739381962360253116934243E-3), |
141 | L(-6.3731336597098858051938306767880719015261E-3), |
142 | L(-7.2169643436165454612058905294782949315193E-3), |
143 | L(-8.1090214990427641365934846191367315083867E-3), |
144 | L(-9.0486422112807274112838713105168375482480E-3), |
145 | L(-1.0035177140880864314674126398350812606841E-2), |
146 | L(-1.1067990155502102718064936259435676477423E-2), |
147 | L(-1.2146457974158024928196575103115488672416E-2), |
148 | L(-1.3269969823361415906628825374158424754308E-2), |
149 | L(-1.4437927104692837124388550722759686270765E-2), |
150 | L(-1.5649743073340777659901053944852735064621E-2), |
151 | L(-1.6904842527181702880599758489058031645317E-2), |
152 | L(-1.8202661505988007336096407340750378994209E-2), |
153 | L(-1.9542647000370545390701192438691126552961E-2), |
154 | L(-2.0924256670080119637427928803038530924742E-2), |
155 | L(-2.2346958571309108496179613803760727786257E-2), |
156 | L(-2.3810230892650362330447187267648486279460E-2), |
157 | L(-2.5313561699385640380910474255652501521033E-2), |
158 | L(-2.6856448685790244233704909690165496625399E-2), |
159 | L(-2.8438398935154170008519274953860128449036E-2), |
160 | L(-3.0058928687233090922411781058956589863039E-2), |
161 | L(-3.1717563112854831855692484086486099896614E-2), |
162 | L(-3.3413836095418743219397234253475252001090E-2), |
163 | L(-3.5147290019036555862676702093393332533702E-2), |
164 | L(-3.6917475563073933027920505457688955423688E-2), |
165 | L(-3.8723951502862058660874073462456610731178E-2), |
166 | L(-4.0566284516358241168330505467000838017425E-2), |
167 | L(-4.2444048996543693813649967076598766917965E-2), |
168 | L(-4.4356826869355401653098777649745233339196E-2), |
169 | L(-4.6304207416957323121106944474331029996141E-2), |
170 | L(-4.8285787106164123613318093945035804818364E-2), |
171 | L(-5.0301169421838218987124461766244507342648E-2), |
172 | L(-5.2349964705088137924875459464622098310997E-2), |
173 | L(-5.4431789996103111613753440311680967840214E-2), |
174 | L(-5.6546268881465384189752786409400404404794E-2), |
175 | L(-5.8693031345788023909329239565012647817664E-2), |
176 | L(-6.0871713627532018185577188079210189048340E-2), |
177 | L(-6.3081958078862169742820420185833800925568E-2), |
178 | L(-6.5323413029406789694910800219643791556918E-2), |
179 | L(-6.7595732653791419081537811574227049288168E-2) |
180 | }; |
181 | |
182 | /* ln(2) = ln2a + ln2b with extended precision. */ |
183 | static const _Float128 |
184 | ln2a = L(6.93145751953125e-1), |
185 | ln2b = L(1.4286068203094172321214581765680755001344E-6); |
186 | |
187 | _Float128 |
188 | __ieee754_logl(_Float128 x) |
189 | { |
190 | _Float128 z, y, w; |
191 | ieee854_long_double_shape_type u, t; |
192 | unsigned int m; |
193 | int k, e; |
194 | |
195 | u.value = x; |
196 | m = u.parts32.w0; |
197 | |
198 | /* Check for IEEE special cases. */ |
199 | k = m & 0x7fffffff; |
200 | /* log(0) = -infinity. */ |
201 | if ((k | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0) |
202 | { |
203 | return L(-0.5) / ZERO; |
204 | } |
205 | /* log ( x < 0 ) = NaN */ |
206 | if (m & 0x80000000) |
207 | { |
208 | return (x - x) / ZERO; |
209 | } |
210 | /* log (infinity or NaN) */ |
211 | if (k >= 0x7fff0000) |
212 | { |
213 | return x + x; |
214 | } |
215 | |
216 | /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ |
217 | u.value = __frexpl (x, &e); |
218 | m = u.parts32.w0 & 0xffff; |
219 | m |= 0x10000; |
220 | /* Find lookup table index k from high order bits of the significand. */ |
221 | if (m < 0x16800) |
222 | { |
223 | k = (m - 0xff00) >> 9; |
224 | /* t is the argument 0.5 + (k+26)/128 |
225 | of the nearest item to u in the lookup table. */ |
226 | t.parts32.w0 = 0x3fff0000 + (k << 9); |
227 | t.parts32.w1 = 0; |
228 | t.parts32.w2 = 0; |
229 | t.parts32.w3 = 0; |
230 | u.parts32.w0 += 0x10000; |
231 | e -= 1; |
232 | k += 64; |
233 | } |
234 | else |
235 | { |
236 | k = (m - 0xfe00) >> 10; |
237 | t.parts32.w0 = 0x3ffe0000 + (k << 10); |
238 | t.parts32.w1 = 0; |
239 | t.parts32.w2 = 0; |
240 | t.parts32.w3 = 0; |
241 | } |
242 | /* On this interval the table is not used due to cancellation error. */ |
243 | if ((x <= L(1.0078125)) && (x >= L(0.9921875))) |
244 | { |
245 | if (x == 1) |
246 | return 0; |
247 | z = x - 1; |
248 | k = 64; |
249 | t.value = 1; |
250 | e = 0; |
251 | } |
252 | else |
253 | { |
254 | /* log(u) = log( t u/t ) = log(t) + log(u/t) |
255 | log(t) is tabulated in the lookup table. |
256 | Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t. |
257 | cf. Cody & Waite. */ |
258 | z = (u.value - t.value) / t.value; |
259 | } |
260 | /* Series expansion of log(1+z). */ |
261 | w = z * z; |
262 | y = ((((((((((((l15 * z |
263 | + l14) * z |
264 | + l13) * z |
265 | + l12) * z |
266 | + l11) * z |
267 | + l10) * z |
268 | + l9) * z |
269 | + l8) * z |
270 | + l7) * z |
271 | + l6) * z |
272 | + l5) * z |
273 | + l4) * z |
274 | + l3) * z * w; |
275 | y -= 0.5 * w; |
276 | y += e * ln2b; /* Base 2 exponent offset times ln(2). */ |
277 | y += z; |
278 | y += logtbl[k-26]; /* log(t) - (t-1) */ |
279 | y += (t.value - 1); |
280 | y += e * ln2a; |
281 | return y; |
282 | } |
283 | libm_alias_finite (__ieee754_logl, __logl) |
284 | |