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