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