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