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 double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv; 176 #else 177 (double 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(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) double d; 295 #else 296 mantbits(double 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 317 #ifdef KR_headers 318 (s00, se, fpi, exp, bits) 319 CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits; 320 #else 321 (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits) 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 adj, adj0, rv, tol; 331 Long L; 332 ULong *b, *be, y, z; 333 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; 334 #ifdef USE_LOCALE 335 #ifdef NO_LOCALE_CACHE 336 char *decimalpoint = localeconv()->decimal_point; 337 #else 338 char *decimalpoint; 339 static char *decimalpoint_cache; 340 if (!(s0 = decimalpoint_cache)) { 341 size_t len; 342 s0 = localeconv()->decimal_point; 343 len = strlen(s0) + 1; 344 if ((decimalpoint_cache = (char*)malloc(len))) { 345 strlcpy(decimalpoint_cache, s0, len); 346 s0 = decimalpoint_cache; 347 } 348 } 349 decimalpoint = (char*)s0; 350 #endif 351 #endif 352 353 irv = STRTOG_Zero; 354 denorm = sign = nz0 = nz = 0; 355 dval(rv) = 0.; 356 rvb = 0; 357 nbits = fpi->nbits; 358 for(s = s00;;s++) switch(*s) { 359 case '-': 360 sign = 1; 361 /* no break */ 362 case '+': 363 if (*++s) 364 goto break2; 365 /* no break */ 366 case 0: 367 sign = 0; 368 irv = STRTOG_NoNumber; 369 s = s00; 370 goto ret; 371 case '\t': 372 case '\n': 373 case '\v': 374 case '\f': 375 case '\r': 376 case ' ': 377 continue; 378 default: 379 goto break2; 380 } 381 break2: 382 if (*s == '0') { 383 #ifndef NO_HEX_FP 384 switch(s[1]) { 385 case 'x': 386 case 'X': 387 irv = gethex(&s, fpi, exp, &rvb, sign); 388 if (irv == STRTOG_NoNumber) { 389 s = s00; 390 sign = 0; 391 } 392 goto ret; 393 } 394 #endif 395 nz0 = 1; 396 while(*++s == '0') ; 397 if (!*s) 398 goto ret; 399 } 400 sudden_underflow = fpi->sudden_underflow; 401 s0 = s; 402 y = z = 0; 403 for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 404 if (nd < 9) 405 y = 10*y + c - '0'; 406 else if (nd < 16) 407 z = 10*z + c - '0'; 408 nd0 = nd; 409 #ifdef USE_LOCALE 410 if (c == *decimalpoint) { 411 for(i = 1; decimalpoint[i]; ++i) 412 if (s[i] != decimalpoint[i]) 413 goto dig_done; 414 s += i; 415 c = *s; 416 #else 417 if (c == '.') { 418 c = *++s; 419 #endif 420 decpt = 1; 421 if (!nd) { 422 for(; c == '0'; c = *++s) 423 nz++; 424 if (c > '0' && c <= '9') { 425 s0 = s; 426 nf += nz; 427 nz = 0; 428 goto have_dig; 429 } 430 goto dig_done; 431 } 432 for(; c >= '0' && c <= '9'; c = *++s) { 433 have_dig: 434 nz++; 435 if (c -= '0') { 436 nf += nz; 437 for(i = 1; i < nz; i++) 438 if (nd++ < 9) 439 y *= 10; 440 else if (nd <= DBL_DIG + 1) 441 z *= 10; 442 if (nd++ < 9) 443 y = 10*y + c; 444 else if (nd <= DBL_DIG + 1) 445 z = 10*z + c; 446 nz = 0; 447 } 448 } 449 }/*}*/ 450 dig_done: 451 e = 0; 452 if (c == 'e' || c == 'E') { 453 if (!nd && !nz && !nz0) { 454 irv = STRTOG_NoNumber; 455 s = s00; 456 goto ret; 457 } 458 s00 = s; 459 esign = 0; 460 switch(c = *++s) { 461 case '-': 462 esign = 1; 463 case '+': 464 c = *++s; 465 } 466 if (c >= '0' && c <= '9') { 467 while(c == '0') 468 c = *++s; 469 if (c > '0' && c <= '9') { 470 L = c - '0'; 471 s1 = s; 472 while((c = *++s) >= '0' && c <= '9') 473 L = 10*L + c - '0'; 474 if (s - s1 > 8 || L > 19999) 475 /* Avoid confusion from exponents 476 * so large that e might overflow. 477 */ 478 e = 19999; /* safe for 16 bit ints */ 479 else 480 e = (int)L; 481 if (esign) 482 e = -e; 483 } 484 else 485 e = 0; 486 } 487 else 488 s = s00; 489 } 490 if (!nd) { 491 if (!nz && !nz0) { 492 #ifdef INFNAN_CHECK 493 /* Check for Nan and Infinity */ 494 if (!decpt) 495 switch(c) { 496 case 'i': 497 case 'I': 498 if (match(&s,"nf")) { 499 --s; 500 if (!match(&s,"inity")) 501 ++s; 502 irv = STRTOG_Infinite; 503 goto infnanexp; 504 } 505 break; 506 case 'n': 507 case 'N': 508 if (match(&s, "an")) { 509 irv = STRTOG_NaN; 510 *exp = fpi->emax + 1; 511 #ifndef No_Hex_NaN 512 if (*s == '(') /*)*/ 513 irv = hexnan(&s, fpi, bits); 514 #endif 515 goto infnanexp; 516 } 517 } 518 #endif /* INFNAN_CHECK */ 519 irv = STRTOG_NoNumber; 520 s = s00; 521 } 522 goto ret; 523 } 524 525 irv = STRTOG_Normal; 526 e1 = e -= nf; 527 rd = 0; 528 switch(fpi->rounding & 3) { 529 case FPI_Round_up: 530 rd = 2 - sign; 531 break; 532 case FPI_Round_zero: 533 rd = 1; 534 break; 535 case FPI_Round_down: 536 rd = 1 + sign; 537 } 538 539 /* Now we have nd0 digits, starting at s0, followed by a 540 * decimal point, followed by nd-nd0 digits. The number we're 541 * after is the integer represented by those digits times 542 * 10**e */ 543 544 if (!nd0) 545 nd0 = nd; 546 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 547 dval(rv) = y; 548 if (k > 9) 549 dval(rv) = tens[k - 9] * dval(rv) + z; 550 bd0 = 0; 551 if (nbits <= P && nd <= DBL_DIG) { 552 if (!e) { 553 if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv)) 554 goto ret; 555 } 556 else if (e > 0) { 557 if (e <= Ten_pmax) { 558 #ifdef VAX 559 goto vax_ovfl_check; 560 #else 561 i = fivesbits[e] + mantbits(dval(rv)) <= P; 562 /* rv = */ rounded_product(dval(rv), tens[e]); 563 if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv)) 564 goto ret; 565 e1 -= e; 566 goto rv_notOK; 567 #endif 568 } 569 i = DBL_DIG - nd; 570 if (e <= Ten_pmax + i) { 571 /* A fancier test would sometimes let us do 572 * this for larger i values. 573 */ 574 e2 = e - i; 575 e1 -= i; 576 dval(rv) *= tens[i]; 577 #ifdef VAX 578 /* VAX exponent range is so narrow we must 579 * worry about overflow here... 580 */ 581 vax_ovfl_check: 582 dval(adj) = dval(rv); 583 word0(adj) -= P*Exp_msk1; 584 /* adj = */ rounded_product(dval(adj), tens[e2]); 585 if ((word0(adj) & Exp_mask) 586 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 587 goto rv_notOK; 588 word0(adj) += P*Exp_msk1; 589 dval(rv) = dval(adj); 590 #else 591 /* rv = */ rounded_product(dval(rv), tens[e2]); 592 #endif 593 if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv)) 594 goto ret; 595 e1 -= e2; 596 } 597 } 598 #ifndef Inaccurate_Divide 599 else if (e >= -Ten_pmax) { 600 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 601 if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv)) 602 goto ret; 603 e1 -= e; 604 } 605 #endif 606 } 607 rv_notOK: 608 e1 += nd - k; 609 610 /* Get starting approximation = rv * 10**e1 */ 611 612 e2 = 0; 613 if (e1 > 0) { 614 if ( (i = e1 & 15) !=0) 615 dval(rv) *= tens[i]; 616 if (e1 &= ~15) { 617 e1 >>= 4; 618 while(e1 >= (1 << n_bigtens-1)) { 619 e2 += ((word0(rv) & Exp_mask) 620 >> Exp_shift1) - Bias; 621 word0(rv) &= ~Exp_mask; 622 word0(rv) |= Bias << Exp_shift1; 623 dval(rv) *= bigtens[n_bigtens-1]; 624 e1 -= 1 << n_bigtens-1; 625 } 626 e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; 627 word0(rv) &= ~Exp_mask; 628 word0(rv) |= Bias << Exp_shift1; 629 for(j = 0; e1 > 0; j++, e1 >>= 1) 630 if (e1 & 1) 631 dval(rv) *= bigtens[j]; 632 } 633 } 634 else if (e1 < 0) { 635 e1 = -e1; 636 if ( (i = e1 & 15) !=0) 637 dval(rv) /= tens[i]; 638 if (e1 &= ~15) { 639 e1 >>= 4; 640 while(e1 >= (1 << n_bigtens-1)) { 641 e2 += ((word0(rv) & Exp_mask) 642 >> Exp_shift1) - Bias; 643 word0(rv) &= ~Exp_mask; 644 word0(rv) |= Bias << Exp_shift1; 645 dval(rv) *= tinytens[n_bigtens-1]; 646 e1 -= 1 << n_bigtens-1; 647 } 648 e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; 649 word0(rv) &= ~Exp_mask; 650 word0(rv) |= Bias << Exp_shift1; 651 for(j = 0; e1 > 0; j++, e1 >>= 1) 652 if (e1 & 1) 653 dval(rv) *= tinytens[j]; 654 } 655 } 656 #ifdef IBM 657 /* e2 is a correction to the (base 2) exponent of the return 658 * value, reflecting adjustments above to avoid overflow in the 659 * native arithmetic. For native IBM (base 16) arithmetic, we 660 * must multiply e2 by 4 to change from base 16 to 2. 661 */ 662 e2 <<= 2; 663 #endif 664 rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */ 665 rve += e2; 666 if ((j = rvbits - nbits) > 0) { 667 rshift(rvb, j); 668 rvbits = nbits; 669 rve += j; 670 } 671 bb0 = 0; /* trailing zero bits in rvb */ 672 e2 = rve + rvbits - nbits; 673 if (e2 > fpi->emax + 1) 674 goto huge; 675 rve1 = rve + rvbits - nbits; 676 if (e2 < (emin = fpi->emin)) { 677 denorm = 1; 678 j = rve - emin; 679 if (j > 0) { 680 rvb = lshift(rvb, j); 681 rvbits += j; 682 } 683 else if (j < 0) { 684 rvbits += j; 685 if (rvbits <= 0) { 686 if (rvbits < -1) { 687 ufl: 688 rvb->wds = 0; 689 rvb->x[0] = 0; 690 *exp = emin; 691 irv = STRTOG_Underflow | STRTOG_Inexlo; 692 goto ret; 693 } 694 rvb->x[0] = rvb->wds = rvbits = 1; 695 } 696 else 697 rshift(rvb, -j); 698 } 699 rve = rve1 = emin; 700 if (sudden_underflow && e2 + 1 < emin) 701 goto ufl; 702 } 703 704 /* Now the hard part -- adjusting rv to the correct value.*/ 705 706 /* Put digits into bd: true value = bd * 10^e */ 707 708 bd0 = s2b(s0, nd0, nd, y); 709 710 for(;;) { 711 bd = Balloc(bd0->k); 712 Bcopy(bd, bd0); 713 bb = Balloc(rvb->k); 714 Bcopy(bb, rvb); 715 bbbits = rvbits - bb0; 716 bbe = rve + bb0; 717 bs = i2b(1); 718 719 if (e >= 0) { 720 bb2 = bb5 = 0; 721 bd2 = bd5 = e; 722 } 723 else { 724 bb2 = bb5 = -e; 725 bd2 = bd5 = 0; 726 } 727 if (bbe >= 0) 728 bb2 += bbe; 729 else 730 bd2 -= bbe; 731 bs2 = bb2; 732 j = nbits + 1 - bbbits; 733 i = bbe + bbbits - nbits; 734 if (i < emin) /* denormal */ 735 j += i - emin; 736 bb2 += j; 737 bd2 += j; 738 i = bb2 < bd2 ? bb2 : bd2; 739 if (i > bs2) 740 i = bs2; 741 if (i > 0) { 742 bb2 -= i; 743 bd2 -= i; 744 bs2 -= i; 745 } 746 if (bb5 > 0) { 747 bs = pow5mult(bs, bb5); 748 bb1 = mult(bs, bb); 749 Bfree(bb); 750 bb = bb1; 751 } 752 bb2 -= bb0; 753 if (bb2 > 0) 754 bb = lshift(bb, bb2); 755 else if (bb2 < 0) 756 rshift(bb, -bb2); 757 if (bd5 > 0) 758 bd = pow5mult(bd, bd5); 759 if (bd2 > 0) 760 bd = lshift(bd, bd2); 761 if (bs2 > 0) 762 bs = lshift(bs, bs2); 763 asub = 1; 764 inex = STRTOG_Inexhi; 765 delta = diff(bb, bd); 766 if (delta->wds <= 1 && !delta->x[0]) 767 break; 768 dsign = delta->sign; 769 delta->sign = finished = 0; 770 L = 0; 771 i = cmp(delta, bs); 772 if (rd && i <= 0) { 773 irv = STRTOG_Normal; 774 if ( (finished = dsign ^ (rd&1)) !=0) { 775 if (dsign != 0) { 776 irv |= STRTOG_Inexhi; 777 goto adj1; 778 } 779 irv |= STRTOG_Inexlo; 780 if (rve1 == emin) 781 goto adj1; 782 for(i = 0, j = nbits; j >= ULbits; 783 i++, j -= ULbits) { 784 if (rvb->x[i] & ALL_ON) 785 goto adj1; 786 } 787 if (j > 1 && lo0bits(rvb->x + i) < j - 1) 788 goto adj1; 789 rve = rve1 - 1; 790 rvb = set_ones(rvb, rvbits = nbits); 791 break; 792 } 793 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; 794 break; 795 } 796 if (i < 0) { 797 /* Error is less than half an ulp -- check for 798 * special case of mantissa a power of two. 799 */ 800 irv = dsign 801 ? STRTOG_Normal | STRTOG_Inexlo 802 : STRTOG_Normal | STRTOG_Inexhi; 803 if (dsign || bbbits > 1 || denorm || rve1 == emin) 804 break; 805 delta = lshift(delta,1); 806 if (cmp(delta, bs) > 0) { 807 irv = STRTOG_Normal | STRTOG_Inexlo; 808 goto drop_down; 809 } 810 break; 811 } 812 if (i == 0) { 813 /* exactly half-way between */ 814 if (dsign) { 815 if (denorm && all_on(rvb, rvbits)) { 816 /*boundary case -- increment exponent*/ 817 rvb->wds = 1; 818 rvb->x[0] = 1; 819 rve = emin + nbits - (rvbits = 1); 820 irv = STRTOG_Normal | STRTOG_Inexhi; 821 denorm = 0; 822 break; 823 } 824 irv = STRTOG_Normal | STRTOG_Inexlo; 825 } 826 else if (bbbits == 1) { 827 irv = STRTOG_Normal; 828 drop_down: 829 /* boundary case -- decrement exponent */ 830 if (rve1 == emin) { 831 irv = STRTOG_Normal | STRTOG_Inexhi; 832 if (rvb->wds == 1 && rvb->x[0] == 1) 833 sudden_underflow = 1; 834 break; 835 } 836 rve -= nbits; 837 rvb = set_ones(rvb, rvbits = nbits); 838 break; 839 } 840 else 841 irv = STRTOG_Normal | STRTOG_Inexhi; 842 if (bbbits < nbits && !denorm || !(rvb->x[0] & 1)) 843 break; 844 if (dsign) { 845 rvb = increment(rvb); 846 j = kmask & (ULbits - (rvbits & kmask)); 847 if (hi0bits(rvb->x[rvb->wds - 1]) != j) 848 rvbits++; 849 irv = STRTOG_Normal | STRTOG_Inexhi; 850 } 851 else { 852 if (bbbits == 1) 853 goto undfl; 854 decrement(rvb); 855 irv = STRTOG_Normal | STRTOG_Inexlo; 856 } 857 break; 858 } 859 if ((dval(adj) = ratio(delta, bs)) <= 2.) { 860 adj1: 861 inex = STRTOG_Inexlo; 862 if (dsign) { 863 asub = 0; 864 inex = STRTOG_Inexhi; 865 } 866 else if (denorm && bbbits <= 1) { 867 undfl: 868 rvb->wds = 0; 869 rve = emin; 870 irv = STRTOG_Underflow | STRTOG_Inexlo; 871 break; 872 } 873 adj0 = dval(adj) = 1.; 874 } 875 else { 876 adj0 = dval(adj) *= 0.5; 877 if (dsign) { 878 asub = 0; 879 inex = STRTOG_Inexlo; 880 } 881 if (dval(adj) < 2147483647.) { 882 L = adj0; 883 adj0 -= L; 884 switch(rd) { 885 case 0: 886 if (adj0 >= .5) 887 goto inc_L; 888 break; 889 case 1: 890 if (asub && adj0 > 0.) 891 goto inc_L; 892 break; 893 case 2: 894 if (!asub && adj0 > 0.) { 895 inc_L: 896 L++; 897 inex = STRTOG_Inexact - inex; 898 } 899 } 900 dval(adj) = L; 901 } 902 } 903 y = rve + rvbits; 904 905 /* adj *= ulp(dval(rv)); */ 906 /* if (asub) rv -= adj; else rv += adj; */ 907 908 if (!denorm && rvbits < nbits) { 909 rvb = lshift(rvb, j = nbits - rvbits); 910 rve -= j; 911 rvbits = nbits; 912 } 913 ab = d2b(dval(adj), &abe, &abits); 914 if (abe < 0) 915 rshift(ab, -abe); 916 else if (abe > 0) 917 ab = lshift(ab, abe); 918 rvb0 = rvb; 919 if (asub) { 920 /* rv -= adj; */ 921 j = hi0bits(rvb->x[rvb->wds-1]); 922 rvb = diff(rvb, ab); 923 k = rvb0->wds - 1; 924 if (denorm) 925 /* do nothing */; 926 else if (rvb->wds <= k 927 || hi0bits( rvb->x[k]) > 928 hi0bits(rvb0->x[k])) { 929 /* unlikely; can only have lost 1 high bit */ 930 if (rve1 == emin) { 931 --rvbits; 932 denorm = 1; 933 } 934 else { 935 rvb = lshift(rvb, 1); 936 --rve; 937 --rve1; 938 L = finished = 0; 939 } 940 } 941 } 942 else { 943 rvb = sum(rvb, ab); 944 k = rvb->wds - 1; 945 if (k >= rvb0->wds 946 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) { 947 if (denorm) { 948 if (++rvbits == nbits) 949 denorm = 0; 950 } 951 else { 952 rshift(rvb, 1); 953 rve++; 954 rve1++; 955 L = 0; 956 } 957 } 958 } 959 Bfree(ab); 960 Bfree(rvb0); 961 if (finished) 962 break; 963 964 z = rve + rvbits; 965 if (y == z && L) { 966 /* Can we stop now? */ 967 tol = dval(adj) * 5e-16; /* > max rel error */ 968 dval(adj) = adj0 - .5; 969 if (dval(adj) < -tol) { 970 if (adj0 > tol) { 971 irv |= inex; 972 break; 973 } 974 } 975 else if (dval(adj) > tol && adj0 < 1. - tol) { 976 irv |= inex; 977 break; 978 } 979 } 980 bb0 = denorm ? 0 : trailz(rvb); 981 Bfree(bb); 982 Bfree(bd); 983 Bfree(bs); 984 Bfree(delta); 985 } 986 if (!denorm && (j = nbits - rvbits)) { 987 if (j > 0) 988 rvb = lshift(rvb, j); 989 else 990 rshift(rvb, -j); 991 rve -= j; 992 } 993 *exp = rve; 994 Bfree(bb); 995 Bfree(bd); 996 Bfree(bs); 997 Bfree(bd0); 998 Bfree(delta); 999 if (rve > fpi->emax) { 1000 switch(fpi->rounding & 3) { 1001 case FPI_Round_near: 1002 goto huge; 1003 case FPI_Round_up: 1004 if (!sign) 1005 goto huge; 1006 break; 1007 case FPI_Round_down: 1008 if (sign) 1009 goto huge; 1010 } 1011 /* Round to largest representable magnitude */ 1012 Bfree(rvb); 1013 rvb = 0; 1014 irv = STRTOG_Normal | STRTOG_Inexlo; 1015 *exp = fpi->emax; 1016 b = bits; 1017 be = b + ((fpi->nbits + 31) >> 5); 1018 while(b < be) 1019 *b++ = -1; 1020 if ((j = fpi->nbits & 0x1f)) 1021 *--be >>= (32 - j); 1022 goto ret; 1023 huge: 1024 rvb->wds = 0; 1025 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 1026 #ifndef NO_ERRNO 1027 errno = ERANGE; 1028 #endif 1029 infnanexp: 1030 *exp = fpi->emax + 1; 1031 } 1032 ret: 1033 if (denorm) { 1034 if (sudden_underflow) { 1035 rvb->wds = 0; 1036 irv = STRTOG_Underflow | STRTOG_Inexlo; 1037 #ifndef NO_ERRNO 1038 errno = ERANGE; 1039 #endif 1040 } 1041 else { 1042 irv = (irv & ~STRTOG_Retmask) | 1043 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); 1044 if (irv & STRTOG_Inexact) { 1045 irv |= STRTOG_Underflow; 1046 #ifndef NO_ERRNO 1047 errno = ERANGE; 1048 #endif 1049 } 1050 } 1051 } 1052 if (se) 1053 *se = (char *)s; 1054 if (sign) 1055 irv |= STRTOG_Neg; 1056 if (rvb) { 1057 copybits(bits, nbits, rvb); 1058 Bfree(rvb); 1059 } 1060 return irv; 1061 } 1062