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