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