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