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.1 (Berkeley) 06/04/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 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 double logb(), ldexp(); 352 volatile double u1; 353 354 /* Catch special cases */ 355 if (x <= 0) 356 if (_IEEE && x == zero) /* log(0) = -Inf */ 357 return (-one/zero); 358 else if (_IEEE) /* log(neg) = NaN */ 359 return (zero/zero); 360 else if (x == zero) /* NOT REACHED IF _IEEE */ 361 return (infnan(-ERANGE)); 362 else 363 return (infnan(EDOM)); 364 else if (!finite(x)) 365 if (_IEEE) /* x = NaN, Inf */ 366 return (x+x); 367 else 368 return (infnan(ERANGE)); 369 370 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 371 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 372 373 m = logb(x); 374 g = ldexp(x, -m); 375 if (_IEEE && m == -1022) { 376 j = logb(g), m += j; 377 g = ldexp(g, -j); 378 } 379 j = N*(g-1) + .5; 380 F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */ 381 f = g - F; 382 383 /* Approximate expansion for log(1+f/F) ~= u + q */ 384 g = 1/(2*F+f); 385 u = 2*f*g; 386 v = u*u; 387 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 388 389 /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, 390 * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits. 391 * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750 392 */ 393 if (m | j) 394 u1 = u + 513, u1 -= 513; 395 396 /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero; 397 * u1 = u to 24 bits. 398 */ 399 else 400 u1 = u, TRUNC(u1); 401 u2 = (2.0*(f - F*u1) - u1*f) * g; 402 /* u1 + u2 = 2f/(2F+f) to extra precision. */ 403 404 /* log(x) = log(2^m*F*(1+f/F)) = */ 405 /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */ 406 /* (exact) + (tiny) */ 407 408 u1 += m*logF_head[N] + logF_head[j]; /* exact */ 409 u2 = (u2 + logF_tail[j]) + q; /* tiny */ 410 u2 += logF_tail[N]*m; 411 return (u1 + u2); 412 } 413 414 /* 415 * Extra precision variant, returning struct {double a, b;}; 416 * log(x) = a+b to 63 bits, with a is rounded to 26 bits. 417 */ 418 struct Double 419 #ifdef _ANSI_SOURCE 420 __log__D(double x) 421 #else 422 __log__D(x) double x; 423 #endif 424 { 425 int m, j; 426 double F, f, g, q, u, v, u2, one = 1.0; 427 double logb(), ldexp(); 428 volatile double u1; 429 struct Double r; 430 431 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 432 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 433 434 m = logb(x); 435 g = ldexp(x, -m); 436 if (_IEEE && m == -1022) { 437 j = logb(g), m += j; 438 g = ldexp(g, -j); 439 } 440 j = N*(g-1) + .5; 441 F = (1.0/N) * j + 1; 442 f = g - F; 443 444 g = 1/(2*F+f); 445 u = 2*f*g; 446 v = u*u; 447 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 448 if (m | j) 449 u1 = u + 513, u1 -= 513; 450 else 451 u1 = u, TRUNC(u1); 452 u2 = (2.0*(f - F*u1) - u1*f) * g; 453 454 u1 += m*logF_head[N] + logF_head[j]; 455 456 u2 += logF_tail[j]; u2 += q; 457 u2 += logF_tail[N]*m; 458 r.a = u1 + u2; /* Only difference is here */ 459 TRUNC(r.a); 460 r.b = (u1 - r.a) + u2; 461 return (r); 462 } 463