Branch data Line data Source code
1 : : /*
2 : : * Copyright (c) 1992, 1993
3 : : * The Regents of the University of California. All rights reserved.
4 : : *
5 : : * Redistribution and use in source and binary forms, with or without
6 : : * modification, are permitted provided that the following conditions
7 : : * are met:
8 : : * 1. Redistributions of source code must retain the above copyright
9 : : * notice, this list of conditions and the following disclaimer.
10 : : * 2. Redistributions in binary form must reproduce the above copyright
11 : : * notice, this list of conditions and the following disclaimer in the
12 : : * documentation and/or other materials provided with the distribution.
13 : : * 3. Neither the name of the University nor the names of its contributors
14 : : * may be used to endorse or promote products derived from this software
15 : : * without specific prior written permission.
16 : : *
17 : : * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
18 : : * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 : : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 : : * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
21 : : * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22 : : * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23 : : * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 : : * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25 : : * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26 : : * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27 : : * SUCH DAMAGE.
28 : : */
29 : :
30 : : /* @(#)log.c 8.2 (Berkeley) 11/30/93 */
31 : : #include "cdefs-compat.h"
32 : : //__FBSDID("$FreeBSD: src/lib/msun/bsdsrc/b_log.c,v 1.9 2008/02/22 02:26:51 das Exp $");
33 : :
34 : : #include <openlibm_math.h>
35 : :
36 : : #include "mathimpl.h"
37 : :
38 : : /* Table-driven natural logarithm.
39 : : *
40 : : * This code was derived, with minor modifications, from:
41 : : * Peter Tang, "Table-Driven Implementation of the
42 : : * Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
43 : : * Math Software, vol 16. no 4, pp 378-400, Dec 1990).
44 : : *
45 : : * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256,
46 : : * where F = j/128 for j an integer in [0, 128].
47 : : *
48 : : * log(2^m) = log2_hi*m + log2_tail*m
49 : : * since m is an integer, the dominant term is exact.
50 : : * m has at most 10 digits (for subnormal numbers),
51 : : * and log2_hi has 11 trailing zero bits.
52 : : *
53 : : * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
54 : : * logF_hi[] + 512 is exact.
55 : : *
56 : : * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ...
57 : : * the leading term is calculated to extra precision in two
58 : : * parts, the larger of which adds exactly to the dominant
59 : : * m and F terms.
60 : : * There are two cases:
61 : : * 1. when m, j are non-zero (m | j), use absolute
62 : : * precision for the leading term.
63 : : * 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1).
64 : : * In this case, use a relative precision of 24 bits.
65 : : * (This is done differently in the original paper)
66 : : *
67 : : * Special cases:
68 : : * 0 return signalling -Inf
69 : : * neg return signalling NaN
70 : : * +Inf return +Inf
71 : : */
72 : :
73 : : #define N 128
74 : :
75 : : /* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
76 : : * Used for generation of extend precision logarithms.
77 : : * The constant 35184372088832 is 2^45, so the divide is exact.
78 : : * It ensures correct reading of logF_head, even for inaccurate
79 : : * decimal-to-binary conversion routines. (Everybody gets the
80 : : * right answer for integers less than 2^53.)
81 : : * Values for log(F) were generated using error < 10^-57 absolute
82 : : * with the bc -l package.
83 : : */
84 : : static double A1 = .08333333333333178827;
85 : : static double A2 = .01250000000377174923;
86 : : static double A3 = .002232139987919447809;
87 : : static double A4 = .0004348877777076145742;
88 : :
89 : : static double logF_head[N+1] = {
90 : : 0.,
91 : : .007782140442060381246,
92 : : .015504186535963526694,
93 : : .023167059281547608406,
94 : : .030771658666765233647,
95 : : .038318864302141264488,
96 : : .045809536031242714670,
97 : : .053244514518837604555,
98 : : .060624621816486978786,
99 : : .067950661908525944454,
100 : : .075223421237524235039,
101 : : .082443669210988446138,
102 : : .089612158689760690322,
103 : : .096729626458454731618,
104 : : .103796793681567578460,
105 : : .110814366340264314203,
106 : : .117783035656430001836,
107 : : .124703478501032805070,
108 : : .131576357788617315236,
109 : : .138402322859292326029,
110 : : .145182009844575077295,
111 : : .151916042025732167530,
112 : : .158605030176659056451,
113 : : .165249572895390883786,
114 : : .171850256926518341060,
115 : : .178407657472689606947,
116 : : .184922338493834104156,
117 : : .191394852999565046047,
118 : : .197825743329758552135,
119 : : .204215541428766300668,
120 : : .210564769107350002741,
121 : : .216873938300523150246,
122 : : .223143551314024080056,
123 : : .229374101064877322642,
124 : : .235566071312860003672,
125 : : .241719936886966024758,
126 : : .247836163904594286577,
127 : : .253915209980732470285,
128 : : .259957524436686071567,
129 : : .265963548496984003577,
130 : : .271933715484010463114,
131 : : .277868451003087102435,
132 : : .283768173130738432519,
133 : : .289633292582948342896,
134 : : .295464212893421063199,
135 : : .301261330578199704177,
136 : : .307025035294827830512,
137 : : .312755710004239517729,
138 : : .318453731118097493890,
139 : : .324119468654316733591,
140 : : .329753286372579168528,
141 : : .335355541920762334484,
142 : : .340926586970454081892,
143 : : .346466767346100823488,
144 : : .351976423156884266063,
145 : : .357455888922231679316,
146 : : .362905493689140712376,
147 : : .368325561158599157352,
148 : : .373716409793814818840,
149 : : .379078352934811846353,
150 : : .384411698910298582632,
151 : : .389716751140440464951,
152 : : .394993808240542421117,
153 : : .400243164127459749579,
154 : : .405465108107819105498,
155 : : .410659924985338875558,
156 : : .415827895143593195825,
157 : : .420969294644237379543,
158 : : .426084395310681429691,
159 : : .431173464818130014464,
160 : : .436236766774527495726,
161 : : .441274560805140936281,
162 : : .446287102628048160113,
163 : : .451274644139630254358,
164 : : .456237433481874177232,
165 : : .461175715122408291790,
166 : : .466089729924533457960,
167 : : .470979715219073113985,
168 : : .475845904869856894947,
169 : : .480688529345570714212,
170 : : .485507815781602403149,
171 : : .490303988045525329653,
172 : : .495077266798034543171,
173 : : .499827869556611403822,
174 : : .504556010751912253908,
175 : : .509261901790523552335,
176 : : .513945751101346104405,
177 : : .518607764208354637958,
178 : : .523248143765158602036,
179 : : .527867089620485785417,
180 : : .532464798869114019908,
181 : : .537041465897345915436,
182 : : .541597282432121573947,
183 : : .546132437597407260909,
184 : : .550647117952394182793,
185 : : .555141507540611200965,
186 : : .559615787935399566777,
187 : : .564070138285387656651,
188 : : .568504735352689749561,
189 : : .572919753562018740922,
190 : : .577315365035246941260,
191 : : .581691739635061821900,
192 : : .586049045003164792433,
193 : : .590387446602107957005,
194 : : .594707107746216934174,
195 : : .599008189645246602594,
196 : : .603290851438941899687,
197 : : .607555250224322662688,
198 : : .611801541106615331955,
199 : : .616029877215623855590,
200 : : .620240409751204424537,
201 : : .624433288012369303032,
202 : : .628608659422752680256,
203 : : .632766669570628437213,
204 : : .636907462236194987781,
205 : : .641031179420679109171,
206 : : .645137961373620782978,
207 : : .649227946625615004450,
208 : : .653301272011958644725,
209 : : .657358072709030238911,
210 : : .661398482245203922502,
211 : : .665422632544505177065,
212 : : .669430653942981734871,
213 : : .673422675212350441142,
214 : : .677398823590920073911,
215 : : .681359224807238206267,
216 : : .685304003098281100392,
217 : : .689233281238557538017,
218 : : .693147180560117703862
219 : : };
220 : :
221 : : static double logF_tail[N+1] = {
222 : : 0.,
223 : : -.00000000000000543229938420049,
224 : : .00000000000000172745674997061,
225 : : -.00000000000001323017818229233,
226 : : -.00000000000001154527628289872,
227 : : -.00000000000000466529469958300,
228 : : .00000000000005148849572685810,
229 : : -.00000000000002532168943117445,
230 : : -.00000000000005213620639136504,
231 : : -.00000000000001819506003016881,
232 : : .00000000000006329065958724544,
233 : : .00000000000008614512936087814,
234 : : -.00000000000007355770219435028,
235 : : .00000000000009638067658552277,
236 : : .00000000000007598636597194141,
237 : : .00000000000002579999128306990,
238 : : -.00000000000004654729747598444,
239 : : -.00000000000007556920687451336,
240 : : .00000000000010195735223708472,
241 : : -.00000000000017319034406422306,
242 : : -.00000000000007718001336828098,
243 : : .00000000000010980754099855238,
244 : : -.00000000000002047235780046195,
245 : : -.00000000000008372091099235912,
246 : : .00000000000014088127937111135,
247 : : .00000000000012869017157588257,
248 : : .00000000000017788850778198106,
249 : : .00000000000006440856150696891,
250 : : .00000000000016132822667240822,
251 : : -.00000000000007540916511956188,
252 : : -.00000000000000036507188831790,
253 : : .00000000000009120937249914984,
254 : : .00000000000018567570959796010,
255 : : -.00000000000003149265065191483,
256 : : -.00000000000009309459495196889,
257 : : .00000000000017914338601329117,
258 : : -.00000000000001302979717330866,
259 : : .00000000000023097385217586939,
260 : : .00000000000023999540484211737,
261 : : .00000000000015393776174455408,
262 : : -.00000000000036870428315837678,
263 : : .00000000000036920375082080089,
264 : : -.00000000000009383417223663699,
265 : : .00000000000009433398189512690,
266 : : .00000000000041481318704258568,
267 : : -.00000000000003792316480209314,
268 : : .00000000000008403156304792424,
269 : : -.00000000000034262934348285429,
270 : : .00000000000043712191957429145,
271 : : -.00000000000010475750058776541,
272 : : -.00000000000011118671389559323,
273 : : .00000000000037549577257259853,
274 : : .00000000000013912841212197565,
275 : : .00000000000010775743037572640,
276 : : .00000000000029391859187648000,
277 : : -.00000000000042790509060060774,
278 : : .00000000000022774076114039555,
279 : : .00000000000010849569622967912,
280 : : -.00000000000023073801945705758,
281 : : .00000000000015761203773969435,
282 : : .00000000000003345710269544082,
283 : : -.00000000000041525158063436123,
284 : : .00000000000032655698896907146,
285 : : -.00000000000044704265010452446,
286 : : .00000000000034527647952039772,
287 : : -.00000000000007048962392109746,
288 : : .00000000000011776978751369214,
289 : : -.00000000000010774341461609578,
290 : : .00000000000021863343293215910,
291 : : .00000000000024132639491333131,
292 : : .00000000000039057462209830700,
293 : : -.00000000000026570679203560751,
294 : : .00000000000037135141919592021,
295 : : -.00000000000017166921336082431,
296 : : -.00000000000028658285157914353,
297 : : -.00000000000023812542263446809,
298 : : .00000000000006576659768580062,
299 : : -.00000000000028210143846181267,
300 : : .00000000000010701931762114254,
301 : : .00000000000018119346366441110,
302 : : .00000000000009840465278232627,
303 : : -.00000000000033149150282752542,
304 : : -.00000000000018302857356041668,
305 : : -.00000000000016207400156744949,
306 : : .00000000000048303314949553201,
307 : : -.00000000000071560553172382115,
308 : : .00000000000088821239518571855,
309 : : -.00000000000030900580513238244,
310 : : -.00000000000061076551972851496,
311 : : .00000000000035659969663347830,
312 : : .00000000000035782396591276383,
313 : : -.00000000000046226087001544578,
314 : : .00000000000062279762917225156,
315 : : .00000000000072838947272065741,
316 : : .00000000000026809646615211673,
317 : : -.00000000000010960825046059278,
318 : : .00000000000002311949383800537,
319 : : -.00000000000058469058005299247,
320 : : -.00000000000002103748251144494,
321 : : -.00000000000023323182945587408,
322 : : -.00000000000042333694288141916,
323 : : -.00000000000043933937969737844,
324 : : .00000000000041341647073835565,
325 : : .00000000000006841763641591466,
326 : : .00000000000047585534004430641,
327 : : .00000000000083679678674757695,
328 : : -.00000000000085763734646658640,
329 : : .00000000000021913281229340092,
330 : : -.00000000000062242842536431148,
331 : : -.00000000000010983594325438430,
332 : : .00000000000065310431377633651,
333 : : -.00000000000047580199021710769,
334 : : -.00000000000037854251265457040,
335 : : .00000000000040939233218678664,
336 : : .00000000000087424383914858291,
337 : : .00000000000025218188456842882,
338 : : -.00000000000003608131360422557,
339 : : -.00000000000050518555924280902,
340 : : .00000000000078699403323355317,
341 : : -.00000000000067020876961949060,
342 : : .00000000000016108575753932458,
343 : : .00000000000058527188436251509,
344 : : -.00000000000035246757297904791,
345 : : -.00000000000018372084495629058,
346 : : .00000000000088606689813494916,
347 : : .00000000000066486268071468700,
348 : : .00000000000063831615170646519,
349 : : .00000000000025144230728376072,
350 : : -.00000000000017239444525614834
351 : : };
352 : :
353 : : #if 0
354 : : OLM_DLLEXPORT double
355 : : #ifdef _ANSI_SOURCE
356 : : log(double x)
357 : : #else
358 : : log(x) double x;
359 : : #endif
360 : : {
361 : : int m, j;
362 : : double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
363 : : volatile double u1;
364 : :
365 : : /* Catch special cases */
366 : : if (x <= 0)
367 : : if (x == zero) /* log(0) = -Inf */
368 : : return (-one/zero);
369 : : else /* log(neg) = NaN */
370 : : return (zero/zero);
371 : : else if (!finite(x))
372 : : return (x+x); /* x = NaN, Inf */
373 : :
374 : : /* Argument reduction: 1 <= g < 2; x/2^m = g; */
375 : : /* y = F*(1 + f/F) for |f| <= 2^-8 */
376 : :
377 : : m = logb(x);
378 : : g = ldexp(x, -m);
379 : : if (m == -1022) {
380 : : j = logb(g), m += j;
381 : : g = ldexp(g, -j);
382 : : }
383 : : j = N*(g-1) + .5;
384 : : F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */
385 : : f = g - F;
386 : :
387 : : /* Approximate expansion for log(1+f/F) ~= u + q */
388 : : g = 1/(2*F+f);
389 : : u = 2*f*g;
390 : : v = u*u;
391 : : q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
392 : :
393 : : /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8,
394 : : * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits.
395 : : * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750
396 : : */
397 : : if (m | j)
398 : : u1 = u + 513, u1 -= 513;
399 : :
400 : : /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero;
401 : : * u1 = u to 24 bits.
402 : : */
403 : : else
404 : : u1 = u, TRUNC(u1);
405 : : u2 = (2.0*(f - F*u1) - u1*f) * g;
406 : : /* u1 + u2 = 2f/(2F+f) to extra precision. */
407 : :
408 : : /* log(x) = log(2^m*F*(1+f/F)) = */
409 : : /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */
410 : : /* (exact) + (tiny) */
411 : :
412 : : u1 += m*logF_head[N] + logF_head[j]; /* exact */
413 : : u2 = (u2 + logF_tail[j]) + q; /* tiny */
414 : : u2 += logF_tail[N]*m;
415 : : return (u1 + u2);
416 : : }
417 : : #endif
418 : :
419 : : /*
420 : : * Extra precision variant, returning struct {double a, b;};
421 : : * log(x) = a+b to 63 bits, with a rounded to 26 bits.
422 : : */
423 : : struct Double
424 : : #ifdef _ANSI_SOURCE
425 : : __log__D(double x)
426 : : #else
427 : 0 : __log__D(x) double x;
428 : : #endif
429 : : {
430 : : int m, j;
431 : : double F, f, g, q, u, v, u2;
432 : : volatile double u1;
433 : : struct Double r;
434 : :
435 : : /* Argument reduction: 1 <= g < 2; x/2^m = g; */
436 : : /* y = F*(1 + f/F) for |f| <= 2^-8 */
437 : :
438 : 0 : m = logb(x);
439 : 0 : g = ldexp(x, -m);
440 [ # # ]: 0 : if (m == -1022) {
441 : 0 : j = logb(g), m += j;
442 : 0 : g = ldexp(g, -j);
443 : : }
444 : 0 : j = N*(g-1) + .5;
445 : 0 : F = (1.0/N) * j + 1;
446 : 0 : f = g - F;
447 : :
448 : 0 : g = 1/(2*F+f);
449 : 0 : u = 2*f*g;
450 : 0 : v = u*u;
451 : 0 : q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
452 [ # # ]: 0 : if (m | j)
453 : 0 : u1 = u + 513, u1 -= 513;
454 : : else
455 : 0 : u1 = u, TRUNC(u1);
456 : 0 : u2 = (2.0*(f - F*u1) - u1*f) * g;
457 : :
458 : 0 : u1 += m*logF_head[N] + logF_head[j];
459 : :
460 : 0 : u2 += logF_tail[j]; u2 += q;
461 : 0 : u2 += logF_tail[N]*m;
462 : 0 : r.a = u1 + u2; /* Only difference is here */
463 : 0 : TRUNC(r.a);
464 : 0 : r.b = (u1 - r.a) + u2;
465 : 0 : return (r);
466 : : }
|
