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