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