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