1 /**************************************************************** 2 3 The author of this software is David M. Gay. 4 5 Copyright (C) 1998-2001 by Lucent Technologies 6 All Rights Reserved 7 8 Permission to use, copy, modify, and distribute this software and 9 its documentation for any purpose and without fee is hereby 10 granted, provided that the above copyright notice appear in all 11 copies and that both that the copyright notice and this 12 permission notice and warranty disclaimer appear in supporting 13 documentation, and that the name of Lucent or any of its entities 14 not be used in advertising or publicity pertaining to 15 distribution of the software without specific, written prior 16 permission. 17 18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25 THIS SOFTWARE. 26 27 ****************************************************************/ 28 29 /* Please send bug reports to David M. Gay (dmg at acm dot org, 30 * with " at " changed at "@" and " dot " changed to "."). */ 31 32 #include "gdtoaimp.h" 33 #ifndef NO_FENV_H 34 #include <fenv.h> 35 #endif 36 37 #ifdef USE_LOCALE 38 #include "locale.h" 39 #endif 40 41 #ifdef IEEE_Arith 42 #ifndef NO_IEEE_Scale 43 #define Avoid_Underflow 44 #undef tinytens 45 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */ 46 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 47 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 48 9007199254740992.*9007199254740992.e-256 49 }; 50 #endif 51 #endif 52 53 #ifdef Honor_FLT_ROUNDS 54 #undef Check_FLT_ROUNDS 55 #define Check_FLT_ROUNDS 56 #else 57 #define Rounding Flt_Rounds 58 #endif 59 60 #ifdef Avoid_Underflow /*{*/ 61 static double 62 sulp 63 #ifdef KR_headers 64 (x, scale) U *x; int scale; 65 #else 66 (U *x, int scale) 67 #endif 68 { 69 U u; 70 double rv; 71 int i; 72 73 rv = ulp(x); 74 if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0) 75 return rv; /* Is there an example where i <= 0 ? */ 76 word0(&u) = Exp_1 + (i << Exp_shift); 77 word1(&u) = 0; 78 return rv * u.d; 79 } 80 #endif /*}*/ 81 82 double 83 strtod 84 #ifdef KR_headers 85 (s00, se) CONST char *s00; char **se; 86 #else 87 (CONST char *s00, char **se) 88 #endif 89 { 90 #ifdef Avoid_Underflow 91 int scale; 92 #endif 93 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, 94 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 95 CONST char *s, *s0, *s1; 96 double aadj; 97 Long L; 98 U adj, aadj1, rv, rv0; 99 ULong y, z; 100 Bigint *bb = NULL, *bb1, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL; 101 #ifdef Avoid_Underflow 102 ULong Lsb, Lsb1; 103 #endif 104 #ifdef SET_INEXACT 105 int inexact, oldinexact; 106 #endif 107 #ifdef USE_LOCALE /*{{*/ 108 #ifdef NO_LOCALE_CACHE 109 char *decimalpoint = localeconv()->decimal_point; 110 int dplen = strlen(decimalpoint); 111 #else 112 char *decimalpoint; 113 static char *decimalpoint_cache; 114 static int dplen; 115 if (!(s0 = decimalpoint_cache)) { 116 s0 = localeconv()->decimal_point; 117 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { 118 strlcpy(decimalpoint_cache, s0, strlen(s0) + 1); 119 s0 = decimalpoint_cache; 120 } 121 dplen = strlen(s0); 122 } 123 decimalpoint = (char*)s0; 124 #endif /*NO_LOCALE_CACHE*/ 125 #else /*USE_LOCALE}{*/ 126 #define dplen 1 127 #endif /*USE_LOCALE}}*/ 128 129 #ifdef Honor_FLT_ROUNDS /*{*/ 130 int Rounding; 131 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ 132 Rounding = Flt_Rounds; 133 #else /*}{*/ 134 Rounding = 1; 135 switch(fegetround()) { 136 case FE_TOWARDZERO: Rounding = 0; break; 137 case FE_UPWARD: Rounding = 2; break; 138 case FE_DOWNWARD: Rounding = 3; 139 } 140 #endif /*}}*/ 141 #endif /*}*/ 142 143 sign = nz0 = nz = decpt = 0; 144 dval(&rv) = 0.; 145 for(s = s00;;s++) switch(*s) { 146 case '-': 147 sign = 1; 148 /* no break */ 149 case '+': 150 if (*++s) 151 goto break2; 152 /* no break */ 153 case 0: 154 goto ret0; 155 case '\t': 156 case '\n': 157 case '\v': 158 case '\f': 159 case '\r': 160 case ' ': 161 continue; 162 default: 163 goto break2; 164 } 165 break2: 166 if (*s == '0') { 167 #ifndef NO_HEX_FP /*{*/ 168 { 169 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 170 Long exp; 171 ULong bits[2]; 172 switch(s[1]) { 173 case 'x': 174 case 'X': 175 { 176 #ifdef Honor_FLT_ROUNDS 177 FPI fpi1 = fpi; 178 fpi1.rounding = Rounding; 179 #else 180 #define fpi1 fpi 181 #endif 182 switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { 183 case STRTOG_NoMemory: 184 goto ovfl; 185 case STRTOG_NoNumber: 186 s = s00; 187 sign = 0; 188 case STRTOG_Zero: 189 break; 190 default: 191 if (bb) { 192 copybits(bits, fpi.nbits, bb); 193 Bfree(bb); 194 } 195 ULtod(((U*)&rv)->L, bits, exp, i); 196 }} 197 goto ret; 198 } 199 } 200 #endif /*}*/ 201 nz0 = 1; 202 while(*++s == '0') ; 203 if (!*s) 204 goto ret; 205 } 206 s0 = s; 207 y = z = 0; 208 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 209 if (nd < 9) 210 y = 10*y + c - '0'; 211 else if (nd < 16) 212 z = 10*z + c - '0'; 213 nd0 = nd; 214 #ifdef USE_LOCALE 215 if (c == *decimalpoint) { 216 for(i = 1; decimalpoint[i]; ++i) 217 if (s[i] != decimalpoint[i]) 218 goto dig_done; 219 s += i; 220 c = *s; 221 #else 222 if (c == '.') { 223 c = *++s; 224 #endif 225 decpt = 1; 226 if (!nd) { 227 for(; c == '0'; c = *++s) 228 nz++; 229 if (c > '0' && c <= '9') { 230 s0 = s; 231 nf += nz; 232 nz = 0; 233 goto have_dig; 234 } 235 goto dig_done; 236 } 237 for(; c >= '0' && c <= '9'; c = *++s) { 238 have_dig: 239 nz++; 240 if (c -= '0') { 241 nf += nz; 242 for(i = 1; i < nz; i++) 243 if (nd++ < 9) 244 y *= 10; 245 else if (nd <= DBL_DIG + 1) 246 z *= 10; 247 if (nd++ < 9) 248 y = 10*y + c; 249 else if (nd <= DBL_DIG + 1) 250 z = 10*z + c; 251 nz = 0; 252 } 253 } 254 }/*}*/ 255 dig_done: 256 e = 0; 257 if (c == 'e' || c == 'E') { 258 if (!nd && !nz && !nz0) { 259 goto ret0; 260 } 261 s00 = s; 262 esign = 0; 263 switch(c = *++s) { 264 case '-': 265 esign = 1; 266 case '+': 267 c = *++s; 268 } 269 if (c >= '0' && c <= '9') { 270 while(c == '0') 271 c = *++s; 272 if (c > '0' && c <= '9') { 273 L = c - '0'; 274 s1 = s; 275 while((c = *++s) >= '0' && c <= '9') 276 L = 10*L + c - '0'; 277 if (s - s1 > 8 || L > 19999) 278 /* Avoid confusion from exponents 279 * so large that e might overflow. 280 */ 281 e = 19999; /* safe for 16 bit ints */ 282 else 283 e = (int)L; 284 if (esign) 285 e = -e; 286 } 287 else 288 e = 0; 289 } 290 else 291 s = s00; 292 } 293 if (!nd) { 294 if (!nz && !nz0) { 295 #ifdef INFNAN_CHECK 296 /* Check for Nan and Infinity */ 297 ULong bits[2]; 298 static FPI fpinan = /* only 52 explicit bits */ 299 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 300 if (!decpt) 301 switch(c) { 302 case 'i': 303 case 'I': 304 if (match(&s,"nf")) { 305 --s; 306 if (!match(&s,"inity")) 307 ++s; 308 word0(&rv) = 0x7ff00000; 309 word1(&rv) = 0; 310 goto ret; 311 } 312 break; 313 case 'n': 314 case 'N': 315 if (match(&s, "an")) { 316 #ifndef No_Hex_NaN 317 if (*s == '(' /*)*/ 318 && hexnan(&s, &fpinan, bits) 319 == STRTOG_NaNbits) { 320 word0(&rv) = 0x7ff00000 | bits[1]; 321 word1(&rv) = bits[0]; 322 } 323 else { 324 #endif 325 word0(&rv) = NAN_WORD0; 326 word1(&rv) = NAN_WORD1; 327 #ifndef No_Hex_NaN 328 } 329 #endif 330 goto ret; 331 } 332 } 333 #endif /* INFNAN_CHECK */ 334 ret0: 335 s = s00; 336 sign = 0; 337 } 338 goto ret; 339 } 340 e1 = e -= nf; 341 342 /* Now we have nd0 digits, starting at s0, followed by a 343 * decimal point, followed by nd-nd0 digits. The number we're 344 * after is the integer represented by those digits times 345 * 10**e */ 346 347 if (!nd0) 348 nd0 = nd; 349 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 350 dval(&rv) = y; 351 if (k > 9) { 352 #ifdef SET_INEXACT 353 if (k > DBL_DIG) 354 oldinexact = get_inexact(); 355 #endif 356 dval(&rv) = tens[k - 9] * dval(&rv) + z; 357 } 358 if (nd <= DBL_DIG 359 #ifndef RND_PRODQUOT 360 #ifndef Honor_FLT_ROUNDS 361 && Flt_Rounds == 1 362 #endif 363 #endif 364 ) { 365 if (!e) 366 goto ret; 367 #ifndef ROUND_BIASED_without_Round_Up 368 if (e > 0) { 369 if (e <= Ten_pmax) { 370 #ifdef VAX 371 goto vax_ovfl_check; 372 #else 373 #ifdef Honor_FLT_ROUNDS 374 /* round correctly FLT_ROUNDS = 2 or 3 */ 375 if (sign) { 376 rv.d = -rv.d; 377 sign = 0; 378 } 379 #endif 380 /* rv = */ rounded_product(dval(&rv), tens[e]); 381 goto ret; 382 #endif 383 } 384 i = DBL_DIG - nd; 385 if (e <= Ten_pmax + i) { 386 /* A fancier test would sometimes let us do 387 * this for larger i values. 388 */ 389 #ifdef Honor_FLT_ROUNDS 390 /* round correctly FLT_ROUNDS = 2 or 3 */ 391 if (sign) { 392 rv.d = -rv.d; 393 sign = 0; 394 } 395 #endif 396 e -= i; 397 dval(&rv) *= tens[i]; 398 #ifdef VAX 399 /* VAX exponent range is so narrow we must 400 * worry about overflow here... 401 */ 402 vax_ovfl_check: 403 word0(&rv) -= P*Exp_msk1; 404 /* rv = */ rounded_product(dval(&rv), tens[e]); 405 if ((word0(&rv) & Exp_mask) 406 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 407 goto ovfl; 408 word0(&rv) += P*Exp_msk1; 409 #else 410 /* rv = */ rounded_product(dval(&rv), tens[e]); 411 #endif 412 goto ret; 413 } 414 } 415 #ifndef Inaccurate_Divide 416 else if (e >= -Ten_pmax) { 417 #ifdef Honor_FLT_ROUNDS 418 /* round correctly FLT_ROUNDS = 2 or 3 */ 419 if (sign) { 420 rv.d = -rv.d; 421 sign = 0; 422 } 423 #endif 424 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 425 goto ret; 426 } 427 #endif 428 #endif /* ROUND_BIASED_without_Round_Up */ 429 } 430 e1 += nd - k; 431 432 #ifdef IEEE_Arith 433 #ifdef SET_INEXACT 434 inexact = 1; 435 if (k <= DBL_DIG) 436 oldinexact = get_inexact(); 437 #endif 438 #ifdef Avoid_Underflow 439 scale = 0; 440 #endif 441 #ifdef Honor_FLT_ROUNDS 442 if (Rounding >= 2) { 443 if (sign) 444 Rounding = Rounding == 2 ? 0 : 2; 445 else 446 if (Rounding != 2) 447 Rounding = 0; 448 } 449 #endif 450 #endif /*IEEE_Arith*/ 451 452 /* Get starting approximation = rv * 10**e1 */ 453 454 if (e1 > 0) { 455 if ( (i = e1 & 15) !=0) 456 dval(&rv) *= tens[i]; 457 if (e1 &= ~15) { 458 if (e1 > DBL_MAX_10_EXP) { 459 ovfl: 460 /* Can't trust HUGE_VAL */ 461 #ifdef IEEE_Arith 462 #ifdef Honor_FLT_ROUNDS 463 switch(Rounding) { 464 case 0: /* toward 0 */ 465 case 3: /* toward -infinity */ 466 word0(&rv) = Big0; 467 word1(&rv) = Big1; 468 break; 469 default: 470 word0(&rv) = Exp_mask; 471 word1(&rv) = 0; 472 } 473 #else /*Honor_FLT_ROUNDS*/ 474 word0(&rv) = Exp_mask; 475 word1(&rv) = 0; 476 #endif /*Honor_FLT_ROUNDS*/ 477 #ifdef SET_INEXACT 478 /* set overflow bit */ 479 dval(&rv0) = 1e300; 480 dval(&rv0) *= dval(&rv0); 481 #endif 482 #else /*IEEE_Arith*/ 483 word0(&rv) = Big0; 484 word1(&rv) = Big1; 485 #endif /*IEEE_Arith*/ 486 range_err: 487 if (bd0) { 488 Bfree(bb); 489 Bfree(bd); 490 Bfree(bs); 491 Bfree(bd0); 492 Bfree(delta); 493 } 494 #ifndef NO_ERRNO 495 errno = ERANGE; 496 #endif 497 goto ret; 498 } 499 e1 >>= 4; 500 for(j = 0; e1 > 1; j++, e1 >>= 1) 501 if (e1 & 1) 502 dval(&rv) *= bigtens[j]; 503 /* The last multiplication could overflow. */ 504 word0(&rv) -= P*Exp_msk1; 505 dval(&rv) *= bigtens[j]; 506 if ((z = word0(&rv) & Exp_mask) 507 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 508 goto ovfl; 509 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 510 /* set to largest number */ 511 /* (Can't trust DBL_MAX) */ 512 word0(&rv) = Big0; 513 word1(&rv) = Big1; 514 } 515 else 516 word0(&rv) += P*Exp_msk1; 517 } 518 } 519 else if (e1 < 0) { 520 e1 = -e1; 521 if ( (i = e1 & 15) !=0) 522 dval(&rv) /= tens[i]; 523 if (e1 >>= 4) { 524 if (e1 >= 1 << n_bigtens) 525 goto undfl; 526 #ifdef Avoid_Underflow 527 if (e1 & Scale_Bit) 528 scale = 2*P; 529 for(j = 0; e1 > 0; j++, e1 >>= 1) 530 if (e1 & 1) 531 dval(&rv) *= tinytens[j]; 532 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 533 >> Exp_shift)) > 0) { 534 /* scaled rv is denormal; zap j low bits */ 535 if (j >= 32) { 536 word1(&rv) = 0; 537 if (j >= 53) 538 word0(&rv) = (P+2)*Exp_msk1; 539 else 540 word0(&rv) &= 0xffffffff << (j-32); 541 } 542 else 543 word1(&rv) &= 0xffffffff << j; 544 } 545 #else 546 for(j = 0; e1 > 1; j++, e1 >>= 1) 547 if (e1 & 1) 548 dval(&rv) *= tinytens[j]; 549 /* The last multiplication could underflow. */ 550 dval(&rv0) = dval(&rv); 551 dval(&rv) *= tinytens[j]; 552 if (!dval(&rv)) { 553 dval(&rv) = 2.*dval(&rv0); 554 dval(&rv) *= tinytens[j]; 555 #endif 556 if (!dval(&rv)) { 557 undfl: 558 dval(&rv) = 0.; 559 goto range_err; 560 } 561 #ifndef Avoid_Underflow 562 word0(&rv) = Tiny0; 563 word1(&rv) = Tiny1; 564 /* The refinement below will clean 565 * this approximation up. 566 */ 567 } 568 #endif 569 } 570 } 571 572 /* Now the hard part -- adjusting rv to the correct value.*/ 573 574 /* Put digits into bd: true value = bd * 10^e */ 575 576 bd0 = s2b(s0, nd0, nd, y, dplen); 577 if (bd0 == NULL) 578 goto ovfl; 579 580 for(;;) { 581 bd = Balloc(bd0->k); 582 if (bd == NULL) 583 goto ovfl; 584 Bcopy(bd, bd0); 585 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 586 if (bb == NULL) 587 goto ovfl; 588 bs = i2b(1); 589 if (bs == NULL) 590 goto ovfl; 591 592 if (e >= 0) { 593 bb2 = bb5 = 0; 594 bd2 = bd5 = e; 595 } 596 else { 597 bb2 = bb5 = -e; 598 bd2 = bd5 = 0; 599 } 600 if (bbe >= 0) 601 bb2 += bbe; 602 else 603 bd2 -= bbe; 604 bs2 = bb2; 605 #ifdef Honor_FLT_ROUNDS 606 if (Rounding != 1) 607 bs2++; 608 #endif 609 #ifdef Avoid_Underflow 610 Lsb = LSB; 611 Lsb1 = 0; 612 j = bbe - scale; 613 i = j + bbbits - 1; /* logb(rv) */ 614 j = P + 1 - bbbits; 615 if (i < Emin) { /* denormal */ 616 i = Emin - i; 617 j -= i; 618 if (i < 32) 619 Lsb <<= i; 620 else 621 Lsb1 = Lsb << (i-32); 622 } 623 #else /*Avoid_Underflow*/ 624 #ifdef Sudden_Underflow 625 #ifdef IBM 626 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 627 #else 628 j = P + 1 - bbbits; 629 #endif 630 #else /*Sudden_Underflow*/ 631 j = bbe; 632 i = j + bbbits - 1; /* logb(&rv) */ 633 if (i < Emin) /* denormal */ 634 j += P - Emin; 635 else 636 j = P + 1 - bbbits; 637 #endif /*Sudden_Underflow*/ 638 #endif /*Avoid_Underflow*/ 639 bb2 += j; 640 bd2 += j; 641 #ifdef Avoid_Underflow 642 bd2 += scale; 643 #endif 644 i = bb2 < bd2 ? bb2 : bd2; 645 if (i > bs2) 646 i = bs2; 647 if (i > 0) { 648 bb2 -= i; 649 bd2 -= i; 650 bs2 -= i; 651 } 652 if (bb5 > 0) { 653 bs = pow5mult(bs, bb5); 654 if (bs == NULL) 655 goto ovfl; 656 bb1 = mult(bs, bb); 657 if (bb1 == NULL) 658 goto ovfl; 659 Bfree(bb); 660 bb = bb1; 661 } 662 if (bb2 > 0) { 663 bb = lshift(bb, bb2); 664 if (bb == NULL) 665 goto ovfl; 666 } 667 if (bd5 > 0) { 668 bd = pow5mult(bd, bd5); 669 if (bd == NULL) 670 goto ovfl; 671 } 672 if (bd2 > 0) { 673 bd = lshift(bd, bd2); 674 if (bd == NULL) 675 goto ovfl; 676 } 677 if (bs2 > 0) { 678 bs = lshift(bs, bs2); 679 if (bs == NULL) 680 goto ovfl; 681 } 682 delta = diff(bb, bd); 683 if (delta == NULL) 684 goto ovfl; 685 dsign = delta->sign; 686 delta->sign = 0; 687 i = cmp(delta, bs); 688 #ifdef Honor_FLT_ROUNDS 689 if (Rounding != 1) { 690 if (i < 0) { 691 /* Error is less than an ulp */ 692 if (!delta->x[0] && delta->wds <= 1) { 693 /* exact */ 694 #ifdef SET_INEXACT 695 inexact = 0; 696 #endif 697 break; 698 } 699 if (Rounding) { 700 if (dsign) { 701 dval(&adj) = 1.; 702 goto apply_adj; 703 } 704 } 705 else if (!dsign) { 706 dval(&adj) = -1.; 707 if (!word1(&rv) 708 && !(word0(&rv) & Frac_mask)) { 709 y = word0(&rv) & Exp_mask; 710 #ifdef Avoid_Underflow 711 if (!scale || y > 2*P*Exp_msk1) 712 #else 713 if (y) 714 #endif 715 { 716 delta = lshift(delta,Log2P); 717 if (delta == NULL) 718 goto ovfl; 719 if (cmp(delta, bs) <= 0) 720 dval(&adj) = -0.5; 721 } 722 } 723 apply_adj: 724 #ifdef Avoid_Underflow 725 if (scale && (y = word0(&rv) & Exp_mask) 726 <= 2*P*Exp_msk1) 727 word0(&adj) += (2*P+1)*Exp_msk1 - y; 728 #else 729 #ifdef Sudden_Underflow 730 if ((word0(&rv) & Exp_mask) <= 731 P*Exp_msk1) { 732 word0(&rv) += P*Exp_msk1; 733 dval(&rv) += adj*ulp(&rv); 734 word0(&rv) -= P*Exp_msk1; 735 } 736 else 737 #endif /*Sudden_Underflow*/ 738 #endif /*Avoid_Underflow*/ 739 dval(&rv) += adj.d*ulp(&rv); 740 } 741 break; 742 } 743 dval(&adj) = ratio(delta, bs); 744 if (adj.d < 1.) 745 dval(&adj) = 1.; 746 if (adj.d <= 0x7ffffffe) { 747 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */ 748 y = adj.d; 749 if (y != adj.d) { 750 if (!((Rounding>>1) ^ dsign)) 751 y++; 752 dval(&adj) = y; 753 } 754 } 755 #ifdef Avoid_Underflow 756 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 757 word0(&adj) += (2*P+1)*Exp_msk1 - y; 758 #else 759 #ifdef Sudden_Underflow 760 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 761 word0(&rv) += P*Exp_msk1; 762 dval(&adj) *= ulp(&rv); 763 if (dsign) 764 dval(&rv) += adj; 765 else 766 dval(&rv) -= adj; 767 word0(&rv) -= P*Exp_msk1; 768 goto cont; 769 } 770 #endif /*Sudden_Underflow*/ 771 #endif /*Avoid_Underflow*/ 772 dval(&adj) *= ulp(&rv); 773 if (dsign) { 774 if (word0(&rv) == Big0 && word1(&rv) == Big1) 775 goto ovfl; 776 dval(&rv) += adj.d; 777 } 778 else 779 dval(&rv) -= adj.d; 780 goto cont; 781 } 782 #endif /*Honor_FLT_ROUNDS*/ 783 784 if (i < 0) { 785 /* Error is less than half an ulp -- check for 786 * special case of mantissa a power of two. 787 */ 788 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask 789 #ifdef IEEE_Arith 790 #ifdef Avoid_Underflow 791 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 792 #else 793 || (word0(&rv) & Exp_mask) <= Exp_msk1 794 #endif 795 #endif 796 ) { 797 #ifdef SET_INEXACT 798 if (!delta->x[0] && delta->wds <= 1) 799 inexact = 0; 800 #endif 801 break; 802 } 803 if (!delta->x[0] && delta->wds <= 1) { 804 /* exact result */ 805 #ifdef SET_INEXACT 806 inexact = 0; 807 #endif 808 break; 809 } 810 delta = lshift(delta,Log2P); 811 if (delta == NULL) 812 goto ovfl; 813 if (cmp(delta, bs) > 0) 814 goto drop_down; 815 break; 816 } 817 if (i == 0) { 818 /* exactly half-way between */ 819 if (dsign) { 820 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 821 && word1(&rv) == ( 822 #ifdef Avoid_Underflow 823 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 824 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 825 #endif 826 0xffffffff)) { 827 /*boundary case -- increment exponent*/ 828 if (word0(&rv) == Big0 && word1(&rv) == Big1) 829 goto ovfl; 830 word0(&rv) = (word0(&rv) & Exp_mask) 831 + Exp_msk1 832 #ifdef IBM 833 | Exp_msk1 >> 4 834 #endif 835 ; 836 word1(&rv) = 0; 837 #ifdef Avoid_Underflow 838 dsign = 0; 839 #endif 840 break; 841 } 842 } 843 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 844 drop_down: 845 /* boundary case -- decrement exponent */ 846 #ifdef Sudden_Underflow /*{{*/ 847 L = word0(&rv) & Exp_mask; 848 #ifdef IBM 849 if (L < Exp_msk1) 850 #else 851 #ifdef Avoid_Underflow 852 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 853 #else 854 if (L <= Exp_msk1) 855 #endif /*Avoid_Underflow*/ 856 #endif /*IBM*/ 857 goto undfl; 858 L -= Exp_msk1; 859 #else /*Sudden_Underflow}{*/ 860 #ifdef Avoid_Underflow 861 if (scale) { 862 L = word0(&rv) & Exp_mask; 863 if (L <= (2*P+1)*Exp_msk1) { 864 if (L > (P+2)*Exp_msk1) 865 /* round even ==> */ 866 /* accept rv */ 867 break; 868 /* rv = smallest denormal */ 869 goto undfl; 870 } 871 } 872 #endif /*Avoid_Underflow*/ 873 L = (word0(&rv) & Exp_mask) - Exp_msk1; 874 #endif /*Sudden_Underflow}}*/ 875 word0(&rv) = L | Bndry_mask1; 876 word1(&rv) = 0xffffffff; 877 #ifdef IBM 878 goto cont; 879 #else 880 break; 881 #endif 882 } 883 #ifndef ROUND_BIASED 884 #ifdef Avoid_Underflow 885 if (Lsb1) { 886 if (!(word0(&rv) & Lsb1)) 887 break; 888 } 889 else if (!(word1(&rv) & Lsb)) 890 break; 891 #else 892 if (!(word1(&rv) & LSB)) 893 break; 894 #endif 895 #endif 896 if (dsign) 897 #ifdef Avoid_Underflow 898 dval(&rv) += sulp(&rv, scale); 899 #else 900 dval(&rv) += ulp(&rv); 901 #endif 902 #ifndef ROUND_BIASED 903 else { 904 #ifdef Avoid_Underflow 905 dval(&rv) -= sulp(&rv, scale); 906 #else 907 dval(&rv) -= ulp(&rv); 908 #endif 909 #ifndef Sudden_Underflow 910 if (!dval(&rv)) 911 goto undfl; 912 #endif 913 } 914 #ifdef Avoid_Underflow 915 dsign = 1 - dsign; 916 #endif 917 #endif 918 break; 919 } 920 if ((aadj = ratio(delta, bs)) <= 2.) { 921 if (dsign) 922 aadj = dval(&aadj1) = 1.; 923 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 924 #ifndef Sudden_Underflow 925 if (word1(&rv) == Tiny1 && !word0(&rv)) 926 goto undfl; 927 #endif 928 aadj = 1.; 929 dval(&aadj1) = -1.; 930 } 931 else { 932 /* special case -- power of FLT_RADIX to be */ 933 /* rounded down... */ 934 935 if (aadj < 2./FLT_RADIX) 936 aadj = 1./FLT_RADIX; 937 else 938 aadj *= 0.5; 939 dval(&aadj1) = -aadj; 940 } 941 } 942 else { 943 aadj *= 0.5; 944 dval(&aadj1) = dsign ? aadj : -aadj; 945 #ifdef Check_FLT_ROUNDS 946 switch(Rounding) { 947 case 2: /* towards +infinity */ 948 dval(&aadj1) -= 0.5; 949 break; 950 case 0: /* towards 0 */ 951 case 3: /* towards -infinity */ 952 dval(&aadj1) += 0.5; 953 } 954 #else 955 if (Flt_Rounds == 0) 956 dval(&aadj1) += 0.5; 957 #endif /*Check_FLT_ROUNDS*/ 958 } 959 y = word0(&rv) & Exp_mask; 960 961 /* Check for overflow */ 962 963 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 964 dval(&rv0) = dval(&rv); 965 word0(&rv) -= P*Exp_msk1; 966 dval(&adj) = dval(&aadj1) * ulp(&rv); 967 dval(&rv) += dval(&adj); 968 if ((word0(&rv) & Exp_mask) >= 969 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 970 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 971 goto ovfl; 972 word0(&rv) = Big0; 973 word1(&rv) = Big1; 974 goto cont; 975 } 976 else 977 word0(&rv) += P*Exp_msk1; 978 } 979 else { 980 #ifdef Avoid_Underflow 981 if (scale && y <= 2*P*Exp_msk1) { 982 if (aadj <= 0x7fffffff) { 983 if ((z = aadj) <= 0) 984 z = 1; 985 aadj = z; 986 dval(&aadj1) = dsign ? aadj : -aadj; 987 } 988 word0(&aadj1) += (2*P+1)*Exp_msk1 - y; 989 } 990 dval(&adj) = dval(&aadj1) * ulp(&rv); 991 dval(&rv) += dval(&adj); 992 #else 993 #ifdef Sudden_Underflow 994 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 995 dval(&rv0) = dval(&rv); 996 word0(&rv) += P*Exp_msk1; 997 dval(&adj) = dval(&aadj1) * ulp(&rv); 998 dval(&rv) += dval(&adj); 999 #ifdef IBM 1000 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) 1001 #else 1002 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) 1003 #endif 1004 { 1005 if (word0(&rv0) == Tiny0 1006 && word1(&rv0) == Tiny1) 1007 goto undfl; 1008 word0(&rv) = Tiny0; 1009 word1(&rv) = Tiny1; 1010 goto cont; 1011 } 1012 else 1013 word0(&rv) -= P*Exp_msk1; 1014 } 1015 else { 1016 dval(&adj) = dval(&aadj1) * ulp(&rv); 1017 dval(&rv) += dval(&adj); 1018 } 1019 #else /*Sudden_Underflow*/ 1020 /* Compute dval(&adj) so that the IEEE rounding rules will 1021 * correctly round rv + dval(&adj) in some half-way cases. 1022 * If rv * ulp(&rv) is denormalized (i.e., 1023 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 1024 * trouble from bits lost to denormalization; 1025 * example: 1.2e-307 . 1026 */ 1027 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 1028 dval(&aadj1) = (double)(int)(aadj + 0.5); 1029 if (!dsign) 1030 dval(&aadj1) = -dval(&aadj1); 1031 } 1032 dval(&adj) = dval(&aadj1) * ulp(&rv); 1033 dval(&rv) += adj; 1034 #endif /*Sudden_Underflow*/ 1035 #endif /*Avoid_Underflow*/ 1036 } 1037 z = word0(&rv) & Exp_mask; 1038 #ifndef SET_INEXACT 1039 #ifdef Avoid_Underflow 1040 if (!scale) 1041 #endif 1042 if (y == z) { 1043 /* Can we stop now? */ 1044 L = (Long)aadj; 1045 aadj -= L; 1046 /* The tolerances below are conservative. */ 1047 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 1048 if (aadj < .4999999 || aadj > .5000001) 1049 break; 1050 } 1051 else if (aadj < .4999999/FLT_RADIX) 1052 break; 1053 } 1054 #endif 1055 cont: 1056 Bfree(bb); 1057 Bfree(bd); 1058 Bfree(bs); 1059 Bfree(delta); 1060 } 1061 Bfree(bb); 1062 Bfree(bd); 1063 Bfree(bs); 1064 Bfree(bd0); 1065 Bfree(delta); 1066 #ifdef SET_INEXACT 1067 if (inexact) { 1068 if (!oldinexact) { 1069 word0(&rv0) = Exp_1 + (70 << Exp_shift); 1070 word1(&rv0) = 0; 1071 dval(&rv0) += 1.; 1072 } 1073 } 1074 else if (!oldinexact) 1075 clear_inexact(); 1076 #endif 1077 #ifdef Avoid_Underflow 1078 if (scale) { 1079 word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 1080 word1(&rv0) = 0; 1081 dval(&rv) *= dval(&rv0); 1082 #ifndef NO_ERRNO 1083 /* try to avoid the bug of testing an 8087 register value */ 1084 #ifdef IEEE_Arith 1085 if (!(word0(&rv) & Exp_mask)) 1086 #else 1087 if (word0(&rv) == 0 && word1(&rv) == 0) 1088 #endif 1089 errno = ERANGE; 1090 #endif 1091 } 1092 #endif /* Avoid_Underflow */ 1093 #ifdef SET_INEXACT 1094 if (inexact && !(word0(&rv) & Exp_mask)) { 1095 /* set underflow bit */ 1096 dval(&rv0) = 1e-300; 1097 dval(&rv0) *= dval(&rv0); 1098 } 1099 #endif 1100 ret: 1101 if (se) 1102 *se = (char *)s; 1103 return sign ? -dval(&rv) : dval(&rv); 1104 } 1105 DEF_STRONG(strtod); 1106
