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