strtodg.c revision 112158
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 30 David M. Gay 31 Bell Laboratories, Room 2C-463 32 600 Mountain Avenue 33 Murray Hill, NJ 07974-0636 34 U.S.A. 35 dmg@bell-labs.com 36 */ 37 38#include "gdtoaimp.h" 39 40#ifdef USE_LOCALE 41#include "locale.h" 42#endif 43 44 static CONST int 45fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21, 46 24, 26, 28, 31, 33, 35, 38, 40, 42, 45, 47 47, 49, 52 48#ifdef VAX 49 , 54, 56 50#endif 51 }; 52 53 Bigint * 54#ifdef KR_headers 55increment(b) Bigint *b; 56#else 57increment(Bigint *b) 58#endif 59{ 60 ULong *x, *xe; 61 Bigint *b1; 62#ifdef USE_LOCALE 63 CONST char *s2; 64#endif 65#ifdef Pack_16 66 ULong carry = 1, y; 67#endif 68 69 x = b->x; 70 xe = x + b->wds; 71#ifdef Pack_32 72 do { 73 if (*x < (ULong)0xffffffffL) { 74 ++*x; 75 return b; 76 } 77 *x++ = 0; 78 } while(x < xe); 79#else 80 do { 81 y = *x + carry; 82 carry = y >> 16; 83 *x++ = y & 0xffff; 84 if (!carry) 85 return b; 86 } while(x < xe); 87 if (carry) 88#endif 89 { 90 if (b->wds >= b->maxwds) { 91 b1 = Balloc(b->k+1); 92 Bcopy(b1,b); 93 Bfree(b); 94 b = b1; 95 } 96 b->x[b->wds++] = 1; 97 } 98 return b; 99 } 100 101 int 102#ifdef KR_headers 103decrement(b) Bigint *b; 104#else 105decrement(Bigint *b) 106#endif 107{ 108 ULong *x, *xe; 109#ifdef Pack_16 110 ULong borrow = 1, y; 111#endif 112 113 x = b->x; 114 xe = x + b->wds; 115#ifdef Pack_32 116 do { 117 if (*x) { 118 --*x; 119 break; 120 } 121 *x++ = 0xffffffffL; 122 } 123 while(x < xe); 124#else 125 do { 126 y = *x - borrow; 127 borrow = (y & 0x10000) >> 16; 128 *x++ = y & 0xffff; 129 } while(borrow && x < xe); 130#endif 131 return STRTOG_Inexlo; 132 } 133 134 static int 135#ifdef KR_headers 136all_on(b, n) Bigint *b; int n; 137#else 138all_on(Bigint *b, int n) 139#endif 140{ 141 ULong *x, *xe; 142 143 x = b->x; 144 xe = x + (n >> kshift); 145 while(x < xe) 146 if ((*x++ & ALL_ON) != ALL_ON) 147 return 0; 148 if (n &= kmask) 149 return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON; 150 return 1; 151 } 152 153 Bigint * 154#ifdef KR_headers 155set_ones(b, n) Bigint *b; int n; 156#else 157set_ones(Bigint *b, int n) 158#endif 159{ 160 int k; 161 ULong *x, *xe; 162 163 k = (n + ((1 << kshift) - 1)) >> kshift; 164 if (b->k < k) { 165 Bfree(b); 166 b = Balloc(k); 167 } 168 k = n >> kshift; 169 if (n &= kmask) 170 k++; 171 b->wds = k; 172 x = b->x; 173 xe = x + k; 174 while(x < xe) 175 *x++ = ALL_ON; 176 if (n) 177 x[-1] >>= ULbits - n; 178 return b; 179 } 180 181 static int 182rvOK 183#ifdef KR_headers 184 (d, fpi, exp, bits, exact, rd, irv) 185 double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv; 186#else 187 (double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv) 188#endif 189{ 190 Bigint *b; 191 ULong carry, inex, lostbits; 192 int bdif, e, j, k, k1, nb, rv; 193 194 carry = rv = 0; 195 b = d2b(d, &e, &bdif); 196 bdif -= nb = fpi->nbits; 197 e += bdif; 198 if (bdif <= 0) { 199 if (exact) 200 goto trunc; 201 goto ret; 202 } 203 if (P == nb) { 204 if ( 205#ifndef IMPRECISE_INEXACT 206 exact && 207#endif 208 fpi->rounding == 209#ifdef RND_PRODQUOT 210 FPI_Round_near 211#else 212 Flt_Rounds 213#endif 214 ) goto trunc; 215 goto ret; 216 } 217 switch(rd) { 218 case 1: 219 goto trunc; 220 case 2: 221 break; 222 default: /* round near */ 223 k = bdif - 1; 224 if (k < 0) 225 goto trunc; 226 if (!k) { 227 if (!exact) 228 goto ret; 229 if (b->x[0] & 2) 230 break; 231 goto trunc; 232 } 233 if (b->x[k>>kshift] & ((ULong)1 << (k & kmask))) 234 break; 235 goto trunc; 236 } 237 /* "break" cases: round up 1 bit, then truncate; bdif > 0 */ 238 carry = 1; 239 trunc: 240 inex = lostbits = 0; 241 if (bdif > 0) { 242 if ( (lostbits = any_on(b, bdif)) !=0) 243 inex = STRTOG_Inexlo; 244 rshift(b, bdif); 245 if (carry) { 246 inex = STRTOG_Inexhi; 247 b = increment(b); 248 if ( (j = nb & kmask) !=0) 249 j = 32 - j; 250 if (hi0bits(b->x[b->wds - 1]) != j) { 251 if (!lostbits) 252 lostbits = b->x[0] & 1; 253 rshift(b, 1); 254 e++; 255 } 256 } 257 } 258 else if (bdif < 0) 259 b = lshift(b, -bdif); 260 if (e < fpi->emin) { 261 k = fpi->emin - e; 262 e = fpi->emin; 263 if (k > nb || fpi->sudden_underflow) { 264 b->wds = inex = 0; 265 *irv = STRTOG_Underflow | STRTOG_Inexlo; 266 } 267 else { 268 k1 = k - 1; 269 if (k1 > 0 && !lostbits) 270 lostbits = any_on(b, k1); 271 if (!lostbits && !exact) 272 goto ret; 273 lostbits |= 274 carry = b->x[k1>>kshift] & (1 << (k1 & kmask)); 275 rshift(b, k); 276 *irv = STRTOG_Denormal; 277 if (carry) { 278 b = increment(b); 279 inex = STRTOG_Inexhi | STRTOG_Underflow; 280 } 281 else if (lostbits) 282 inex = STRTOG_Inexlo | STRTOG_Underflow; 283 } 284 } 285 else if (e > fpi->emax) { 286 e = fpi->emax + 1; 287 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 288#ifndef NO_ERRNO 289 errno = ERANGE; 290#endif 291 b->wds = inex = 0; 292 } 293 *exp = e; 294 copybits(bits, nb, b); 295 *irv |= inex; 296 rv = 1; 297 ret: 298 Bfree(b); 299 return rv; 300 } 301 302 static int 303#ifdef KR_headers 304mantbits(d) double d; 305#else 306mantbits(double d) 307#endif 308{ 309 ULong L; 310#ifdef VAX 311 L = word1(d) << 16 | word1(d) >> 16; 312 if (L) 313#else 314 if ( (L = word1(d)) !=0) 315#endif 316 return P - lo0bits(&L); 317#ifdef VAX 318 L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11; 319#else 320 L = word0(d) | Exp_msk1; 321#endif 322 return P - 32 - lo0bits(&L); 323 } 324 325 int 326strtodg 327#ifdef KR_headers 328 (s00, se, fpi, exp, bits) 329 CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits; 330#else 331 (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits) 332#endif 333{ 334 int abe, abits, asub; 335 int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2; 336 int c, denorm, dsign, e, e1, e2, emin, esign, finished, i, inex, irv; 337 int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign; 338 int sudden_underflow; 339 CONST char *s, *s0, *s1; 340 double adj, adj0, rv, tol; 341 Long L; 342 ULong y, z; 343 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; 344 345 irv = STRTOG_Zero; 346 denorm = sign = nz0 = nz = 0; 347 dval(rv) = 0.; 348 rvb = 0; 349 nbits = fpi->nbits; 350 for(s = s00;;s++) switch(*s) { 351 case '-': 352 sign = 1; 353 /* no break */ 354 case '+': 355 if (*++s) 356 goto break2; 357 /* no break */ 358 case 0: 359 sign = 0; 360 irv = STRTOG_NoNumber; 361 s = s00; 362 goto ret; 363 case '\t': 364 case '\n': 365 case '\v': 366 case '\f': 367 case '\r': 368 case ' ': 369 continue; 370 default: 371 goto break2; 372 } 373 break2: 374 if (*s == '0') { 375#ifndef NO_HEX_FP 376 switch(s[1]) { 377 case 'x': 378 case 'X': 379 irv = gethex(&s, fpi, exp, &rvb, sign); 380 if (irv == STRTOG_NoNumber) { 381 s = s00; 382 sign = 0; 383 } 384 goto ret; 385 } 386#endif 387 nz0 = 1; 388 while(*++s == '0') ; 389 if (!*s) 390 goto ret; 391 } 392 sudden_underflow = fpi->sudden_underflow; 393 s0 = s; 394 y = z = 0; 395 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 396 if (nd < 9) 397 y = 10*y + c - '0'; 398 else if (nd < 16) 399 z = 10*z + c - '0'; 400 nd0 = nd; 401#ifdef USE_LOCALE 402 s1 = localeconv()->decimal_point; 403 if (c == *s1) { 404 c = '.'; 405 if (*++s1) { 406 s2 = s; 407 for(;;) { 408 if (*++s2 != *s1) { 409 c = 0; 410 break; 411 } 412 if (!*++s1) { 413 s = s2; 414 break; 415 } 416 } 417 } 418 } 419#endif 420 if (c == '.') { 421 c = *++s; 422 if (!nd) { 423 for(; c == '0'; c = *++s) 424 nz++; 425 if (c > '0' && c <= '9') { 426 s0 = s; 427 nf += nz; 428 nz = 0; 429 goto have_dig; 430 } 431 goto dig_done; 432 } 433 for(; c >= '0' && c <= '9'; c = *++s) { 434 have_dig: 435 nz++; 436 if (c -= '0') { 437 nf += nz; 438 for(i = 1; i < nz; i++) 439 if (nd++ < 9) 440 y *= 10; 441 else if (nd <= DBL_DIG + 1) 442 z *= 10; 443 if (nd++ < 9) 444 y = 10*y + c; 445 else if (nd <= DBL_DIG + 1) 446 z = 10*z + c; 447 nz = 0; 448 } 449 } 450 } 451 dig_done: 452 e = 0; 453 if (c == 'e' || c == 'E') { 454 if (!nd && !nz && !nz0) { 455 irv = STRTOG_NoNumber; 456 s = s00; 457 goto ret; 458 } 459 s00 = s; 460 esign = 0; 461 switch(c = *++s) { 462 case '-': 463 esign = 1; 464 case '+': 465 c = *++s; 466 } 467 if (c >= '0' && c <= '9') { 468 while(c == '0') 469 c = *++s; 470 if (c > '0' && c <= '9') { 471 L = c - '0'; 472 s1 = s; 473 while((c = *++s) >= '0' && c <= '9') 474 L = 10*L + c - '0'; 475 if (s - s1 > 8 || L > 19999) 476 /* Avoid confusion from exponents 477 * so large that e might overflow. 478 */ 479 e = 19999; /* safe for 16 bit ints */ 480 else 481 e = (int)L; 482 if (esign) 483 e = -e; 484 } 485 else 486 e = 0; 487 } 488 else 489 s = s00; 490 } 491 if (!nd) { 492 if (!nz && !nz0) { 493#ifdef INFNAN_CHECK 494 /* Check for Nan and Infinity */ 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 657 rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */ 658 rve += e2; 659 if ((j = rvbits - nbits) > 0) { 660 rshift(rvb, j); 661 rvbits = nbits; 662 rve += j; 663 } 664 bb0 = 0; /* trailing zero bits in rvb */ 665 e2 = rve + rvbits - nbits; 666 if (e2 > fpi->emax) { 667 rvb->wds = 0; 668 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 669#ifndef NO_ERRNO 670 errno = ERANGE; 671#endif 672 infnanexp: 673 *exp = fpi->emax + 1; 674 goto ret; 675 } 676 rve1 = rve + rvbits - nbits; 677 if (e2 < (emin = fpi->emin)) { 678 denorm = 1; 679 j = rve - emin; 680 if (j > 0) { 681 rvb = lshift(rvb, j); 682 rvbits += j; 683 } 684 else if (j < 0) { 685 rvbits += j; 686 if (rvbits <= 0) { 687 if (rvbits < -1) { 688 ufl: 689 rvb->wds = 0; 690 rvb->x[0] = 0; 691 *exp = emin; 692 irv = STRTOG_Underflow | STRTOG_Inexlo; 693 goto ret; 694 } 695 rvb->x[0] = rvb->wds = rvbits = 1; 696 } 697 else 698 rshift(rvb, -j); 699 } 700 rve = rve1 = emin; 701 if (sudden_underflow && e2 + 1 < emin) 702 goto ufl; 703 } 704 705 /* Now the hard part -- adjusting rv to the correct value.*/ 706 707 /* Put digits into bd: true value = bd * 10^e */ 708 709 bd0 = s2b(s0, nd0, nd, y); 710 711 for(;;) { 712 bd = Balloc(bd0->k); 713 Bcopy(bd, bd0); 714 bb = Balloc(rvb->k); 715 Bcopy(bb, rvb); 716 bbbits = rvbits - bb0; 717 bbe = rve + bb0; 718 bs = i2b(1); 719 720 if (e >= 0) { 721 bb2 = bb5 = 0; 722 bd2 = bd5 = e; 723 } 724 else { 725 bb2 = bb5 = -e; 726 bd2 = bd5 = 0; 727 } 728 if (bbe >= 0) 729 bb2 += bbe; 730 else 731 bd2 -= bbe; 732 bs2 = bb2; 733 j = nbits + 1 - bbbits; 734 i = bbe + bbbits - nbits; 735 if (i < emin) /* denormal */ 736 j += i - emin; 737 bb2 += j; 738 bd2 += j; 739 i = bb2 < bd2 ? bb2 : bd2; 740 if (i > bs2) 741 i = bs2; 742 if (i > 0) { 743 bb2 -= i; 744 bd2 -= i; 745 bs2 -= i; 746 } 747 if (bb5 > 0) { 748 bs = pow5mult(bs, bb5); 749 bb1 = mult(bs, bb); 750 Bfree(bb); 751 bb = bb1; 752 } 753 bb2 -= bb0; 754 if (bb2 > 0) 755 bb = lshift(bb, bb2); 756 else if (bb2 < 0) 757 rshift(bb, -bb2); 758 if (bd5 > 0) 759 bd = pow5mult(bd, bd5); 760 if (bd2 > 0) 761 bd = lshift(bd, bd2); 762 if (bs2 > 0) 763 bs = lshift(bs, bs2); 764 asub = 1; 765 inex = STRTOG_Inexhi; 766 delta = diff(bb, bd); 767 if (delta->wds <= 1 && !delta->x[0]) 768 break; 769 dsign = delta->sign; 770 delta->sign = finished = 0; 771 L = 0; 772 i = cmp(delta, bs); 773 if (rd && i <= 0) { 774 irv = STRTOG_Normal; 775 if ( (finished = dsign ^ (rd&1)) !=0) { 776 if (dsign != 0) { 777 irv |= STRTOG_Inexhi; 778 goto adj1; 779 } 780 irv |= STRTOG_Inexlo; 781 if (rve1 == emin) 782 goto adj1; 783 for(i = 0, j = nbits; j >= ULbits; 784 i++, j -= ULbits) { 785 if (rvb->x[i] & ALL_ON) 786 goto adj1; 787 } 788 if (j > 1 && lo0bits(rvb->x + i) < j - 1) 789 goto adj1; 790 rve = rve1 - 1; 791 rvb = set_ones(rvb, rvbits = nbits); 792 break; 793 } 794 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; 795 break; 796 } 797 if (i < 0) { 798 /* Error is less than half an ulp -- check for 799 * special case of mantissa a power of two. 800 */ 801 irv = dsign 802 ? STRTOG_Normal | STRTOG_Inexlo 803 : STRTOG_Normal | STRTOG_Inexhi; 804 if (dsign || bbbits > 1 || denorm || rve1 == emin) 805 break; 806 delta = lshift(delta,1); 807 if (cmp(delta, bs) > 0) { 808 irv = STRTOG_Normal | STRTOG_Inexlo; 809 goto drop_down; 810 } 811 break; 812 } 813 if (i == 0) { 814 /* exactly half-way between */ 815 if (dsign) { 816 if (denorm && all_on(rvb, rvbits)) { 817 /*boundary case -- increment exponent*/ 818 rvb->wds = 1; 819 rvb->x[0] = 1; 820 rve = emin + nbits - (rvbits = 1); 821 irv = STRTOG_Normal | STRTOG_Inexhi; 822 denorm = 0; 823 break; 824 } 825 irv = STRTOG_Normal | STRTOG_Inexlo; 826 } 827 else if (bbbits == 1) { 828 irv = STRTOG_Normal; 829 drop_down: 830 /* boundary case -- decrement exponent */ 831 if (rve1 == emin) { 832 irv = STRTOG_Normal | STRTOG_Inexhi; 833 if (rvb->wds == 1 && rvb->x[0] == 1) 834 sudden_underflow = 1; 835 break; 836 } 837 rve -= nbits; 838 rvb = set_ones(rvb, rvbits = nbits); 839 break; 840 } 841 else 842 irv = STRTOG_Normal | STRTOG_Inexhi; 843 if (bbbits < nbits && !denorm || !(rvb->x[0] & 1)) 844 break; 845 if (dsign) { 846 rvb = increment(rvb); 847 if ( (j = rvbits >> kshift) !=0) 848 j = 32 - j; 849 if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift]) 850 != j) 851 rvbits++; 852 irv = STRTOG_Normal | STRTOG_Inexhi; 853 } 854 else { 855 if (bbbits == 1) 856 goto undfl; 857 decrement(rvb); 858 irv = STRTOG_Normal | STRTOG_Inexlo; 859 } 860 break; 861 } 862 if ((dval(adj) = ratio(delta, bs)) <= 2.) { 863 adj1: 864 inex = STRTOG_Inexlo; 865 if (dsign) { 866 asub = 0; 867 inex = STRTOG_Inexhi; 868 } 869 else if (denorm && bbbits <= 1) { 870 undfl: 871 rvb->wds = 0; 872 rve = emin; 873 irv = STRTOG_Underflow | STRTOG_Inexlo; 874 break; 875 } 876 adj0 = dval(adj) = 1.; 877 } 878 else { 879 adj0 = dval(adj) *= 0.5; 880 if (dsign) { 881 asub = 0; 882 inex = STRTOG_Inexlo; 883 } 884 if (dval(adj) < 2147483647.) { 885 L = adj0; 886 adj0 -= L; 887 switch(rd) { 888 case 0: 889 if (adj0 >= .5) 890 goto inc_L; 891 break; 892 case 1: 893 if (asub && adj0 > 0.) 894 goto inc_L; 895 break; 896 case 2: 897 if (!asub && adj0 > 0.) { 898 inc_L: 899 L++; 900 inex = STRTOG_Inexact - inex; 901 } 902 } 903 dval(adj) = L; 904 } 905 } 906 y = rve + rvbits; 907 908 /* adj *= ulp(dval(rv)); */ 909 /* if (asub) rv -= adj; else rv += adj; */ 910 911 if (!denorm && rvbits < nbits) { 912 rvb = lshift(rvb, j = nbits - rvbits); 913 rve -= j; 914 rvbits = nbits; 915 } 916 ab = d2b(dval(adj), &abe, &abits); 917 if (abe < 0) 918 rshift(ab, -abe); 919 else if (abe > 0) 920 ab = lshift(ab, abe); 921 rvb0 = rvb; 922 if (asub) { 923 /* rv -= adj; */ 924 j = hi0bits(rvb->x[rvb->wds-1]); 925 rvb = diff(rvb, ab); 926 k = rvb0->wds - 1; 927 if (denorm) 928 /* do nothing */; 929 else if (rvb->wds <= k 930 || hi0bits( rvb->x[k]) > 931 hi0bits(rvb0->x[k])) { 932 /* unlikely; can only have lost 1 high bit */ 933 if (rve1 == emin) { 934 --rvbits; 935 denorm = 1; 936 } 937 else { 938 rvb = lshift(rvb, 1); 939 --rve; 940 --rve1; 941 L = finished = 0; 942 } 943 } 944 } 945 else { 946 rvb = sum(rvb, ab); 947 k = rvb->wds - 1; 948 if (k >= rvb0->wds 949 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) { 950 if (denorm) { 951 if (++rvbits == nbits) 952 denorm = 0; 953 } 954 else { 955 rshift(rvb, 1); 956 rve++; 957 rve1++; 958 L = 0; 959 } 960 } 961 } 962 Bfree(ab); 963 Bfree(rvb0); 964 if (finished) 965 break; 966 967 z = rve + rvbits; 968 if (y == z && L) { 969 /* Can we stop now? */ 970 tol = dval(adj) * 5e-16; /* > max rel error */ 971 dval(adj) = adj0 - .5; 972 if (dval(adj) < -tol) { 973 if (adj0 > tol) { 974 irv |= inex; 975 break; 976 } 977 } 978 else if (dval(adj) > tol && adj0 < 1. - tol) { 979 irv |= inex; 980 break; 981 } 982 } 983 bb0 = denorm ? 0 : trailz(rvb); 984 Bfree(bb); 985 Bfree(bd); 986 Bfree(bs); 987 Bfree(delta); 988 } 989 if (!denorm && rvbits < nbits) { 990 j = nbits - rvbits; 991 rvb = lshift(rvb, j); 992 rve -= j; 993 } 994 *exp = rve; 995 Bfree(bb); 996 Bfree(bd); 997 Bfree(bs); 998 Bfree(bd0); 999 Bfree(delta); 1000 ret: 1001 if (denorm) { 1002 if (sudden_underflow) { 1003 rvb->wds = 0; 1004 irv = STRTOG_Underflow | STRTOG_Inexlo; 1005 } 1006 else { 1007 irv = (irv & ~STRTOG_Retmask) | 1008 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); 1009 if (irv & STRTOG_Inexact) 1010 irv |= STRTOG_Underflow; 1011 } 1012 } 1013 if (se) 1014 *se = (char *)s; 1015 if (sign) 1016 irv |= STRTOG_Neg; 1017 if (rvb) { 1018 copybits(bits, nbits, rvb); 1019 Bfree(rvb); 1020 } 1021 return irv; 1022 } 1023