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