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