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 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 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