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