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