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