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 bd0 = 0; 347 if (nd <= DBL_DIG 348 #ifndef RND_PRODQUOT 349 && FLT_ROUNDS == 1 350 #endif 351 ) { 352 if (!e) 353 goto ret; 354 if (e > 0) { 355 if (e <= Ten_pmax) { 356 #ifdef VAX 357 goto vax_ovfl_check; 358 #else 359 /* rv = */ rounded_product(rv.d, tens[e]); 360 goto ret; 361 #endif 362 } 363 i = DBL_DIG - nd; 364 if (e <= Ten_pmax + i) { 365 /* A fancier test would sometimes let us do 366 * this for larger i values. 367 */ 368 e -= i; 369 rv.d *= tens[i]; 370 #ifdef VAX 371 /* VAX exponent range is so narrow we must 372 * worry about overflow here... 373 */ 374 vax_ovfl_check: 375 word0(rv) -= P*Exp_msk1; 376 /* rv = */ rounded_product(rv.d, tens[e]); 377 if ((word0(rv) & Exp_mask) 378 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 379 goto ovfl; 380 word0(rv) += P*Exp_msk1; 381 #else 382 /* rv = */ rounded_product(rv.d, tens[e]); 383 #endif 384 goto ret; 385 } 386 } 387 else if (e >= -Ten_pmax) { 388 /* rv = */ rounded_quotient(rv.d, tens[-e]); 389 goto ret; 390 } 391 } 392 e1 += nd - k; 393 394 /* Get starting approximation = rv * 10**e1 */ 395 396 if (e1 > 0) { 397 if (i = e1 & 15) 398 rv.d *= tens[i]; 399 if (e1 &= ~15) { 400 if (e1 > DBL_MAX_10_EXP) { 401 ovfl: 402 errno = ERANGE; 403 rv.d = HUGE_VAL; 404 if (bd0) 405 goto retfree; 406 goto ret; 407 } 408 if (e1 >>= 4) { 409 for(j = 0; e1 > 1; j++, e1 >>= 1) 410 if (e1 & 1) 411 rv.d *= bigtens[j]; 412 /* The last multiplication could overflow. */ 413 word0(rv) -= P*Exp_msk1; 414 rv.d *= bigtens[j]; 415 if ((z = word0(rv) & Exp_mask) 416 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 417 goto ovfl; 418 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 419 /* set to largest number */ 420 /* (Can't trust DBL_MAX) */ 421 word0(rv) = Big0; 422 word1(rv) = Big1; 423 } 424 else 425 word0(rv) += P*Exp_msk1; 426 } 427 428 } 429 } 430 else if (e1 < 0) { 431 e1 = -e1; 432 if (i = e1 & 15) 433 rv.d /= tens[i]; 434 if (e1 &= ~15) { 435 e1 >>= 4; 436 if (e1 >= 1 << n_bigtens) 437 goto undfl; 438 for(j = 0; e1 > 1; j++, e1 >>= 1) 439 if (e1 & 1) 440 rv.d *= tinytens[j]; 441 /* The last multiplication could underflow. */ 442 rv0.d = rv.d; 443 rv.d *= tinytens[j]; 444 if (rv.d == 0) { 445 rv.d = 2.*rv0.d; 446 rv.d *= tinytens[j]; 447 if (rv.d == 0) { 448 undfl: 449 rv.d = 0.; 450 errno = ERANGE; 451 if (bd0) 452 goto retfree; 453 goto ret; 454 } 455 word0(rv) = Tiny0; 456 word1(rv) = Tiny1; 457 /* The refinement below will clean 458 * this approximation up. 459 */ 460 } 461 } 462 } 463 464 /* Now the hard part -- adjusting rv to the correct value.*/ 465 466 /* Put digits into bd: true value = bd * 10^e */ 467 468 bd0 = s2b(s0, nd0, nd, y); 469 470 for(;;) { 471 bd = Balloc(bd0->k); 472 Bcopy(bd, bd0); 473 bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */ 474 bs = i2b(1); 475 476 if (e >= 0) { 477 bb2 = bb5 = 0; 478 bd2 = bd5 = e; 479 } 480 else { 481 bb2 = bb5 = -e; 482 bd2 = bd5 = 0; 483 } 484 if (bbe >= 0) 485 bb2 += bbe; 486 else 487 bd2 -= bbe; 488 bs2 = bb2; 489 #ifdef Sudden_Underflow 490 #ifdef IBM 491 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 492 #else 493 j = P + 1 - bbbits; 494 #endif 495 #else 496 i = bbe + bbbits - 1; /* logb(rv) */ 497 if (i < Emin) /* denormal */ 498 j = bbe + (P-Emin); 499 else 500 j = P + 1 - bbbits; 501 #endif 502 bb2 += j; 503 bd2 += j; 504 i = bb2 < bd2 ? bb2 : bd2; 505 if (i > bs2) 506 i = bs2; 507 if (i > 0) { 508 bb2 -= i; 509 bd2 -= i; 510 bs2 -= i; 511 } 512 if (bb5 > 0) { 513 bs = pow5mult(bs, bb5); 514 bb1 = mult(bs, bb); 515 Bfree(bb); 516 bb = bb1; 517 } 518 if (bb2 > 0) 519 bb = lshift(bb, bb2); 520 if (bd5 > 0) 521 bd = pow5mult(bd, bd5); 522 if (bd2 > 0) 523 bd = lshift(bd, bd2); 524 if (bs2 > 0) 525 bs = lshift(bs, bs2); 526 delta = diff(bb, bd); 527 dsign = delta->sign; 528 delta->sign = 0; 529 i = cmp(delta, bs); 530 if (i < 0) { 531 /* Error is less than half an ulp -- check for 532 * special case of mantissa a power of two. 533 */ 534 if (dsign || word1(rv) || word0(rv) & Bndry_mask) 535 break; 536 delta = lshift(delta,Log2P); 537 if (cmp(delta, bs) > 0) 538 goto drop_down; 539 break; 540 } 541 if (i == 0) { 542 /* exactly half-way between */ 543 if (dsign) { 544 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 545 && word1(rv) == 0xffffffff) { 546 /*boundary case -- increment exponent*/ 547 word0(rv) = (word0(rv) & Exp_mask) 548 + Exp_msk1 549 #ifdef IBM 550 | Exp_msk1 >> 4 551 #endif 552 ; 553 word1(rv) = 0; 554 break; 555 } 556 } 557 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 558 drop_down: 559 /* boundary case -- decrement exponent */ 560 #ifdef Sudden_Underflow 561 L = word0(rv) & Exp_mask; 562 #ifdef IBM 563 if (L < Exp_msk1) 564 #else 565 if (L <= Exp_msk1) 566 #endif 567 goto undfl; 568 L -= Exp_msk1; 569 #else 570 L = (word0(rv) & Exp_mask) - Exp_msk1; 571 #endif 572 word0(rv) = L | Bndry_mask1; 573 word1(rv) = 0xffffffff; 574 #ifdef IBM 575 goto cont; 576 #else 577 break; 578 #endif 579 } 580 #ifndef ROUND_BIASED 581 if (!(word1(rv) & LSB)) 582 break; 583 #endif 584 if (dsign) 585 rv.d += ulp(rv.d); 586 #ifndef ROUND_BIASED 587 else { 588 rv.d -= ulp(rv.d); 589 #ifndef Sudden_Underflow 590 if (rv.d == 0) 591 goto undfl; 592 #endif 593 } 594 #endif 595 break; 596 } 597 if ((aadj = ratio(delta, bs)) <= 2.) { 598 if (dsign) 599 aadj = aadj1 = 1.; 600 else if (word1(rv) || word0(rv) & Bndry_mask) { 601 #ifndef Sudden_Underflow 602 if (word1(rv) == Tiny1 && !word0(rv)) 603 goto undfl; 604 #endif 605 aadj = 1.; 606 aadj1 = -1.; 607 } 608 else { 609 /* special case -- power of FLT_RADIX to be */ 610 /* rounded down... */ 611 612 if (aadj < 2./FLT_RADIX) 613 aadj = 1./FLT_RADIX; 614 else 615 aadj *= 0.5; 616 aadj1 = -aadj; 617 } 618 } 619 else { 620 aadj *= 0.5; 621 aadj1 = dsign ? aadj : -aadj; 622 #ifdef Check_FLT_ROUNDS 623 switch(FLT_ROUNDS) { 624 case 2: /* towards +infinity */ 625 aadj1 -= 0.5; 626 break; 627 case 0: /* towards 0 */ 628 case 3: /* towards -infinity */ 629 aadj1 += 0.5; 630 } 631 #else 632 if (FLT_ROUNDS == 0) 633 aadj1 += 0.5; 634 #endif 635 } 636 y = word0(rv) & Exp_mask; 637 638 /* Check for overflow */ 639 640 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 641 rv0.d = rv.d; 642 word0(rv) -= P*Exp_msk1; 643 adj = aadj1 * ulp(rv.d); 644 rv.d += adj; 645 if ((word0(rv) & Exp_mask) >= 646 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 647 if (word0(rv0) == Big0 && word1(rv0) == Big1) 648 goto ovfl; 649 word0(rv) = Big0; 650 word1(rv) = Big1; 651 goto cont; 652 } 653 else 654 word0(rv) += P*Exp_msk1; 655 } 656 else { 657 #ifdef Sudden_Underflow 658 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 659 rv0.d = rv.d; 660 word0(rv) += P*Exp_msk1; 661 adj = aadj1 * ulp(rv.d); 662 rv.d += adj; 663 #ifdef IBM 664 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 665 #else 666 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 667 #endif 668 { 669 if (word0(rv0) == Tiny0 670 && word1(rv0) == Tiny1) 671 goto undfl; 672 word0(rv) = Tiny0; 673 word1(rv) = Tiny1; 674 goto cont; 675 } 676 else 677 word0(rv) -= P*Exp_msk1; 678 } 679 else { 680 adj = aadj1 * ulp(rv.d); 681 rv.d += adj; 682 } 683 #else 684 /* Compute adj so that the IEEE rounding rules will 685 * correctly round rv + adj in some half-way cases. 686 * If rv * ulp(rv) is denormalized (i.e., 687 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 688 * trouble from bits lost to denormalization; 689 * example: 1.2e-307 . 690 */ 691 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { 692 aadj1 = (double)(int)(aadj + 0.5); 693 if (!dsign) 694 aadj1 = -aadj1; 695 } 696 adj = aadj1 * ulp(rv.d); 697 rv.d += adj; 698 #endif 699 } 700 z = word0(rv) & Exp_mask; 701 if (y == z) { 702 /* Can we stop now? */ 703 L = aadj; 704 aadj -= L; 705 /* The tolerances below are conservative. */ 706 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 707 if (aadj < .4999999 || aadj > .5000001) 708 break; 709 } 710 else if (aadj < .4999999/FLT_RADIX) 711 break; 712 } 713 cont: 714 Bfree(bb); 715 Bfree(bd); 716 Bfree(bs); 717 Bfree(delta); 718 } 719 retfree: 720 Bfree(bb); 721 Bfree(bd); 722 Bfree(bs); 723 Bfree(bd0); 724 Bfree(delta); 725 ret: 726 if (se) 727 *se = (char *)s; 728 return sign ? -rv.d : rv.d; 729 } 730