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