1 /* $NetBSD: strtod.c,v 1.5 2008/03/21 23:13:48 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 "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 = NULL, *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 if (bd0 == NULL) 527 goto ovfl; 528 529 for(;;) { 530 bd = Balloc(bd0->k); 531 if (bd == NULL) 532 goto ovfl; 533 Bcopy(bd, bd0); 534 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ 535 if (bb == NULL) 536 goto ovfl; 537 bs = i2b(1); 538 if (bs == NULL) 539 goto ovfl; 540 541 if (e >= 0) { 542 bb2 = bb5 = 0; 543 bd2 = bd5 = e; 544 } 545 else { 546 bb2 = bb5 = -e; 547 bd2 = bd5 = 0; 548 } 549 if (bbe >= 0) 550 bb2 += bbe; 551 else 552 bd2 -= bbe; 553 bs2 = bb2; 554 #ifdef Honor_FLT_ROUNDS 555 if (rounding != 1) 556 bs2++; 557 #endif 558 #ifdef Avoid_Underflow 559 j = bbe - scale; 560 i = j + bbbits - 1; /* logb(rv) */ 561 if (i < Emin) /* denormal */ 562 j += P - Emin; 563 else 564 j = P + 1 - bbbits; 565 #else /*Avoid_Underflow*/ 566 #ifdef Sudden_Underflow 567 #ifdef IBM 568 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 569 #else 570 j = P + 1 - bbbits; 571 #endif 572 #else /*Sudden_Underflow*/ 573 j = bbe; 574 i = j + bbbits - 1; /* logb(rv) */ 575 if (i < Emin) /* denormal */ 576 j += P - Emin; 577 else 578 j = P + 1 - bbbits; 579 #endif /*Sudden_Underflow*/ 580 #endif /*Avoid_Underflow*/ 581 bb2 += j; 582 bd2 += j; 583 #ifdef Avoid_Underflow 584 bd2 += scale; 585 #endif 586 i = bb2 < bd2 ? bb2 : bd2; 587 if (i > bs2) 588 i = bs2; 589 if (i > 0) { 590 bb2 -= i; 591 bd2 -= i; 592 bs2 -= i; 593 } 594 if (bb5 > 0) { 595 bs = pow5mult(bs, bb5); 596 if (bs == NULL) 597 goto ovfl; 598 bb1 = mult(bs, bb); 599 if (bb1 == NULL) 600 goto ovfl; 601 Bfree(bb); 602 bb = bb1; 603 } 604 if (bb2 > 0) { 605 bb = lshift(bb, bb2); 606 if (bb == NULL) 607 goto ovfl; 608 } 609 if (bd5 > 0) { 610 bd = pow5mult(bd, bd5); 611 if (bd == NULL) 612 goto ovfl; 613 } 614 if (bd2 > 0) { 615 bd = lshift(bd, bd2); 616 if (bd == NULL) 617 goto ovfl; 618 } 619 if (bs2 > 0) { 620 bs = lshift(bs, bs2); 621 if (bs == NULL) 622 goto ovfl; 623 } 624 delta = diff(bb, bd); 625 if (delta == NULL) 626 goto ovfl; 627 dsign = delta->sign; 628 delta->sign = 0; 629 i = cmp(delta, bs); 630 #ifdef Honor_FLT_ROUNDS 631 if (rounding != 1) { 632 if (i < 0) { 633 /* Error is less than an ulp */ 634 if (!delta->x[0] && delta->wds <= 1) { 635 /* exact */ 636 #ifdef SET_INEXACT 637 inexact = 0; 638 #endif 639 break; 640 } 641 if (rounding) { 642 if (dsign) { 643 adj = 1.; 644 goto apply_adj; 645 } 646 } 647 else if (!dsign) { 648 adj = -1.; 649 if (!word1(rv) 650 && !(word0(rv) & Frac_mask)) { 651 y = word0(rv) & Exp_mask; 652 #ifdef Avoid_Underflow 653 if (!scale || y > 2*P*Exp_msk1) 654 #else 655 if (y) 656 #endif 657 { 658 delta = lshift(delta,Log2P); 659 if (cmp(delta, bs) <= 0) 660 adj = -0.5; 661 } 662 } 663 apply_adj: 664 #ifdef Avoid_Underflow 665 if (scale && (y = word0(rv) & Exp_mask) 666 <= 2*P*Exp_msk1) 667 word0(adj) += (2*P+1)*Exp_msk1 - y; 668 #else 669 #ifdef Sudden_Underflow 670 if ((word0(rv) & Exp_mask) <= 671 P*Exp_msk1) { 672 word0(rv) += P*Exp_msk1; 673 dval(rv) += adj*ulp(dval(rv)); 674 word0(rv) -= P*Exp_msk1; 675 } 676 else 677 #endif /*Sudden_Underflow*/ 678 #endif /*Avoid_Underflow*/ 679 dval(rv) += adj*ulp(dval(rv)); 680 } 681 break; 682 } 683 adj = ratio(delta, bs); 684 if (adj < 1.) 685 adj = 1.; 686 if (adj <= 0x7ffffffe) { 687 /* adj = rounding ? ceil(adj) : floor(adj); */ 688 y = adj; 689 if (y != adj) { 690 if (!((rounding>>1) ^ dsign)) 691 y++; 692 adj = y; 693 } 694 } 695 #ifdef Avoid_Underflow 696 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 697 word0(adj) += (2*P+1)*Exp_msk1 - y; 698 #else 699 #ifdef Sudden_Underflow 700 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 701 word0(rv) += P*Exp_msk1; 702 adj *= ulp(dval(rv)); 703 if (dsign) 704 dval(rv) += adj; 705 else 706 dval(rv) -= adj; 707 word0(rv) -= P*Exp_msk1; 708 goto cont; 709 } 710 #endif /*Sudden_Underflow*/ 711 #endif /*Avoid_Underflow*/ 712 adj *= ulp(dval(rv)); 713 if (dsign) 714 dval(rv) += adj; 715 else 716 dval(rv) -= adj; 717 goto cont; 718 } 719 #endif /*Honor_FLT_ROUNDS*/ 720 721 if (i < 0) { 722 /* Error is less than half an ulp -- check for 723 * special case of mantissa a power of two. 724 */ 725 if (dsign || word1(rv) || word0(rv) & Bndry_mask 726 #ifdef IEEE_Arith 727 #ifdef Avoid_Underflow 728 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 729 #else 730 || (word0(rv) & Exp_mask) <= Exp_msk1 731 #endif 732 #endif 733 ) { 734 #ifdef SET_INEXACT 735 if (!delta->x[0] && delta->wds <= 1) 736 inexact = 0; 737 #endif 738 break; 739 } 740 if (!delta->x[0] && delta->wds <= 1) { 741 /* exact result */ 742 #ifdef SET_INEXACT 743 inexact = 0; 744 #endif 745 break; 746 } 747 delta = lshift(delta,Log2P); 748 if (cmp(delta, bs) > 0) 749 goto drop_down; 750 break; 751 } 752 if (i == 0) { 753 /* exactly half-way between */ 754 if (dsign) { 755 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 756 && word1(rv) == ( 757 #ifdef Avoid_Underflow 758 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1) 759 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 760 #endif 761 0xffffffff)) { 762 /*boundary case -- increment exponent*/ 763 word0(rv) = (word0(rv) & Exp_mask) 764 + Exp_msk1 765 #ifdef IBM 766 | Exp_msk1 >> 4 767 #endif 768 ; 769 word1(rv) = 0; 770 #ifdef Avoid_Underflow 771 dsign = 0; 772 #endif 773 break; 774 } 775 } 776 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 777 drop_down: 778 /* boundary case -- decrement exponent */ 779 #ifdef Sudden_Underflow /*{{*/ 780 L = word0(rv) & Exp_mask; 781 #ifdef IBM 782 if (L < Exp_msk1) 783 #else 784 #ifdef Avoid_Underflow 785 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) 786 #else 787 if (L <= Exp_msk1) 788 #endif /*Avoid_Underflow*/ 789 #endif /*IBM*/ 790 goto undfl; 791 L -= Exp_msk1; 792 #else /*Sudden_Underflow}{*/ 793 #ifdef Avoid_Underflow 794 if (scale) { 795 L = word0(rv) & Exp_mask; 796 if (L <= (2*P+1)*Exp_msk1) { 797 if (L > (P+2)*Exp_msk1) 798 /* round even ==> */ 799 /* accept rv */ 800 break; 801 /* rv = smallest denormal */ 802 goto undfl; 803 } 804 } 805 #endif /*Avoid_Underflow*/ 806 L = (word0(rv) & Exp_mask) - Exp_msk1; 807 #endif /*Sudden_Underflow}*/ 808 word0(rv) = L | Bndry_mask1; 809 word1(rv) = 0xffffffff; 810 #ifdef IBM 811 goto cont; 812 #else 813 break; 814 #endif 815 } 816 #ifndef ROUND_BIASED 817 if (!(word1(rv) & LSB)) 818 break; 819 #endif 820 if (dsign) 821 dval(rv) += ulp(dval(rv)); 822 #ifndef ROUND_BIASED 823 else { 824 dval(rv) -= ulp(dval(rv)); 825 #ifndef Sudden_Underflow 826 if (!dval(rv)) 827 goto undfl; 828 #endif 829 } 830 #ifdef Avoid_Underflow 831 dsign = 1 - dsign; 832 #endif 833 #endif 834 break; 835 } 836 if ((aadj = ratio(delta, bs)) <= 2.) { 837 if (dsign) 838 aadj = aadj1 = 1.; 839 else if (word1(rv) || word0(rv) & Bndry_mask) { 840 #ifndef Sudden_Underflow 841 if (word1(rv) == Tiny1 && !word0(rv)) 842 goto undfl; 843 #endif 844 aadj = 1.; 845 aadj1 = -1.; 846 } 847 else { 848 /* special case -- power of FLT_RADIX to be */ 849 /* rounded down... */ 850 851 if (aadj < 2./FLT_RADIX) 852 aadj = 1./FLT_RADIX; 853 else 854 aadj *= 0.5; 855 aadj1 = -aadj; 856 } 857 } 858 else { 859 aadj *= 0.5; 860 aadj1 = dsign ? aadj : -aadj; 861 #ifdef Check_FLT_ROUNDS 862 switch(Rounding) { 863 case 2: /* towards +infinity */ 864 aadj1 -= 0.5; 865 break; 866 case 0: /* towards 0 */ 867 case 3: /* towards -infinity */ 868 aadj1 += 0.5; 869 } 870 #else 871 if (Flt_Rounds == 0) 872 aadj1 += 0.5; 873 #endif /*Check_FLT_ROUNDS*/ 874 } 875 y = word0(rv) & Exp_mask; 876 877 /* Check for overflow */ 878 879 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 880 dval(rv0) = dval(rv); 881 word0(rv) -= P*Exp_msk1; 882 adj = aadj1 * ulp(dval(rv)); 883 dval(rv) += adj; 884 if ((word0(rv) & Exp_mask) >= 885 Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 886 if (word0(rv0) == Big0 && word1(rv0) == Big1) 887 goto ovfl; 888 word0(rv) = Big0; 889 word1(rv) = Big1; 890 goto cont; 891 } 892 else 893 word0(rv) += P*Exp_msk1; 894 } 895 else { 896 #ifdef Avoid_Underflow 897 if (scale && y <= 2*P*Exp_msk1) { 898 if (aadj <= 0x7fffffff) { 899 if ((z = aadj) == 0) 900 z = 1; 901 aadj = z; 902 aadj1 = dsign ? aadj : -aadj; 903 } 904 word0(aadj1) += (2*P+1)*Exp_msk1 - y; 905 } 906 adj = aadj1 * ulp(dval(rv)); 907 dval(rv) += adj; 908 #else 909 #ifdef Sudden_Underflow 910 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { 911 dval(rv0) = dval(rv); 912 word0(rv) += P*Exp_msk1; 913 adj = aadj1 * ulp(dval(rv)); 914 dval(rv) += adj; 915 #ifdef IBM 916 if ((word0(rv) & Exp_mask) < P*Exp_msk1) 917 #else 918 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) 919 #endif 920 { 921 if (word0(rv0) == Tiny0 922 && word1(rv0) == Tiny1) 923 goto undfl; 924 word0(rv) = Tiny0; 925 word1(rv) = Tiny1; 926 goto cont; 927 } 928 else 929 word0(rv) -= P*Exp_msk1; 930 } 931 else { 932 adj = aadj1 * ulp(dval(rv)); 933 dval(rv) += adj; 934 } 935 #else /*Sudden_Underflow*/ 936 /* Compute adj so that the IEEE rounding rules will 937 * correctly round rv + adj in some half-way cases. 938 * If rv * ulp(rv) is denormalized (i.e., 939 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 940 * trouble from bits lost to denormalization; 941 * example: 1.2e-307 . 942 */ 943 if (y <= (P-1)*Exp_msk1 && aadj > 1.) { 944 aadj1 = (double)(int)(aadj + 0.5); 945 if (!dsign) 946 aadj1 = -aadj1; 947 } 948 adj = aadj1 * ulp(dval(rv)); 949 dval(rv) += adj; 950 #endif /*Sudden_Underflow*/ 951 #endif /*Avoid_Underflow*/ 952 } 953 z = word0(rv) & Exp_mask; 954 #ifndef SET_INEXACT 955 #ifdef Avoid_Underflow 956 if (!scale) 957 #endif 958 if (y == z) { 959 /* Can we stop now? */ 960 L = (Long)aadj; 961 aadj -= L; 962 /* The tolerances below are conservative. */ 963 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 964 if (aadj < .4999999 || aadj > .5000001) 965 break; 966 } 967 else if (aadj < .4999999/FLT_RADIX) 968 break; 969 } 970 #endif 971 cont: 972 Bfree(bb); 973 Bfree(bd); 974 Bfree(bs); 975 Bfree(delta); 976 } 977 #ifdef SET_INEXACT 978 if (inexact) { 979 if (!oldinexact) { 980 word0(rv0) = Exp_1 + (70 << Exp_shift); 981 word1(rv0) = 0; 982 dval(rv0) += 1.; 983 } 984 } 985 else if (!oldinexact) 986 clear_inexact(); 987 #endif 988 #ifdef Avoid_Underflow 989 if (scale) { 990 word0(rv0) = Exp_1 - 2*P*Exp_msk1; 991 word1(rv0) = 0; 992 dval(rv) *= dval(rv0); 993 #ifndef NO_ERRNO 994 /* try to avoid the bug of testing an 8087 register value */ 995 if (word0(rv) == 0 && word1(rv) == 0) 996 errno = ERANGE; 997 #endif 998 } 999 #endif /* Avoid_Underflow */ 1000 #ifdef SET_INEXACT 1001 if (inexact && !(word0(rv) & Exp_mask)) { 1002 /* set underflow bit */ 1003 dval(rv0) = 1e-300; 1004 dval(rv0) *= dval(rv0); 1005 } 1006 #endif 1007 retfree: 1008 Bfree(bb); 1009 Bfree(bd); 1010 Bfree(bs); 1011 Bfree(bd0); 1012 Bfree(delta); 1013 ret: 1014 if (se) 1015 *se = __UNCONST(s); 1016 return sign ? -dval(rv) : dval(rv); 1017 } 1018 1019