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