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