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