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