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