1 /**************************************************************** 2 3 The author of this software is David M. Gay. 4 5 Copyright (C) 1998-2001 by Lucent Technologies 6 All Rights Reserved 7 8 Permission to use, copy, modify, and distribute this software and 9 its documentation for any purpose and without fee is hereby 10 granted, provided that the above copyright notice appear in all 11 copies and that both that the copyright notice and this 12 permission notice and warranty disclaimer appear in supporting 13 documentation, and that the name of Lucent or any of its entities 14 not be used in advertising or publicity pertaining to 15 distribution of the software without specific, written prior 16 permission. 17 18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25 THIS 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 39 fivesbits[] = { 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 49 increment(b) Bigint *b; 50 #else 51 increment(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 96 decrement(b) Bigint *b; 97 #else 98 decrement(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 128 all_on(b, n) Bigint *b; int n; 129 #else 130 all_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 147 set_ones(b, n) Bigint *b; int n; 148 #else 149 set_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 176 rvOK 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 315 mantbits(d) U *d; 316 #else 317 mantbits(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 337 strtodg 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 decimalpoint_cache = strdup(s0); 367 dplen = strlen(s0); 368 } 369 decimalpoint = (char*)s0; 370 #endif /*NO_LOCALE_CACHE*/ 371 #else /*USE_LOCALE}{*/ 372 #define dplen 1 373 #endif /*USE_LOCALE}}*/ 374 375 irv = STRTOG_Zero; 376 denorm = sign = nz0 = nz = 0; 377 dval(&rv) = 0.; 378 rvb = 0; 379 nbits = fpi->nbits; 380 for(s = s00;;s++) switch(*s) { 381 case '-': 382 sign = 1; 383 /* no break */ 384 case '+': 385 if (*++s) 386 goto break2; 387 /* no break */ 388 case 0: 389 sign = 0; 390 irv = STRTOG_NoNumber; 391 s = s00; 392 goto ret; 393 case '\t': 394 case '\n': 395 case '\v': 396 case '\f': 397 case '\r': 398 case ' ': 399 continue; 400 default: 401 goto break2; 402 } 403 break2: 404 if (*s == '0') { 405 #ifndef NO_HEX_FP 406 switch(s[1]) { 407 case 'x': 408 case 'X': 409 irv = gethex(&s, fpi, exp, &rvb, sign); 410 if (irv == STRTOG_NoMemory) 411 return (STRTOG_NoMemory); 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 549 irv = STRTOG_Normal; 550 e1 = e -= nf; 551 rd = 0; 552 switch(fpi->rounding & 3) { 553 case FPI_Round_up: 554 rd = 2 - sign; 555 break; 556 case FPI_Round_zero: 557 rd = 1; 558 break; 559 case FPI_Round_down: 560 rd = 1 + sign; 561 } 562 563 /* Now we have nd0 digits, starting at s0, followed by a 564 * decimal point, followed by nd-nd0 digits. The number we're 565 * after is the integer represented by those digits times 566 * 10**e */ 567 568 if (!nd0) 569 nd0 = nd; 570 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 571 dval(&rv) = y; 572 if (k > 9) 573 dval(&rv) = tens[k - 9] * dval(&rv) + z; 574 bd0 = 0; 575 if (nbits <= P && nd <= DBL_DIG) { 576 if (!e) { 577 if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv)) { 578 if (irv == STRTOG_NoMemory) 579 return (STRTOG_NoMemory); 580 goto ret; 581 } 582 } 583 else if (e > 0) { 584 if (e <= Ten_pmax) { 585 #ifdef VAX 586 goto vax_ovfl_check; 587 #else 588 i = fivesbits[e] + mantbits(&rv) <= P; 589 /* rv = */ rounded_product(dval(&rv), tens[e]); 590 if (rvOK(&rv, fpi, exp, bits, i, rd, &irv)) { 591 if (irv == STRTOG_NoMemory) 592 return (STRTOG_NoMemory); 593 goto ret; 594 } 595 e1 -= e; 596 goto rv_notOK; 597 #endif 598 } 599 i = DBL_DIG - nd; 600 if (e <= Ten_pmax + i) { 601 /* A fancier test would sometimes let us do 602 * this for larger i values. 603 */ 604 e2 = e - i; 605 e1 -= i; 606 dval(&rv) *= tens[i]; 607 #ifdef VAX 608 /* VAX exponent range is so narrow we must 609 * worry about overflow here... 610 */ 611 vax_ovfl_check: 612 dval(&adj) = dval(&rv); 613 word0(&adj) -= P*Exp_msk1; 614 /* adj = */ rounded_product(dval(&adj), tens[e2]); 615 if ((word0(&adj) & Exp_mask) 616 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 617 goto rv_notOK; 618 word0(&adj) += P*Exp_msk1; 619 dval(&rv) = dval(&adj); 620 #else 621 /* rv = */ rounded_product(dval(&rv), tens[e2]); 622 #endif 623 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv)) { 624 if (irv == STRTOG_NoMemory) 625 return (STRTOG_NoMemory); 626 goto ret; 627 } 628 e1 -= e2; 629 } 630 } 631 #ifndef Inaccurate_Divide 632 else if (e >= -Ten_pmax) { 633 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); 634 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv)) { 635 if (irv == STRTOG_NoMemory) 636 return (STRTOG_NoMemory); 637 goto ret; 638 } 639 e1 -= e; 640 } 641 #endif 642 } 643 rv_notOK: 644 e1 += nd - k; 645 646 /* Get starting approximation = rv * 10**e1 */ 647 648 e2 = 0; 649 if (e1 > 0) { 650 if ( (i = e1 & 15) !=0) 651 dval(&rv) *= tens[i]; 652 if (e1 &= ~15) { 653 e1 >>= 4; 654 while(e1 >= (1 << (n_bigtens-1))) { 655 e2 += ((word0(&rv) & Exp_mask) 656 >> Exp_shift1) - Bias; 657 word0(&rv) &= ~Exp_mask; 658 word0(&rv) |= Bias << Exp_shift1; 659 dval(&rv) *= bigtens[n_bigtens-1]; 660 e1 -= 1 << (n_bigtens-1); 661 } 662 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias; 663 word0(&rv) &= ~Exp_mask; 664 word0(&rv) |= Bias << Exp_shift1; 665 for(j = 0; e1 > 0; j++, e1 >>= 1) 666 if (e1 & 1) 667 dval(&rv) *= bigtens[j]; 668 } 669 } 670 else if (e1 < 0) { 671 e1 = -e1; 672 if ( (i = e1 & 15) !=0) 673 dval(&rv) /= tens[i]; 674 if (e1 &= ~15) { 675 e1 >>= 4; 676 while(e1 >= (1 << (n_bigtens-1))) { 677 e2 += ((word0(&rv) & Exp_mask) 678 >> Exp_shift1) - Bias; 679 word0(&rv) &= ~Exp_mask; 680 word0(&rv) |= Bias << Exp_shift1; 681 dval(&rv) *= tinytens[n_bigtens-1]; 682 e1 -= 1 << (n_bigtens-1); 683 } 684 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias; 685 word0(&rv) &= ~Exp_mask; 686 word0(&rv) |= Bias << Exp_shift1; 687 for(j = 0; e1 > 0; j++, e1 >>= 1) 688 if (e1 & 1) 689 dval(&rv) *= tinytens[j]; 690 } 691 } 692 #ifdef IBM 693 /* e2 is a correction to the (base 2) exponent of the return 694 * value, reflecting adjustments above to avoid overflow in the 695 * native arithmetic. For native IBM (base 16) arithmetic, we 696 * must multiply e2 by 4 to change from base 16 to 2. 697 */ 698 e2 <<= 2; 699 #endif 700 rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */ 701 if (rvb == NULL) 702 return (STRTOG_NoMemory); 703 rve += e2; 704 if ((j = rvbits - nbits) > 0) { 705 rshift(rvb, j); 706 rvbits = nbits; 707 rve += j; 708 } 709 bb0 = 0; /* trailing zero bits in rvb */ 710 e2 = rve + rvbits - nbits; 711 if (e2 > fpi->emax + 1) 712 goto huge; 713 rve1 = rve + rvbits - nbits; 714 if (e2 < (emin = fpi->emin)) { 715 denorm = 1; 716 j = rve - emin; 717 if (j > 0) { 718 rvb = lshift(rvb, j); 719 if (rvb == NULL) 720 return (STRTOG_NoMemory); 721 rvbits += j; 722 } 723 else if (j < 0) { 724 rvbits += j; 725 if (rvbits <= 0) { 726 if (rvbits < -1) { 727 ufl: 728 rvb->wds = 0; 729 rvb->x[0] = 0; 730 *exp = emin; 731 irv = STRTOG_Underflow | STRTOG_Inexlo; 732 goto ret; 733 } 734 rvb->x[0] = rvb->wds = rvbits = 1; 735 } 736 else 737 rshift(rvb, -j); 738 } 739 rve = rve1 = emin; 740 if (sudden_underflow && e2 + 1 < emin) 741 goto ufl; 742 } 743 744 /* Now the hard part -- adjusting rv to the correct value.*/ 745 746 /* Put digits into bd: true value = bd * 10^e */ 747 748 bd0 = s2b(s0, nd0, nd, y, dplen); 749 if (bd0 == NULL) 750 return (STRTOG_NoMemory); 751 752 for(;;) { 753 bd = Balloc(bd0->k); 754 if (bd == NULL) 755 return (STRTOG_NoMemory); 756 Bcopy(bd, bd0); 757 bb = Balloc(rvb->k); 758 if (bb == NULL) 759 return (STRTOG_NoMemory); 760 Bcopy(bb, rvb); 761 bbbits = rvbits - bb0; 762 bbe = rve + bb0; 763 bs = i2b(1); 764 if (bs == NULL) 765 return (STRTOG_NoMemory); 766 767 if (e >= 0) { 768 bb2 = bb5 = 0; 769 bd2 = bd5 = e; 770 } 771 else { 772 bb2 = bb5 = -e; 773 bd2 = bd5 = 0; 774 } 775 if (bbe >= 0) 776 bb2 += bbe; 777 else 778 bd2 -= bbe; 779 bs2 = bb2; 780 j = nbits + 1 - bbbits; 781 i = bbe + bbbits - nbits; 782 if (i < emin) /* denormal */ 783 j += i - emin; 784 bb2 += j; 785 bd2 += j; 786 i = bb2 < bd2 ? bb2 : bd2; 787 if (i > bs2) 788 i = bs2; 789 if (i > 0) { 790 bb2 -= i; 791 bd2 -= i; 792 bs2 -= i; 793 } 794 if (bb5 > 0) { 795 bs = pow5mult(bs, bb5); 796 if (bs == NULL) 797 return (STRTOG_NoMemory); 798 bb1 = mult(bs, bb); 799 if (bb1 == NULL) 800 return (STRTOG_NoMemory); 801 Bfree(bb); 802 bb = bb1; 803 } 804 bb2 -= bb0; 805 if (bb2 > 0) { 806 bb = lshift(bb, bb2); 807 if (bb == NULL) 808 return (STRTOG_NoMemory); 809 } 810 else if (bb2 < 0) 811 rshift(bb, -bb2); 812 if (bd5 > 0) { 813 bd = pow5mult(bd, bd5); 814 if (bd == NULL) 815 return (STRTOG_NoMemory); 816 } 817 if (bd2 > 0) { 818 bd = lshift(bd, bd2); 819 if (bd == NULL) 820 return (STRTOG_NoMemory); 821 } 822 if (bs2 > 0) { 823 bs = lshift(bs, bs2); 824 if (bs == NULL) 825 return (STRTOG_NoMemory); 826 } 827 asub = 1; 828 inex = STRTOG_Inexhi; 829 delta = diff(bb, bd); 830 if (delta == NULL) 831 return (STRTOG_NoMemory); 832 if (delta->wds <= 1 && !delta->x[0]) 833 break; 834 dsign = delta->sign; 835 delta->sign = finished = 0; 836 L = 0; 837 i = cmp(delta, bs); 838 if (rd && i <= 0) { 839 irv = STRTOG_Normal; 840 if ( (finished = dsign ^ (rd&1)) !=0) { 841 if (dsign != 0) { 842 irv |= STRTOG_Inexhi; 843 goto adj1; 844 } 845 irv |= STRTOG_Inexlo; 846 if (rve1 == emin) 847 goto adj1; 848 for(i = 0, j = nbits; j >= ULbits; 849 i++, j -= ULbits) { 850 if (rvb->x[i] & ALL_ON) 851 goto adj1; 852 } 853 if (j > 1 && lo0bits(rvb->x + i) < j - 1) 854 goto adj1; 855 rve = rve1 - 1; 856 rvb = set_ones(rvb, rvbits = nbits); 857 if (rvb == NULL) 858 return (STRTOG_NoMemory); 859 break; 860 } 861 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; 862 break; 863 } 864 if (i < 0) { 865 /* Error is less than half an ulp -- check for 866 * special case of mantissa a power of two. 867 */ 868 irv = dsign 869 ? STRTOG_Normal | STRTOG_Inexlo 870 : STRTOG_Normal | STRTOG_Inexhi; 871 if (dsign || bbbits > 1 || denorm || rve1 == emin) 872 break; 873 delta = lshift(delta,1); 874 if (delta == NULL) 875 return (STRTOG_NoMemory); 876 if (cmp(delta, bs) > 0) { 877 irv = STRTOG_Normal | STRTOG_Inexlo; 878 goto drop_down; 879 } 880 break; 881 } 882 if (i == 0) { 883 /* exactly half-way between */ 884 if (dsign) { 885 if (denorm && all_on(rvb, rvbits)) { 886 /*boundary case -- increment exponent*/ 887 rvb->wds = 1; 888 rvb->x[0] = 1; 889 rve = emin + nbits - (rvbits = 1); 890 irv = STRTOG_Normal | STRTOG_Inexhi; 891 denorm = 0; 892 break; 893 } 894 irv = STRTOG_Normal | STRTOG_Inexlo; 895 } 896 else if (bbbits == 1) { 897 irv = STRTOG_Normal; 898 drop_down: 899 /* boundary case -- decrement exponent */ 900 if (rve1 == emin) { 901 irv = STRTOG_Normal | STRTOG_Inexhi; 902 if (rvb->wds == 1 && rvb->x[0] == 1) 903 sudden_underflow = 1; 904 break; 905 } 906 rve -= nbits; 907 rvb = set_ones(rvb, rvbits = nbits); 908 if (rvb == NULL) 909 return (STRTOG_NoMemory); 910 break; 911 } 912 else 913 irv = STRTOG_Normal | STRTOG_Inexhi; 914 if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1)) 915 break; 916 if (dsign) { 917 rvb = increment(rvb); 918 if (rvb == NULL) 919 return (STRTOG_NoMemory); 920 j = kmask & (ULbits - (rvbits & kmask)); 921 if (hi0bits(rvb->x[rvb->wds - 1]) != j) 922 rvbits++; 923 irv = STRTOG_Normal | STRTOG_Inexhi; 924 } 925 else { 926 if (bbbits == 1) 927 goto undfl; 928 decrement(rvb); 929 irv = STRTOG_Normal | STRTOG_Inexlo; 930 } 931 break; 932 } 933 if ((dval(&adj) = ratio(delta, bs)) <= 2.) { 934 adj1: 935 inex = STRTOG_Inexlo; 936 if (dsign) { 937 asub = 0; 938 inex = STRTOG_Inexhi; 939 } 940 else if (denorm && bbbits <= 1) { 941 undfl: 942 rvb->wds = 0; 943 rve = emin; 944 irv = STRTOG_Underflow | STRTOG_Inexlo; 945 break; 946 } 947 adj0 = dval(&adj) = 1.; 948 } 949 else { 950 adj0 = dval(&adj) *= 0.5; 951 if (dsign) { 952 asub = 0; 953 inex = STRTOG_Inexlo; 954 } 955 if (dval(&adj) < 2147483647.) { 956 L = adj0; 957 adj0 -= L; 958 switch(rd) { 959 case 0: 960 if (adj0 >= .5) 961 goto inc_L; 962 break; 963 case 1: 964 if (asub && adj0 > 0.) 965 goto inc_L; 966 break; 967 case 2: 968 if (!asub && adj0 > 0.) { 969 inc_L: 970 L++; 971 inex = STRTOG_Inexact - inex; 972 } 973 } 974 dval(&adj) = L; 975 } 976 } 977 y = rve + rvbits; 978 979 /* adj *= ulp(dval(&rv)); */ 980 /* if (asub) rv -= adj; else rv += adj; */ 981 982 if (!denorm && rvbits < nbits) { 983 rvb = lshift(rvb, j = nbits - rvbits); 984 if (rvb == NULL) 985 return (STRTOG_NoMemory); 986 rve -= j; 987 rvbits = nbits; 988 } 989 ab = d2b(dval(&adj), &abe, &abits); 990 if (ab == NULL) 991 return (STRTOG_NoMemory); 992 if (abe < 0) 993 rshift(ab, -abe); 994 else if (abe > 0) { 995 ab = lshift(ab, abe); 996 if (ab == NULL) 997 return (STRTOG_NoMemory); 998 } 999 rvb0 = rvb; 1000 if (asub) { 1001 /* rv -= adj; */ 1002 j = hi0bits(rvb->x[rvb->wds-1]); 1003 rvb = diff(rvb, ab); 1004 if (rvb == NULL) 1005 return (STRTOG_NoMemory); 1006 k = rvb0->wds - 1; 1007 if (denorm) 1008 /* do nothing */; 1009 else if (rvb->wds <= k 1010 || hi0bits( rvb->x[k]) > 1011 hi0bits(rvb0->x[k])) { 1012 /* unlikely; can only have lost 1 high bit */ 1013 if (rve1 == emin) { 1014 --rvbits; 1015 denorm = 1; 1016 } 1017 else { 1018 rvb = lshift(rvb, 1); 1019 if (rvb == NULL) 1020 return (STRTOG_NoMemory); 1021 --rve; 1022 --rve1; 1023 L = finished = 0; 1024 } 1025 } 1026 } 1027 else { 1028 rvb = sum(rvb, ab); 1029 if (rvb == NULL) 1030 return (STRTOG_NoMemory); 1031 k = rvb->wds - 1; 1032 if (k >= rvb0->wds 1033 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) { 1034 if (denorm) { 1035 if (++rvbits == nbits) 1036 denorm = 0; 1037 } 1038 else { 1039 rshift(rvb, 1); 1040 rve++; 1041 rve1++; 1042 L = 0; 1043 } 1044 } 1045 } 1046 Bfree(ab); 1047 Bfree(rvb0); 1048 if (finished) 1049 break; 1050 1051 z = rve + rvbits; 1052 if (y == z && L) { 1053 /* Can we stop now? */ 1054 tol = dval(&adj) * 5e-16; /* > max rel error */ 1055 dval(&adj) = adj0 - .5; 1056 if (dval(&adj) < -tol) { 1057 if (adj0 > tol) { 1058 irv |= inex; 1059 break; 1060 } 1061 } 1062 else if (dval(&adj) > tol && adj0 < 1. - tol) { 1063 irv |= inex; 1064 break; 1065 } 1066 } 1067 bb0 = denorm ? 0 : trailz(rvb); 1068 Bfree(bb); 1069 Bfree(bd); 1070 Bfree(bs); 1071 Bfree(delta); 1072 } 1073 if (!denorm && (j = nbits - rvbits)) { 1074 if (j > 0) { 1075 rvb = lshift(rvb, j); 1076 if (rvb == NULL) 1077 return (STRTOG_NoMemory); 1078 } 1079 else 1080 rshift(rvb, -j); 1081 rve -= j; 1082 } 1083 *exp = rve; 1084 Bfree(bb); 1085 Bfree(bd); 1086 Bfree(bs); 1087 Bfree(bd0); 1088 Bfree(delta); 1089 if (rve > fpi->emax) { 1090 switch(fpi->rounding & 3) { 1091 case FPI_Round_near: 1092 goto huge; 1093 case FPI_Round_up: 1094 if (!sign) 1095 goto huge; 1096 break; 1097 case FPI_Round_down: 1098 if (sign) 1099 goto huge; 1100 } 1101 /* Round to largest representable magnitude */ 1102 Bfree(rvb); 1103 rvb = 0; 1104 irv = STRTOG_Normal | STRTOG_Inexlo; 1105 *exp = fpi->emax; 1106 b = bits; 1107 be = b + ((fpi->nbits + 31) >> 5); 1108 while(b < be) 1109 *b++ = -1; 1110 if ((j = fpi->nbits & 0x1f)) 1111 *--be >>= (32 - j); 1112 goto ret; 1113 huge: 1114 rvb->wds = 0; 1115 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 1116 #ifndef NO_ERRNO 1117 errno = ERANGE; 1118 #endif 1119 infnanexp: 1120 *exp = fpi->emax + 1; 1121 } 1122 ret: 1123 if (denorm) { 1124 if (sudden_underflow) { 1125 rvb->wds = 0; 1126 irv = STRTOG_Underflow | STRTOG_Inexlo; 1127 #ifndef NO_ERRNO 1128 errno = ERANGE; 1129 #endif 1130 } 1131 else { 1132 irv = (irv & ~STRTOG_Retmask) | 1133 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); 1134 if (irv & STRTOG_Inexact) { 1135 irv |= STRTOG_Underflow; 1136 #ifndef NO_ERRNO 1137 errno = ERANGE; 1138 #endif 1139 } 1140 } 1141 } 1142 if (se) 1143 *se = (char *)s; 1144 if (sign) 1145 irv |= STRTOG_Neg; 1146 if (rvb) { 1147 copybits(bits, nbits, rvb); 1148 Bfree(rvb); 1149 } 1150 return irv; 1151 } 1152