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