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