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