1 #include "fconv.h" 2 3 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines (dmg). 4 * 5 * This strtod returns a nearest machine number to the input decimal 6 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 7 * broken by the IEEE round-even rule. Otherwise ties are broken by 8 * biased rounding (add half and chop). 9 * 10 * Inspired loosely by William D. Clinger's paper "How to Read Floating 11 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 12 * 13 * Modifications: 14 * 15 * 1. We only require IEEE, IBM, or VAX double-precision 16 * arithmetic (not IEEE double-extended). 17 * 2. We get by with floating-point arithmetic in a case that 18 * Clinger missed -- when we're computing d * 10^n 19 * for a small integer d and the integer n is not too 20 * much larger than 22 (the maximum integer k for which 21 * we can represent 10^k exactly), we may be able to 22 * compute (d*10^k) * 10^(e-k) with just one roundoff. 23 * 3. Rather than a bit-at-a-time adjustment of the binary 24 * result in the hard case, we use floating-point 25 * arithmetic to determine the adjustment to within 26 * one bit; only in really hard cases do we need to 27 * compute a second residual. 28 * 4. Because of 3., we don't need a large table of powers of 10 29 * for ten-to-e (just some small tables, e.g. of 10^k 30 * for 0 <= k <= 22). 31 */ 32 33 #ifdef RND_PRODQUOT 34 #define rounded_product(a,b) a = rnd_prod(a, b) 35 #define rounded_quotient(a,b) a = rnd_quot(a, b) 36 extern double rnd_prod(double, double), rnd_quot(double, double); 37 #else 38 #define rounded_product(a,b) a *= b 39 #define rounded_quotient(a,b) a /= b 40 #endif 41 42 static double 43 ulp(double xarg) 44 { 45 register long L; 46 Dul a; 47 Dul x; 48 49 x.d = xarg; 50 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 51 #ifndef Sudden_Underflow 52 if (L > 0) { 53 #endif 54 #ifdef IBM 55 L |= Exp_msk1 >> 4; 56 #endif 57 word0(a) = L; 58 word1(a) = 0; 59 #ifndef Sudden_Underflow 60 } 61 else { 62 L = -L >> Exp_shift; 63 if (L < Exp_shift) { 64 word0(a) = 0x80000 >> L; 65 word1(a) = 0; 66 } 67 else { 68 word0(a) = 0; 69 L -= Exp_shift; 70 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 71 } 72 } 73 #endif 74 return a.d; 75 } 76 77 static Bigint * 78 s2b(CONST char *s, int nd0, int nd, unsigned long y9) 79 { 80 Bigint *b; 81 int i, k; 82 long x, y; 83 84 x = (nd + 8) / 9; 85 for(k = 0, y = 1; x > y; y <<= 1, k++) ; 86 #ifdef Pack_32 87 b = Balloc(k); 88 b->x[0] = y9; 89 b->wds = 1; 90 #else 91 b = Balloc(k+1); 92 b->x[0] = y9 & 0xffff; 93 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; 94 #endif 95 96 i = 9; 97 if (9 < nd0) { 98 s += 9; 99 do b = multadd(b, 10, *s++ - '0'); 100 while(++i < nd0); 101 s++; 102 } 103 else 104 s += 10; 105 for(; i < nd; i++) 106 b = multadd(b, 10, *s++ - '0'); 107 return b; 108 } 109 110 static double 111 b2d(Bigint *a, int *e) 112 { 113 unsigned long *xa, *xa0, w, y, z; 114 int k; 115 Dul d; 116 #ifdef VAX 117 unsigned long d0, d1; 118 #else 119 #define d0 word0(d) 120 #define d1 word1(d) 121 #endif 122 123 xa0 = a->x; 124 xa = xa0 + a->wds; 125 y = *--xa; 126 #ifdef DEBUG 127 if (!y) Bug("zero y in b2d"); 128 #endif 129 k = hi0bits(y); 130 *e = 32 - k; 131 #ifdef Pack_32 132 if (k < Ebits) { 133 d0 = Exp_1 | y >> Ebits - k; 134 w = xa > xa0 ? *--xa : 0; 135 d1 = y << (32-Ebits) + k | w >> Ebits - k; 136 goto ret_d; 137 } 138 z = xa > xa0 ? *--xa : 0; 139 if (k -= Ebits) { 140 d0 = Exp_1 | y << k | z >> 32 - k; 141 y = xa > xa0 ? *--xa : 0; 142 d1 = z << k | y >> 32 - k; 143 } 144 else { 145 d0 = Exp_1 | y; 146 d1 = z; 147 } 148 #else 149 if (k < Ebits + 16) { 150 z = xa > xa0 ? *--xa : 0; 151 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 152 w = xa > xa0 ? *--xa : 0; 153 y = xa > xa0 ? *--xa : 0; 154 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 155 goto ret_d; 156 } 157 z = xa > xa0 ? *--xa : 0; 158 w = xa > xa0 ? *--xa : 0; 159 k -= Ebits + 16; 160 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 161 y = xa > xa0 ? *--xa : 0; 162 d1 = w << k + 16 | y << k; 163 #endif 164 ret_d: 165 #ifdef VAX 166 word0(d) = d0 >> 16 | d0 << 16; 167 word1(d) = d1 >> 16 | d1 << 16; 168 #else 169 #undef d0 170 #undef d1 171 #endif 172 return d.d; 173 } 174 175 static double 176 ratio(Bigint *a, Bigint *b) 177 { 178 Dul da, db; 179 int k, ka, kb; 180 181 da.d = b2d(a, &ka); 182 db.d = b2d(b, &kb); 183 #ifdef Pack_32 184 k = ka - kb + 32*(a->wds - b->wds); 185 #else 186 k = ka - kb + 16*(a->wds - b->wds); 187 #endif 188 #ifdef IBM 189 if (k > 0) { 190 word0(da) += (k >> 2)*Exp_msk1; 191 if (k &= 3) 192 da *= 1 << k; 193 } 194 else { 195 k = -k; 196 word0(db) += (k >> 2)*Exp_msk1; 197 if (k &= 3) 198 db *= 1 << k; 199 } 200 #else 201 if (k > 0) 202 word0(da) += k*Exp_msk1; 203 else { 204 k = -k; 205 word0(db) += k*Exp_msk1; 206 } 207 #endif 208 return da.d / db.d; 209 } 210 211 double 212 strtod(CONST char *s00, char **se) 213 { 214 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 215 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 216 CONST char *s, *s0, *s1; 217 double aadj, aadj1, adj; 218 Dul rv, rv0; 219 long L; 220 unsigned long y, z; 221 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 222 sign = nz0 = nz = 0; 223 rv.d = 0.; 224 for(s = s00;;s++) switch(*s) { 225 case '-': 226 sign = 1; 227 /* no break */ 228 case '+': 229 if (*++s) 230 goto break2; 231 /* no break */ 232 case 0: 233 s = s00; 234 goto ret; 235 case '\t': 236 case '\n': 237 case '\v': 238 case '\f': 239 case '\r': 240 case ' ': 241 continue; 242 default: 243 goto break2; 244 } 245 break2: 246 if (*s == '0') { 247 nz0 = 1; 248 while(*++s == '0') ; 249 if (!*s) 250 goto ret; 251 } 252 s0 = s; 253 y = z = 0; 254 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 255 if (nd < 9) 256 y = 10*y + c - '0'; 257 else if (nd < 16) 258 z = 10*z + c - '0'; 259 nd0 = nd; 260 if (c == '.') { 261 c = *++s; 262 if (!nd) { 263 for(; c == '0'; c = *++s) 264 nz++; 265 if (c > '0' && c <= '9') { 266 s0 = s; 267 nf += nz; 268 nz = 0; 269 goto have_dig; 270 } 271 goto dig_done; 272 } 273 for(; c >= '0' && c <= '9'; c = *++s) { 274 have_dig: 275 nz++; 276 if (c -= '0') { 277 nf += nz; 278 for(i = 1; i < nz; i++) 279 if (nd++ < 9) 280 y *= 10; 281 else if (nd <= DBL_DIG + 1) 282 z *= 10; 283 if (nd++ < 9) 284 y = 10*y + c; 285 else if (nd <= DBL_DIG + 1) 286 z = 10*z + c; 287 nz = 0; 288 } 289 } 290 } 291 dig_done: 292 e = 0; 293 if (c == 'e' || c == 'E') { 294 if (!nd && !nz && !nz0) { 295 s = s00; 296 goto ret; 297 } 298 s00 = s; 299 esign = 0; 300 switch(c = *++s) { 301 case '-': 302 esign = 1; 303 case '+': 304 c = *++s; 305 } 306 if (c >= '0' && c <= '9') { 307 while(c == '0') 308 c = *++s; 309 if (c > '0' && c <= '9') { 310 e = c - '0'; 311 s1 = s; 312 while((c = *++s) >= '0' && c <= '9') 313 e = 10*e + c - '0'; 314 if (s - s1 > 8) 315 /* Avoid confusion from exponents 316 * so large that e might overflow. 317 */ 318 e = 9999999; 319 if (esign) 320 e = -e; 321 } 322 else 323 e = 0; 324 } 325 else 326 s = s00; 327 } 328 if (!nd) { 329 if (!nz && !nz0) 330 s = s00; 331 goto ret; 332 } 333 e1 = e -= nf; 334 335 /* Now we have nd0 digits, starting at s0, followed by a 336 * decimal point, followed by nd-nd0 digits. The number we're 337 * after is the integer represented by those digits times 338 * 10**e */ 339 340 if (!nd0) 341 nd0 = nd; 342 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 343 rv.d = y; 344 if (k > 9) 345 rv.d = tens[k - 9] * rv.d + z; 346 if (nd <= DBL_DIG 347 #ifndef RND_PRODQUOT 348 && FLT_ROUNDS == 1 349 #endif 350 ) { 351 if (!e) 352 goto ret; 353 if (e > 0) { 354 if (e <= Ten_pmax) { 355 #ifdef VAX 356 goto vax_ovfl_check; 357 #else 358 /* rv = */ rounded_product(rv.d, tens[e]); 359 goto ret; 360 #endif 361 } 362 i = DBL_DIG - nd; 363 if (e <= Ten_pmax + i) { 364 /* A fancier test would sometimes let us do 365 * this for larger i values. 366 */ 367 e -= i; 368 rv.d *= tens[i]; 369 #ifdef VAX 370 /* VAX exponent range is so narrow we must 371 * worry about overflow here... 372 */ 373 vax_ovfl_check: 374 word0(rv) -= P*Exp_msk1; 375 /* rv = */ rounded_product(rv.d, tens[e]); 376 if ((word0(rv) & Exp_mask) 377 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 378 goto ovfl; 379 word0(rv) += P*Exp_msk1; 380 #else 381 /* rv = */ rounded_product(rv.d, tens[e]); 382 #endif 383 goto ret; 384 } 385 } 386 else if (e >= -Ten_pmax) { 387 /* rv = */ rounded_quotient(rv.d, tens[-e]); 388 goto ret; 389 } 390 } 391 e1 += nd - k; 392 393 /* Get starting approximation = rv * 10**e1 */ 394 395 if (e1 > 0) { 396 if (i = e1 & 15) 397 rv.d *= tens[i]; 398 if (e1 &= ~15) { 399 if (e1 > DBL_MAX_10_EXP) { 400 ovfl: 401 errno = ERANGE; 402 rv.d = HUGE_VAL; 403 goto ret; 404 } 405 if (e1 >>= 4) { 406 for(j = 0; e1 > 1; j++, e1 >>= 1) 407 if (e1 & 1) 408 rv.d *= bigtens[j]; 409 /* The last multiplication could overflow. */ 410 word0(rv) -= P*Exp_msk1; 411 rv.d *= bigtens[j]; 412 if ((z = word0(rv) & Exp_mask) 413 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 414 goto ovfl; 415 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 416 /* set to largest number */ 417 /* (Can't trust DBL_MAX) */ 418 word0(rv) = Big0; 419 word1(rv) = Big1; 420 } 421 else 422 word0(rv) += P*Exp_msk1; 423 } 424 425 } 426 } 427 else if (e1 < 0) { 428 e1 = -e1; 429 if (i = e1 & 15) 430 rv.d /= tens[i]; 431 if (e1 &= ~15) { 432 e1 >>= 4; 433 for(j = 0; e1 > 1; j++, e1 >>= 1) 434 if (e1 & 1) 435 rv.d *= tinytens[j]; 436 /* The last multiplication could underflow. */ 437 rv0.d = rv.d; 438 rv.d *= tinytens[j]; 439 if (rv.d == 0) { 440 rv.d = 2.*rv0.d; 441 rv.d *= tinytens[j]; 442 if (rv.d == 0) { 443 undfl: 444 rv.d = 0.; 445 errno = ERANGE; 446 goto ret; 447 } 448 word0(rv) = Tiny0; 449 word1(rv) = Tiny1; 450 /* The refinement below will clean 451 * this approximation up. 452 */ 453 } 454 } 455 } 456 457 /* Now the hard part -- adjusting rv to the correct value.*/ 458 459 /* Put digits into bd: true value = bd * 10^e */ 460 461 bd0 = s2b(s0, nd0, nd, y); 462 463 for(;;) { 464 bd = Balloc(bd0->k); 465 Bcopy(bd, bd0); 466 bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */ 467 bs = i2b(1); 468 469 if (e >= 0) { 470 bb2 = bb5 = 0; 471 bd2 = bd5 = e; 472 } 473 else { 474 bb2 = bb5 = -e; 475 bd2 = bd5 = 0; 476 } 477 if (bbe >= 0) 478 bb2 += bbe; 479 else 480 bd2 -= bbe; 481 bs2 = bb2; 482 #ifdef Sudden_Underflow 483 #ifdef IBM 484 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 485 #else 486 j = P + 1 - bbbits; 487 #endif 488 #else 489 i = bbe + bbbits - 1; /* logb(rv) */ 490 if (i < Emin) /* denormal */ 491 j = bbe + (P-Emin); 492 else 493 j = P + 1 - bbbits; 494 #endif 495 bb2 += j; 496 bd2 += j; 497 i = bb2 < bd2 ? bb2 : bd2; 498 if (i > bs2) 499 i = bs2; 500 if (i > 0) { 501 bb2 -= i; 502 bd2 -= i; 503 bs2 -= i; 504 } 505 if (bb5 > 0) { 506 bs = pow5mult(bs, bb5); 507 bb1 = mult(bs, bb); 508 Bfree(bb); 509 bb = bb1; 510 } 511 if (bb2 > 0) 512 bb = lshift(bb, bb2); 513 if (bd5 > 0) 514 bd = pow5mult(bd, bd5); 515 if (bd2 > 0) 516 bd = lshift(bd, bd2); 517 if (bs2 > 0) 518 bs = lshift(bs, bs2); 519 delta = diff(bb, bd); 520 dsign = delta->sign; 521 delta->sign = 0; 522 i = cmp(delta, bs); 523 if (i < 0) { 524 /* Error is less than half an ulp -- check for 525 * special case of mantissa a power of two. 526 */ 527 if (dsign || word1(rv) || word0(rv) & Bndry_mask) 528 break; 529 delta = lshift(delta,Log2P); 530 if (cmp(delta, bs) > 0) 531 goto drop_down; 532 break; 533 } 534 if (i == 0) { 535 /* exactly half-way between */ 536 if (dsign) { 537 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 538 && word1(rv) == 0xffffffff) { 539 /*boundary case -- increment exponent*/ 540 word0(rv) = (word0(rv) & Exp_mask) 541 + Exp_msk1 542 #ifdef IBM 543 | Exp_msk1 >> 4 544 #endif 545 ; 546 word1(rv) = 0; 547 break; 548 } 549 } 550 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 551 drop_down: 552 /* boundary case -- decrement exponent */ 553 #ifdef Sudden_Underflow 554 L = word0(rv) & Exp_mask; 555 #ifdef IBM 556 if (L < Exp_msk1) 557 #else 558 if (L <= Exp_msk1) 559 #endif 560 goto undfl; 561 L -= Exp_msk1; 562 #else 563 L = (word0(rv) & Exp_mask) - Exp_msk1; 564 #endif 565 word0(rv) = L | Bndry_mask1; 566 word1(rv) = 0xffffffff; 567 #ifdef IBM 568 goto cont; 569 #else 570 break; 571 #endif 572 } 573 #ifndef ROUND_BIASED 574 if (!(word1(rv) & LSB)) 575 break; 576 #endif 577 if (dsign) 578 rv.d += ulp(rv.d); 579 #ifndef ROUND_BIASED 580 else { 581 rv.d -= ulp(rv.d); 582 #ifndef Sudden_Underflow 583 if (rv.d == 0) 584 goto undfl; 585 #endif 586 } 587 #endif 588 break; 589 } 590 if ((aadj = ratio(delta, bs)) <= 2.) { 591 if (dsign) 592 aadj = aadj1 = 1.; 593 else if (word1(rv) || word0(rv) & Bndry_mask) { 594 #ifndef Sudden_Underflow 595 if (word1(rv) == Tiny1 && !word0(rv)) 596 goto undfl; 597 #endif 598 aadj = 1.; 599 aadj1 = -1.; 600 } 601 else { 602 /* special case -- power of FLT_RADIX to be */ 603 /* rounded down... */ 604 605 if (aadj < 2./FLT_RADIX) 606 aadj = 1./FLT_RADIX; 607 else 608 aadj *= 0.5; 609 aadj1 = -aadj; 610 } 611 } 612 else { 613 aadj *= 0.5; 614 aadj1 = dsign ? aadj : -aadj; 615 #ifdef Check_FLT_ROUNDS 616 switch(FLT_ROUNDS) { 617 case 2: /* towards +infinity */ 618 aadj1 -= 0.5; 619 break; 620 case 0: /* towards 0 */ 621 case 3: /* towards -infinity */ 622 aadj1 += 0.5; 623 } 624 #else 625 if (FLT_ROUNDS == 0) 626 aadj1 += 0.5; 627 #endif 628 } 629 y = word0(rv) & Exp_mask; 630 631 /* Check for overflow */ 632 633 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 634 rv0.d = rv.d; 635 word0(rv) -= P*Exp_msk1; 636 adj = aadj1 * ulp(rv.d); 637 rv.d += adj; 638 if ((word0(rv) & Exp_mask) >= 639 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 640 if (word0(rv0) == Big0 && word1(rv0) == Big1) 641 goto ovfl; 642 word0(rv) = Big0; 643 word1(rv) = Big1; 644 goto cont; 645 } 646 else 647 word0(rv) += P*Exp_msk1; 648 } 649 else { 650 #ifdef Sudden_Underflow 651 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 652 rv0.d = rv.d; 653 word0(rv) += P*Exp_msk1; 654 adj = aadj1 * ulp(rv.d); 655 rv.d += adj; 656 #ifdef IBM 657 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 658 #else 659 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 660 #endif 661 { 662 if (word0(rv0) == Tiny0 663 && word1(rv0) == Tiny1) 664 goto undfl; 665 word0(rv) = Tiny0; 666 word1(rv) = Tiny1; 667 goto cont; 668 } 669 else 670 word0(rv) -= P*Exp_msk1; 671 } 672 else { 673 adj = aadj1 * ulp(rv.d); 674 rv.d += adj; 675 } 676 #else 677 /* Compute adj so that the IEEE rounding rules will 678 * correctly round rv + adj in some half-way cases. 679 * If rv * ulp(rv) is denormalized (i.e., 680 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 681 * trouble from bits lost to denormalization; 682 * example: 1.2e-307 . 683 */ 684 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { 685 aadj1 = (double)(int)(aadj + 0.5); 686 if (!dsign) 687 aadj1 = -aadj1; 688 } 689 adj = aadj1 * ulp(rv.d); 690 rv.d += adj; 691 #endif 692 } 693 z = word0(rv) & Exp_mask; 694 if (y == z) { 695 /* Can we stop now? */ 696 L = aadj; 697 aadj -= L; 698 /* The tolerances below are conservative. */ 699 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 700 if (aadj < .4999999 || aadj > .5000001) 701 break; 702 } 703 else if (aadj < .4999999/FLT_RADIX) 704 break; 705 } 706 cont: 707 Bfree(bb); 708 Bfree(bd); 709 Bfree(bs); 710 Bfree(delta); 711 } 712 Bfree(bb); 713 Bfree(bd); 714 Bfree(bs); 715 Bfree(bd0); 716 Bfree(delta); 717 ret: 718 if (se) 719 *se = (char *)s; 720 return sign ? -rv.d : rv.d; 721 } 722