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