1 /* $NetBSD: strtod.c,v 1.15 2019/08/01 02:27:43 riastradh 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 #ifdef INFNAN_CHECK 100 int decpt; 101 #endif 102 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, 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 #ifdef INFNAN_CHECK 137 decpt = 0; 138 #endif 139 sign = nz0 = nz = 0; 140 dval(&rv) = 0.; 141 for(s = s00;;s++) switch(*s) { 142 case '-': 143 sign = 1; 144 /* FALLTHROUGH */ 145 case '+': 146 if (*++s) 147 goto break2; 148 /* FALLTHROUGH */ 149 case 0: 150 goto ret0; 151 case '\t': 152 case '\n': 153 case '\v': 154 case '\f': 155 case '\r': 156 case ' ': 157 continue; 158 default: 159 goto break2; 160 } 161 break2: 162 if (*s == '0') { 163 #ifndef NO_HEX_FP /*{*/ 164 { 165 static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 166 Long expt; 167 ULong bits[2]; 168 switch(s[1]) { 169 case 'x': 170 case 'X': 171 { 172 #ifdef Honor_FLT_ROUNDS 173 FPI fpi1 = fpi; 174 fpi1.rounding = Rounding; 175 #else 176 #define fpi1 fpi 177 #endif 178 switch((i = gethex(&s, &fpi1, &expt, &bb, sign, loc)) & STRTOG_Retmask) { 179 case STRTOG_NoNumber: 180 s = s00; 181 sign = 0; 182 /* FALLTHROUGH */ 183 case STRTOG_Zero: 184 break; 185 default: 186 if (bb) { 187 copybits(bits, fpi.nbits, bb); 188 Bfree(bb); 189 } 190 ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i); 191 }} 192 goto ret; 193 } 194 } 195 #endif /*}*/ 196 nz0 = 1; 197 while(*++s == '0') ; 198 if (!*s) 199 goto ret; 200 } 201 s0 = s; 202 y = z = 0; 203 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 204 if (nd < 9) 205 y = 10*y + c - '0'; 206 else if (nd < 16) 207 z = 10*z + c - '0'; 208 nd0 = nd; 209 #ifdef USE_LOCALE 210 if (c == *decimalpoint) { 211 for(i = 1; decimalpoint[i]; ++i) 212 if (s[i] != decimalpoint[i]) 213 goto dig_done; 214 s += i; 215 c = *s; 216 #else 217 if (c == '.') { 218 c = *++s; 219 #endif 220 #ifdef INFNAN_CHECK 221 decpt = 1; 222 #endif 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 /* FALLTHROUGH */ 264 case '+': 265 c = *++s; 266 } 267 if (c >= '0' && c <= '9') { 268 while(c == '0') 269 c = *++s; 270 if (c > '0' && c <= '9') { 271 L = c - '0'; 272 s1 = s; 273 while((c = *++s) >= '0' && c <= '9') 274 L = 10*L + c - '0'; 275 if (s - s1 > 8 || L > 19999) 276 /* Avoid confusion from exponents 277 * so large that e might overflow. 278 */ 279 e = 19999; /* safe for 16 bit ints */ 280 else 281 e = (int)L; 282 if (esign) 283 e = -e; 284 } 285 else 286 e = 0; 287 } 288 else 289 s = s00; 290 } 291 if (!nd) { 292 if (!nz && !nz0) { 293 #ifdef INFNAN_CHECK 294 /* Check for Nan and Infinity */ 295 ULong bits[2]; 296 static CONST FPI fpinan = /* only 52 explicit bits */ 297 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 298 if (!decpt) 299 switch(c) { 300 case 'i': 301 case 'I': 302 if (match(&s,"nf")) { 303 --s; 304 if (!match(&s,"inity")) 305 ++s; 306 word0(&rv) = 0x7ff00000; 307 word1(&rv) = 0; 308 goto ret; 309 } 310 break; 311 case 'n': 312 case 'N': 313 if (match(&s, "an")) { 314 #ifndef No_Hex_NaN 315 if (*s == '(' /*)*/ 316 && hexnan(&s, &fpinan, bits) 317 == STRTOG_NaNbits) { 318 word0(&rv) = 0x7ff00000 | bits[1]; 319 word1(&rv) = bits[0]; 320 } 321 else { 322 #endif 323 word0(&rv) = NAN_WORD0; 324 word1(&rv) = NAN_WORD1; 325 #ifndef No_Hex_NaN 326 } 327 #endif 328 goto ret; 329 } 330 } 331 #endif /* INFNAN_CHECK */ 332 ret0: 333 s = s00; 334 sign = 0; 335 } 336 goto ret; 337 } 338 e1 = e -= nf; 339 340 /* Now we have nd0 digits, starting at s0, followed by a 341 * decimal point, followed by nd-nd0 digits. The number we're 342 * after is the integer represented by those digits times 343 * 10**e */ 344 345 if (!nd0) 346 nd0 = nd; 347 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 348 dval(&rv) = y; 349 if (k > 9) { 350 #ifdef SET_INEXACT 351 if (k > DBL_DIG) 352 oldinexact = get_inexact(); 353 #endif 354 dval(&rv) = tens[k - 9] * dval(&rv) + z; 355 } 356 bd0 = 0; 357 if (nd <= DBL_DIG 358 #ifndef RND_PRODQUOT 359 #ifndef Honor_FLT_ROUNDS 360 && Flt_Rounds == 1 361 #endif 362 #endif 363 ) { 364 if (!e) 365 goto ret; 366 #ifndef ROUND_BIASED_without_Round_Up 367 if (e > 0) { 368 if (e <= Ten_pmax) { 369 #ifdef VAX 370 goto vax_ovfl_check; 371 #else 372 #ifdef Honor_FLT_ROUNDS 373 /* round correctly FLT_ROUNDS = 2 or 3 */ 374 if (sign) { 375 rv.d = -rv.d; 376 sign = 0; 377 } 378 #endif 379 /* rv = */ rounded_product(dval(&rv), tens[e]); 380 goto ret; 381 #endif 382 } 383 i = DBL_DIG - nd; 384 if (e <= Ten_pmax + i) { 385 /* A fancier test would sometimes let us do 386 * this for larger i values. 387 */ 388 #ifdef Honor_FLT_ROUNDS 389 /* round correctly FLT_ROUNDS = 2 or 3 */ 390 if (sign) { 391 rv.d = -rv.d; 392 sign = 0; 393 } 394 #endif 395 e -= i; 396 dval(&rv) *= tens[i]; 397 #ifdef VAX 398 /* VAX exponent range is so narrow we must 399 * worry about overflow here... 400 */ 401 vax_ovfl_check: 402 word0(&rv) -= P*Exp_msk1; 403 /* rv = */ rounded_product(dval(&rv), tens[e]); 404 if ((word0(&rv) & Exp_mask) 405 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 406 goto ovfl; 407 word0(&rv) += P*Exp_msk1; 408 #else 409 /* rv = */ rounded_product(dval(&rv), tens[e]); 410 #endif 411 goto ret; 412 } 413 } 414 #ifndef Inaccurate_Divide 415 else if (e >= -Ten_pmax) { 416 #ifdef Honor_FLT_ROUNDS 417 /* round correctly FLT_ROUNDS = 2 or 3 */ 418 if (sign) { 419 rv.d = -rv.d; 420 sign = 0; 421 } 422 #endif 423 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 424 goto ret; 425 } 426 #endif 427 #endif /* ROUND_BIASED_without_Round_Up */ 428 } 429 e1 += nd - k; 430 431 #ifdef IEEE_Arith 432 #ifdef SET_INEXACT 433 inexact = 1; 434 if (k <= DBL_DIG) 435 oldinexact = get_inexact(); 436 #endif 437 #ifdef Avoid_Underflow 438 scale = 0; 439 #endif 440 #ifdef Honor_FLT_ROUNDS 441 if (Rounding >= 2) { 442 if (sign) 443 Rounding = Rounding == 2 ? 0 : 2; 444 else 445 if (Rounding != 2) 446 Rounding = 0; 447 } 448 #endif 449 #endif /*IEEE_Arith*/ 450 451 /* Get starting approximation = rv * 10**e1 */ 452 453 if (e1 > 0) { 454 if ( (i = e1 & 15) !=0) 455 dval(&rv) *= tens[i]; 456 if (e1 &= ~15) { 457 if (e1 > DBL_MAX_10_EXP) { 458 ovfl: 459 /* Can't trust HUGE_VAL */ 460 #ifdef IEEE_Arith 461 #ifdef Honor_FLT_ROUNDS 462 switch(Rounding) { 463 case 0: /* toward 0 */ 464 case 3: /* toward -infinity */ 465 word0(&rv) = Big0; 466 word1(&rv) = Big1; 467 break; 468 default: 469 word0(&rv) = Exp_mask; 470 word1(&rv) = 0; 471 } 472 #else /*Honor_FLT_ROUNDS*/ 473 word0(&rv) = Exp_mask; 474 word1(&rv) = 0; 475 #endif /*Honor_FLT_ROUNDS*/ 476 #ifdef SET_INEXACT 477 /* set overflow bit */ 478 dval(&rv0) = 1e300; 479 dval(&rv0) *= dval(&rv0); 480 #endif 481 #else /*IEEE_Arith*/ 482 word0(&rv) = Big0; 483 word1(&rv) = Big1; 484 #endif /*IEEE_Arith*/ 485 range_err: 486 if (bd0) { 487 Bfree(bb); 488 Bfree(bd); 489 Bfree(bs); 490 Bfree(bd0); 491 Bfree(delta); 492 } 493 #ifndef NO_ERRNO 494 errno = ERANGE; 495 #endif 496 goto ret; 497 } 498 e1 = (unsigned int)e1 >> 4; 499 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1) 500 if (e1 & 1) 501 dval(&rv) *= bigtens[j]; 502 /* The last multiplication could overflow. */ 503 word0(&rv) -= P*Exp_msk1; 504 dval(&rv) *= bigtens[j]; 505 if ((z = word0(&rv) & Exp_mask) 506 > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 507 goto ovfl; 508 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 509 /* set to largest number */ 510 /* (Can't trust DBL_MAX) */ 511 word0(&rv) = Big0; 512 word1(&rv) = Big1; 513 } 514 else 515 word0(&rv) += P*Exp_msk1; 516 } 517 } 518 else if (e1 < 0) { 519 e1 = -e1; 520 if ( (i = e1 & 15) !=0) 521 dval(&rv) /= tens[i]; 522 if (e1 >>= 4) { 523 if (e1 >= 1 << n_bigtens) 524 goto undfl; 525 #ifdef Avoid_Underflow 526 if (e1 & Scale_Bit) 527 scale = 2*P; 528 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1) 529 if (e1 & 1) 530 dval(&rv) *= tinytens[j]; 531 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 532 >> Exp_shift)) > 0) { 533 /* scaled rv is denormal; zap j low bits */ 534 if (j >= 32) { 535 word1(&rv) = 0; 536 if (j >= 53) 537 word0(&rv) = (P+2)*Exp_msk1; 538 else 539 word0(&rv) &= 0xffffffffU << (j-32); 540 } 541 else 542 word1(&rv) &= 0xffffffffU << j; 543 } 544 #else 545 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1) 546 if (e1 & 1) 547 dval(&rv) *= tinytens[j]; 548 /* The last multiplication could underflow. */ 549 dval(&rv0) = dval(&rv); 550 dval(&rv) *= tinytens[j]; 551 if (!dval(&rv)) { 552 dval(&rv) = 2.*dval(&rv0); 553 dval(&rv) *= tinytens[j]; 554 #endif 555 if (!dval(&rv)) { 556 undfl: 557 dval(&rv) = 0.; 558 goto range_err; 559 } 560 #ifndef Avoid_Underflow 561 word0(&rv) = Tiny0; 562 word1(&rv) = Tiny1; 563 /* The refinement below will clean 564 * this approximation up. 565 */ 566 } 567 #endif 568 } 569 } 570 571 /* Now the hard part -- adjusting rv to the correct value.*/ 572 573 /* Put digits into bd: true value = bd * 10^e */ 574 575 bd0 = s2b(s0, nd0, nd, y, dplen); 576 if (bd0 == NULL) 577 goto ovfl; 578 579 for(;;) { 580 bd = Balloc(bd0->k); 581 if (bd == NULL) 582 goto ovfl; 583 Bcopy(bd, bd0); 584 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 585 if (bb == NULL) 586 goto ovfl; 587 bs = i2b(1); 588 if (bs == NULL) 589 goto ovfl; 590 591 if (e >= 0) { 592 bb2 = bb5 = 0; 593 bd2 = bd5 = e; 594 } 595 else { 596 bb2 = bb5 = -e; 597 bd2 = bd5 = 0; 598 } 599 if (bbe >= 0) 600 bb2 += bbe; 601 else 602 bd2 -= bbe; 603 bs2 = bb2; 604 #ifdef Honor_FLT_ROUNDS 605 if (Rounding != 1) 606 bs2++; 607 #endif 608 #ifdef Avoid_Underflow 609 Lsb = LSB; 610 Lsb1 = 0; 611 j = bbe - scale; 612 i = j + bbbits - 1; /* logb(rv) */ 613 j = P + 1 - bbbits; 614 if (i < Emin) { /* denormal */ 615 i = Emin - i; 616 j -= i; 617 if (i < 32) 618 Lsb <<= i; 619 else 620 Lsb1 = Lsb << (i-32); 621 } 622 #else /*Avoid_Underflow*/ 623 #ifdef Sudden_Underflow 624 #ifdef IBM 625 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 626 #else 627 j = P + 1 - bbbits; 628 #endif 629 #else /*Sudden_Underflow*/ 630 j = bbe; 631 i = j + bbbits - 1; /* logb(&rv) */ 632 if (i < Emin) /* denormal */ 633 j += P - Emin; 634 else 635 j = P + 1 - bbbits; 636 #endif /*Sudden_Underflow*/ 637 #endif /*Avoid_Underflow*/ 638 bb2 += j; 639 bd2 += j; 640 #ifdef Avoid_Underflow 641 bd2 += scale; 642 #endif 643 i = bb2 < bd2 ? bb2 : bd2; 644 if (i > bs2) 645 i = bs2; 646 if (i > 0) { 647 bb2 -= i; 648 bd2 -= i; 649 bs2 -= i; 650 } 651 if (bb5 > 0) { 652 bs = pow5mult(bs, bb5); 653 if (bs == NULL) 654 goto ovfl; 655 bb1 = mult(bs, bb); 656 if (bb1 == NULL) 657 goto ovfl; 658 Bfree(bb); 659 bb = bb1; 660 } 661 if (bb2 > 0) { 662 bb = lshift(bb, bb2); 663 if (bb == NULL) 664 goto ovfl; 665 } 666 if (bd5 > 0) { 667 bd = pow5mult(bd, bd5); 668 if (bd == NULL) 669 goto ovfl; 670 } 671 if (bd2 > 0) { 672 bd = lshift(bd, bd2); 673 if (bd == NULL) 674 goto ovfl; 675 } 676 if (bs2 > 0) { 677 bs = lshift(bs, bs2); 678 if (bs == NULL) 679 goto ovfl; 680 } 681 delta = diff(bb, bd); 682 if (delta == NULL) 683 goto ovfl; 684 dsign = delta->sign; 685 delta->sign = 0; 686 i = cmp(delta, bs); 687 #ifdef Honor_FLT_ROUNDS 688 if (Rounding != 1) { 689 if (i < 0) { 690 /* Error is less than an ulp */ 691 if (!delta->x[0] && delta->wds <= 1) { 692 /* exact */ 693 #ifdef SET_INEXACT 694 inexact = 0; 695 #endif 696 break; 697 } 698 if (Rounding) { 699 if (dsign) { 700 dval(&adj) = 1.; 701 goto apply_adj; 702 } 703 } 704 else if (!dsign) { 705 dval(&adj) = -1.; 706 if (!word1(&rv) 707 && !(word0(&rv) & Frac_mask)) { 708 y = word0(&rv) & Exp_mask; 709 #ifdef Avoid_Underflow 710 if (!scale || y > 2*P*Exp_msk1) 711 #else 712 if (y) 713 #endif 714 { 715 delta = lshift(delta,Log2P); 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 (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