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