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