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