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