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