1 /* $NetBSD: strtod.c,v 1.17 2020/09/18 14:19:34 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 (cmp(delta, bs) <= 0) 716 dval(&adj) = -0.5; 717 } 718 } 719 apply_adj: 720 #ifdef Avoid_Underflow 721 if (scale && (y = word0(&rv) & Exp_mask) 722 <= 2*P*Exp_msk1) 723 word0(&adj) += (2*P+1)*Exp_msk1 - y; 724 #else 725 #ifdef Sudden_Underflow 726 if ((word0(&rv) & Exp_mask) <= 727 P*Exp_msk1) { 728 word0(&rv) += P*Exp_msk1; 729 dval(&rv) += adj*ulp(&rv); 730 word0(&rv) -= P*Exp_msk1; 731 } 732 else 733 #endif /*Sudden_Underflow*/ 734 #endif /*Avoid_Underflow*/ 735 dval(&rv) += adj.d*ulp(&rv); 736 } 737 break; 738 } 739 dval(&adj) = ratio(delta, bs); 740 if (adj.d < 1.) 741 dval(&adj) = 1.; 742 if (adj.d <= 0x7ffffffe) { 743 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */ 744 y = adj.d; 745 if (y != adj.d) { 746 if (!(((unsigned int)Rounding>>1) ^ (unsigned int)dsign)) 747 y++; 748 dval(&adj) = y; 749 } 750 } 751 #ifdef Avoid_Underflow 752 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 753 word0(&adj) += (2*P+1)*Exp_msk1 - y; 754 #else 755 #ifdef Sudden_Underflow 756 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 757 word0(&rv) += P*Exp_msk1; 758 dval(&adj) *= ulp(&rv); 759 if (dsign) 760 dval(&rv) += adj; 761 else 762 dval(&rv) -= adj; 763 word0(&rv) -= P*Exp_msk1; 764 goto cont; 765 } 766 #endif /*Sudden_Underflow*/ 767 #endif /*Avoid_Underflow*/ 768 dval(&adj) *= ulp(&rv); 769 if (dsign) { 770 if (word0(&rv) == Big0 && word1(&rv) == Big1) 771 goto ovfl; 772 dval(&rv) += adj.d; 773 } 774 else 775 dval(&rv) -= adj.d; 776 goto cont; 777 } 778 #endif /*Honor_FLT_ROUNDS*/ 779 780 if (i < 0) { 781 /* Error is less than half an ulp -- check for 782 * special case of mantissa a power of two. 783 */ 784 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask 785 #ifdef IEEE_Arith 786 #ifdef Avoid_Underflow 787 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 788 #else 789 || (word0(&rv) & Exp_mask) <= Exp_msk1 790 #endif 791 #endif 792 ) { 793 #ifdef SET_INEXACT 794 if (!delta->x[0] && delta->wds <= 1) 795 inexact = 0; 796 #endif 797 break; 798 } 799 if (!delta->x[0] && delta->wds <= 1) { 800 /* exact result */ 801 #ifdef SET_INEXACT 802 inexact = 0; 803 #endif 804 break; 805 } 806 delta = lshift(delta,Log2P); 807 if (cmp(delta, bs) > 0) 808 goto drop_down; 809 break; 810 } 811 if (i == 0) { 812 /* exactly half-way between */ 813 if (dsign) { 814 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 815 && word1(&rv) == ( 816 #ifdef Avoid_Underflow 817 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) 818 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 819 #endif 820 0xffffffff)) { 821 /*boundary case -- increment exponent*/ 822 if (word0(&rv) == Big0 && word1(&rv) == Big1) 823 goto ovfl; 824 word0(&rv) = (word0(&rv) & Exp_mask) 825 + Exp_msk1 826 #ifdef IBM 827 | Exp_msk1 >> 4 828 #endif 829 ; 830 word1(&rv) = 0; 831 #ifdef Avoid_Underflow 832 dsign = 0; 833 #endif 834 break; 835 } 836 } 837 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 838 drop_down: 839 /* boundary case -- decrement exponent */ 840 #ifdef Sudden_Underflow /*{{*/ 841 L = word0(&rv) & Exp_mask; 842 #ifdef IBM 843 if (L < Exp_msk1) 844 #else 845 #ifdef Avoid_Underflow 846 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 847 #else 848 if (L <= Exp_msk1) 849 #endif /*Avoid_Underflow*/ 850 #endif /*IBM*/ 851 goto undfl; 852 L -= Exp_msk1; 853 #else /*Sudden_Underflow}{*/ 854 #ifdef Avoid_Underflow 855 if (scale) { 856 L = word0(&rv) & Exp_mask; 857 if (L <= (2*P+1)*Exp_msk1) { 858 if (L > (P+2)*Exp_msk1) 859 /* round even ==> */ 860 /* accept rv */ 861 break; 862 /* rv = smallest denormal */ 863 goto undfl; 864 } 865 } 866 #endif /*Avoid_Underflow*/ 867 L = (word0(&rv) & Exp_mask) - Exp_msk1; 868 #endif /*Sudden_Underflow}}*/ 869 word0(&rv) = L | Bndry_mask1; 870 word1(&rv) = 0xffffffff; 871 #ifdef IBM 872 goto cont; 873 #else 874 break; 875 #endif 876 } 877 #ifndef ROUND_BIASED 878 #ifdef Avoid_Underflow 879 if (Lsb1) { 880 if (!(word0(&rv) & Lsb1)) 881 break; 882 } 883 else if (!(word1(&rv) & Lsb)) 884 break; 885 #else 886 if (!(word1(&rv) & LSB)) 887 break; 888 #endif 889 #endif 890 if (dsign) 891 #ifdef Avoid_Underflow 892 dval(&rv) += sulp(&rv, scale); 893 #else 894 dval(&rv) += ulp(&rv); 895 #endif 896 #ifndef ROUND_BIASED 897 else { 898 #ifdef Avoid_Underflow 899 dval(&rv) -= sulp(&rv, scale); 900 #else 901 dval(&rv) -= ulp(&rv); 902 #endif 903 #ifndef Sudden_Underflow 904 if (!dval(&rv)) 905 goto undfl; 906 #endif 907 } 908 #ifdef Avoid_Underflow 909 dsign = 1 - dsign; 910 #endif 911 #endif 912 break; 913 } 914 if ((aadj = ratio(delta, bs)) <= 2.) { 915 if (dsign) 916 aadj = dval(&aadj1) = 1.; 917 else if (word1(&rv) || word0(&rv) & Bndry_mask) { 918 #ifndef Sudden_Underflow 919 if (word1(&rv) == Tiny1 && !word0(&rv)) 920 goto undfl; 921 #endif 922 aadj = 1.; 923 dval(&aadj1) = -1.; 924 } 925 else { 926 /* special case -- power of FLT_RADIX to be */ 927 /* rounded down... */ 928 929 if (aadj < 2./FLT_RADIX) 930 aadj = 1./FLT_RADIX; 931 else 932 aadj *= 0.5; 933 dval(&aadj1) = -aadj; 934 } 935 } 936 else { 937 aadj *= 0.5; 938 dval(&aadj1) = dsign ? aadj : -aadj; 939 #ifdef Check_FLT_ROUNDS 940 /* CONSTCOND */ 941 switch(Rounding) { 942 case 2: /* towards +infinity */ 943 dval(&aadj1) -= 0.5; 944 break; 945 case 0: /* towards 0 */ 946 case 3: /* towards -infinity */ 947 dval(&aadj1) += 0.5; 948 } 949 #else 950 /* CONSTCOND */ 951 if (Flt_Rounds == 0) 952 dval(&aadj1) += 0.5; 953 #endif /*Check_FLT_ROUNDS*/ 954 } 955 y = word0(&rv) & Exp_mask; 956 957 /* Check for overflow */ 958 959 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 960 dval(&rv0) = dval(&rv); 961 word0(&rv) -= P*Exp_msk1; 962 dval(&adj) = dval(&aadj1) * ulp(&rv); 963 dval(&rv) += dval(&adj); 964 if ((word0(&rv) & Exp_mask) >= 965 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 966 if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 967 goto ovfl; 968 word0(&rv) = Big0; 969 word1(&rv) = Big1; 970 goto cont; 971 } 972 else 973 word0(&rv) += P*Exp_msk1; 974 } 975 else { 976 #ifdef Avoid_Underflow 977 if (scale && y <= 2*P*Exp_msk1) { 978 if (aadj <= 0x7fffffff) { 979 if ((z = aadj) == 0) 980 z = 1; 981 aadj = z; 982 dval(&aadj1) = dsign ? aadj : -aadj; 983 } 984 word0(&aadj1) += (2*P+1)*Exp_msk1 - y; 985 } 986 dval(&adj) = dval(&aadj1) * ulp(&rv); 987 dval(&rv) += dval(&adj); 988 #else 989 #ifdef Sudden_Underflow 990 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { 991 dval(&rv0) = dval(&rv); 992 word0(&rv) += P*Exp_msk1; 993 dval(&adj) = dval(&aadj1) * ulp(&rv); 994 dval(&rv) += dval(&adj); 995 #ifdef IBM 996 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) 997 #else 998 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) 999 #endif 1000 { 1001 if (word0(&rv0) == Tiny0 1002 && word1(&rv0) == Tiny1) 1003 goto undfl; 1004 word0(&rv) = Tiny0; 1005 word1(&rv) = Tiny1; 1006 goto cont; 1007 } 1008 else 1009 word0(&rv) -= P*Exp_msk1; 1010 } 1011 else { 1012 dval(&adj) = dval(&aadj1) * ulp(&rv); 1013 dval(&rv) += dval(&adj); 1014 } 1015 #else /*Sudden_Underflow*/ 1016 /* Compute dval(&adj) so that the IEEE rounding rules will 1017 * correctly round rv + dval(&adj) in some half-way cases. 1018 * If rv * ulp(&rv) is denormalized (i.e., 1019 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 1020 * trouble from bits lost to denormalization; 1021 * example: 1.2e-307 . 1022 */ 1023 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 1024 dval(&aadj1) = (double)(int)(aadj + 0.5); 1025 if (!dsign) 1026 dval(&aadj1) = -dval(&aadj1); 1027 } 1028 dval(&adj) = dval(&aadj1) * ulp(&rv); 1029 dval(&rv) += adj; 1030 #endif /*Sudden_Underflow*/ 1031 #endif /*Avoid_Underflow*/ 1032 } 1033 z = word0(&rv) & Exp_mask; 1034 #ifndef SET_INEXACT 1035 #ifdef Avoid_Underflow 1036 if (!scale) 1037 #endif 1038 if (y == z) { 1039 /* Can we stop now? */ 1040 L = (Long)aadj; 1041 aadj -= L; 1042 /* The tolerances below are conservative. */ 1043 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 1044 if (aadj < .4999999 || aadj > .5000001) 1045 break; 1046 } 1047 else if (aadj < .4999999/FLT_RADIX) 1048 break; 1049 } 1050 #endif 1051 cont: 1052 Bfree(bb); 1053 Bfree(bd); 1054 Bfree(bs); 1055 Bfree(delta); 1056 } 1057 Bfree(bb); 1058 Bfree(bd); 1059 Bfree(bs); 1060 Bfree(bd0); 1061 Bfree(delta); 1062 #ifdef SET_INEXACT 1063 if (inexact) { 1064 if (!oldinexact) { 1065 word0(&rv0) = Exp_1 + (70 << Exp_shift); 1066 word1(&rv0) = 0; 1067 dval(&rv0) += 1.; 1068 } 1069 } 1070 else if (!oldinexact) 1071 clear_inexact(); 1072 #endif 1073 #ifdef Avoid_Underflow 1074 if (scale) { 1075 word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 1076 word1(&rv0) = 0; 1077 dval(&rv) *= dval(&rv0); 1078 #ifndef NO_ERRNO 1079 /* try to avoid the bug of testing an 8087 register value */ 1080 #ifdef IEEE_Arith 1081 if (!(word0(&rv) & Exp_mask)) 1082 #else 1083 if (word0(&rv) == 0 && word1(&rv) == 0) 1084 #endif 1085 errno = ERANGE; 1086 #endif 1087 } 1088 #endif /* Avoid_Underflow */ 1089 #ifdef SET_INEXACT 1090 if (inexact && !(word0(&rv) & Exp_mask)) { 1091 /* set underflow bit */ 1092 dval(&rv0) = 1e-300; 1093 dval(&rv0) *= dval(&rv0); 1094 } 1095 #endif 1096 ret: 1097 if (se) 1098 *se = __UNCONST(s); 1099 return sign ? -dval(&rv) : dval(&rv); 1100 } 1101 1102 double 1103 strtod(CONST char *s, char **sp) 1104 { 1105 return _int_strtod_l(s, sp, _current_locale()); 1106 } 1107 1108 #ifdef __weak_alias 1109 __weak_alias(strtod_l, _strtod_l) 1110 #endif 1111 1112 double 1113 strtod_l(CONST char *s, char **sp, locale_t loc) 1114 { 1115 return _int_strtod_l(s, sp, loc); 1116 } 1117