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