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