1 /* $NetBSD: strtodg.c,v 1.5 2006/06/02 19:46:56 mrg Exp $ */ 2 3 /**************************************************************** 4 5 The author of this software is David M. Gay. 6 7 Copyright (C) 1998-2001 by Lucent Technologies 8 All Rights Reserved 9 10 Permission to use, copy, modify, and distribute this software and 11 its documentation for any purpose and without fee is hereby 12 granted, provided that the above copyright notice appear in all 13 copies and that both that the copyright notice and this 14 permission notice and warranty disclaimer appear in supporting 15 documentation, and that the name of Lucent or any of its entities 16 not be used in advertising or publicity pertaining to 17 distribution of the software without specific, written prior 18 permission. 19 20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 27 THIS SOFTWARE. 28 29 ****************************************************************/ 30 31 /* Please send bug reports to David M. Gay (dmg at acm dot org, 32 * with " at " changed at "@" and " dot " changed to "."). */ 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++ = 0xffffffffUL; 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) CONST Bigint *b; int n; 130 #else 131 all_on(CONST Bigint *b, int n) 132 #endif 133 { 134 CONST ULong *x, *xe; 135 136 x = b->x; 137 xe = x + ((unsigned int)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 = (unsigned int)(n + ((1 << kshift) - 1)) >> kshift; 157 if (b->k < k) { 158 Bfree(b); 159 b = Balloc(k); 160 } 161 k = (unsigned int)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, expt, bits, exact, rd, irv) 178 double d; CONST FPI *fpi; Long *expt; ULong *bits; int exact, rd, *irv; 179 #else 180 (double d, CONST FPI *fpi, Long *expt, 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) { 218 if (!exact) 219 goto ret; 220 if (b->x[0] & 2) 221 break; 222 goto trunc; 223 } 224 if (b->x[(unsigned int)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[(unsigned int)k1>>kshift] & 266 (1 << (k1 & kmask)); 267 rshift(b, k); 268 *irv = STRTOG_Denormal; 269 if (carry) { 270 b = increment(b); 271 inex = STRTOG_Inexhi | STRTOG_Underflow; 272 } 273 else if (lostbits) 274 inex = STRTOG_Inexlo | STRTOG_Underflow; 275 } 276 } 277 else if (e > fpi->emax) { 278 e = fpi->emax + 1; 279 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 280 #ifndef NO_ERRNO 281 errno = ERANGE; 282 #endif 283 b->wds = inex = 0; 284 } 285 *expt = e; 286 copybits(bits, nb, b); 287 *irv |= inex; 288 rv = 1; 289 ret: 290 Bfree(b); 291 return rv; 292 } 293 294 #ifndef VAX 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 #endif /* !VAX */ 318 319 int 320 strtodg 321 #ifdef KR_headers 322 (s00, se, fpi, expt, bits) 323 CONST char *s00; char **se; CONST FPI *fpi; Long *expt; ULong *bits; 324 #else 325 (CONST char *s00, char **se, CONST FPI *fpi, Long *expt, ULong *bits) 326 #endif 327 { 328 int abe, abits, asub; 329 int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm; 330 int dsign, e, e1, e2, emin, esign, finished, i, inex, irv; 331 int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign; 332 int sudden_underflow = 0; /* pacify gcc */ 333 CONST char *s, *s0, *s1; 334 double adj, adj0, rv, tol; 335 Long L; 336 ULong y, z; 337 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; 338 339 e2 = 0; /* XXX gcc */ 340 341 irv = STRTOG_Zero; 342 denorm = sign = nz0 = nz = 0; 343 dval(rv) = 0.; 344 rvb = 0; 345 nbits = fpi->nbits; 346 for(s = s00;;s++) switch(*s) { 347 case '-': 348 sign = 1; 349 /* FALLTHROUGH */ 350 case '+': 351 if (*++s) 352 goto break2; 353 /* FALLTHROUGH */ 354 case 0: 355 sign = 0; 356 irv = STRTOG_NoNumber; 357 s = s00; 358 goto ret; 359 case '\t': 360 case '\n': 361 case '\v': 362 case '\f': 363 case '\r': 364 case ' ': 365 continue; 366 default: 367 goto break2; 368 } 369 break2: 370 if (*s == '0') { 371 #ifndef NO_HEX_FP 372 switch(s[1]) { 373 case 'x': 374 case 'X': 375 irv = gethex(&s, fpi, expt, &rvb, sign); 376 if (irv == STRTOG_NoNumber) { 377 s = s00; 378 sign = 0; 379 } 380 goto ret; 381 } 382 #endif 383 nz0 = 1; 384 while(*++s == '0') ; 385 if (!*s) 386 goto ret; 387 } 388 sudden_underflow = fpi->sudden_underflow; 389 s0 = s; 390 y = z = 0; 391 for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 392 if (nd < 9) 393 y = 10*y + c - '0'; 394 else if (nd < 16) 395 z = 10*z + c - '0'; 396 nd0 = nd; 397 #ifdef USE_LOCALE 398 if (c == *localeconv()->decimal_point) 399 #else 400 if (c == '.') 401 #endif 402 { 403 decpt = 1; 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 /* FALLTHROUGH */ 448 case '+': 449 c = *++s; 450 } 451 if (c >= '0' && c <= '9') { 452 while(c == '0') 453 c = *++s; 454 if (c > '0' && c <= '9') { 455 L = c - '0'; 456 s1 = s; 457 while((c = *++s) >= '0' && c <= '9') 458 L = 10*L + c - '0'; 459 if (s - s1 > 8 || L > 19999) 460 /* Avoid confusion from exponents 461 * so large that e might overflow. 462 */ 463 e = 19999; /* safe for 16 bit ints */ 464 else 465 e = (int)L; 466 if (esign) 467 e = -e; 468 } 469 else 470 e = 0; 471 } 472 else 473 s = s00; 474 } 475 if (!nd) { 476 if (!nz && !nz0) { 477 #ifdef INFNAN_CHECK 478 /* Check for Nan and Infinity */ 479 if (!decpt) 480 switch(c) { 481 case 'i': 482 case 'I': 483 if (match(&s,"nf")) { 484 --s; 485 if (!match(&s,"inity")) 486 ++s; 487 irv = STRTOG_Infinite; 488 goto infnanexp; 489 } 490 break; 491 case 'n': 492 case 'N': 493 if (match(&s, "an")) { 494 irv = STRTOG_NaN; 495 *expt = fpi->emax + 1; 496 #ifndef No_Hex_NaN 497 if (*s == '(') /*)*/ 498 irv = hexnan(&s, fpi, bits); 499 #endif 500 goto infnanexp; 501 } 502 } 503 #endif /* INFNAN_CHECK */ 504 irv = STRTOG_NoNumber; 505 s = s00; 506 } 507 goto ret; 508 } 509 510 irv = STRTOG_Normal; 511 e1 = e -= nf; 512 rd = 0; 513 switch(fpi->rounding & 3) { 514 case FPI_Round_up: 515 rd = 2 - sign; 516 break; 517 case FPI_Round_zero: 518 rd = 1; 519 break; 520 case FPI_Round_down: 521 rd = 1 + sign; 522 } 523 524 /* Now we have nd0 digits, starting at s0, followed by a 525 * decimal point, followed by nd-nd0 digits. The number we're 526 * after is the integer represented by those digits times 527 * 10**e */ 528 529 if (!nd0) 530 nd0 = nd; 531 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 532 dval(rv) = y; 533 if (k > 9) 534 dval(rv) = tens[k - 9] * dval(rv) + z; 535 bd0 = 0; 536 if (nbits <= P && nd <= DBL_DIG) { 537 if (!e) { 538 if (rvOK(dval(rv), fpi, expt, bits, 1, rd, &irv)) 539 goto ret; 540 } 541 else if (e > 0) { 542 if (e <= Ten_pmax) { 543 #ifdef VAX 544 goto vax_ovfl_check; 545 #else 546 i = fivesbits[e] + mantbits(dval(rv)) <= P; 547 /* rv = */ rounded_product(dval(rv), tens[e]); 548 if (rvOK(dval(rv), fpi, expt, bits, i, rd, &irv)) 549 goto ret; 550 e1 -= e; 551 goto rv_notOK; 552 #endif 553 } 554 i = DBL_DIG - nd; 555 if (e <= Ten_pmax + i) { 556 /* A fancier test would sometimes let us do 557 * this for larger i values. 558 */ 559 e2 = e - i; 560 e1 -= i; 561 dval(rv) *= tens[i]; 562 #ifdef VAX 563 /* VAX exponent range is so narrow we must 564 * worry about overflow here... 565 */ 566 vax_ovfl_check: 567 dval(adj) = dval(rv); 568 word0(adj) -= P*Exp_msk1; 569 /* adj = */ rounded_product(dval(adj), tens[e2]); 570 if ((word0(adj) & Exp_mask) 571 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) 572 goto rv_notOK; 573 word0(adj) += P*Exp_msk1; 574 dval(rv) = dval(adj); 575 #else 576 /* rv = */ rounded_product(dval(rv), tens[e2]); 577 #endif 578 if (rvOK(dval(rv), fpi, expt, bits, 0, rd, &irv)) 579 goto ret; 580 e1 -= e2; 581 } 582 } 583 #ifndef Inaccurate_Divide 584 else if (e >= -Ten_pmax) { 585 /* rv = */ rounded_quotient(dval(rv), tens[-e]); 586 if (rvOK(dval(rv), fpi, expt, bits, 0, rd, &irv)) 587 goto ret; 588 e1 -= e; 589 } 590 #endif 591 } 592 rv_notOK: 593 e1 += nd - k; 594 595 /* Get starting approximation = rv * 10**e1 */ 596 597 e2 = 0; 598 if (e1 > 0) { 599 if ( (i = e1 & 15) !=0) 600 dval(rv) *= tens[i]; 601 if (e1 &= ~15) { 602 e1 = (unsigned int)e1 >> 4; 603 while(e1 >= (1 << (n_bigtens-1))) { 604 e2 += ((word0(rv) & Exp_mask) 605 >> Exp_shift1) - Bias; 606 word0(rv) &= ~Exp_mask; 607 word0(rv) |= Bias << Exp_shift1; 608 dval(rv) *= bigtens[n_bigtens-1]; 609 e1 -= 1 << (n_bigtens-1); 610 } 611 e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; 612 word0(rv) &= ~Exp_mask; 613 word0(rv) |= Bias << Exp_shift1; 614 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1) 615 if (e1 & 1) 616 dval(rv) *= bigtens[j]; 617 } 618 } 619 else if (e1 < 0) { 620 e1 = -e1; 621 if ( (i = e1 & 15) !=0) 622 dval(rv) /= tens[i]; 623 if (e1 &= ~15) { 624 e1 = (unsigned int)e1 >> 4; 625 while(e1 >= (1 << (n_bigtens-1))) { 626 e2 += ((word0(rv) & Exp_mask) 627 >> Exp_shift1) - Bias; 628 word0(rv) &= ~Exp_mask; 629 word0(rv) |= Bias << Exp_shift1; 630 dval(rv) *= tinytens[n_bigtens-1]; 631 e1 -= 1 << (n_bigtens-1); 632 } 633 e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; 634 word0(rv) &= ~Exp_mask; 635 word0(rv) |= Bias << Exp_shift1; 636 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1) 637 if (e1 & 1) 638 dval(rv) *= tinytens[j]; 639 } 640 } 641 #ifdef IBM 642 /* e2 is a correction to the (base 2) exponent of the return 643 * value, reflecting adjustments above to avoid overflow in the 644 * native arithmetic. For native IBM (base 16) arithmetic, we 645 * must multiply e2 by 4 to change from base 16 to 2. 646 */ 647 e2 <<= 2; 648 #endif 649 rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */ 650 rve += e2; 651 if ((j = rvbits - nbits) > 0) { 652 rshift(rvb, j); 653 rvbits = nbits; 654 rve += j; 655 } 656 bb0 = 0; /* trailing zero bits in rvb */ 657 e2 = rve + rvbits - nbits; 658 if (e2 > fpi->emax + 1) 659 goto huge; 660 rve1 = rve + rvbits - nbits; 661 if (e2 < (emin = fpi->emin)) { 662 denorm = 1; 663 j = rve - emin; 664 if (j > 0) { 665 rvb = lshift(rvb, j); 666 rvbits += j; 667 } 668 else if (j < 0) { 669 rvbits += j; 670 if (rvbits <= 0) { 671 if (rvbits < -1) { 672 ufl: 673 rvb->wds = 0; 674 rvb->x[0] = 0; 675 *expt = emin; 676 irv = STRTOG_Underflow | STRTOG_Inexlo; 677 goto ret; 678 } 679 rvb->x[0] = rvb->wds = rvbits = 1; 680 } 681 else 682 rshift(rvb, -j); 683 } 684 rve = rve1 = emin; 685 if (sudden_underflow && e2 + 1 < emin) 686 goto ufl; 687 } 688 689 /* Now the hard part -- adjusting rv to the correct value.*/ 690 691 /* Put digits into bd: true value = bd * 10^e */ 692 693 bd0 = s2b(s0, nd0, nd, y); 694 695 for(;;) { 696 bd = Balloc(bd0->k); 697 Bcopy(bd, bd0); 698 bb = Balloc(rvb->k); 699 Bcopy(bb, rvb); 700 bbbits = rvbits - bb0; 701 bbe = rve + bb0; 702 bs = i2b(1); 703 704 if (e >= 0) { 705 bb2 = bb5 = 0; 706 bd2 = bd5 = e; 707 } 708 else { 709 bb2 = bb5 = -e; 710 bd2 = bd5 = 0; 711 } 712 if (bbe >= 0) 713 bb2 += bbe; 714 else 715 bd2 -= bbe; 716 bs2 = bb2; 717 j = nbits + 1 - bbbits; 718 i = bbe + bbbits - nbits; 719 if (i < emin) /* denormal */ 720 j += i - emin; 721 bb2 += j; 722 bd2 += j; 723 i = bb2 < bd2 ? bb2 : bd2; 724 if (i > bs2) 725 i = bs2; 726 if (i > 0) { 727 bb2 -= i; 728 bd2 -= i; 729 bs2 -= i; 730 } 731 if (bb5 > 0) { 732 bs = pow5mult(bs, bb5); 733 bb1 = mult(bs, bb); 734 Bfree(bb); 735 bb = bb1; 736 } 737 bb2 -= bb0; 738 if (bb2 > 0) 739 bb = lshift(bb, bb2); 740 else if (bb2 < 0) 741 rshift(bb, -bb2); 742 if (bd5 > 0) 743 bd = pow5mult(bd, bd5); 744 if (bd2 > 0) 745 bd = lshift(bd, bd2); 746 if (bs2 > 0) 747 bs = lshift(bs, bs2); 748 asub = 1; 749 inex = STRTOG_Inexhi; 750 delta = diff(bb, bd); 751 if (delta->wds <= 1 && !delta->x[0]) 752 break; 753 dsign = delta->sign; 754 delta->sign = finished = 0; 755 L = 0; 756 i = cmp(delta, bs); 757 if (rd && i <= 0) { 758 irv = STRTOG_Normal; 759 if ( (finished = dsign ^ (rd&1)) !=0) { 760 if (dsign != 0) { 761 irv |= STRTOG_Inexhi; 762 goto adj1; 763 } 764 irv |= STRTOG_Inexlo; 765 if (rve1 == emin) 766 goto adj1; 767 for(i = 0, j = nbits; j >= ULbits; 768 i++, j -= ULbits) { 769 if (rvb->x[i] & ALL_ON) 770 goto adj1; 771 } 772 if (j > 1 && lo0bits(rvb->x + i) < j - 1) 773 goto adj1; 774 rve = rve1 - 1; 775 rvb = set_ones(rvb, rvbits = nbits); 776 break; 777 } 778 irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; 779 break; 780 } 781 if (i < 0) { 782 /* Error is less than half an ulp -- check for 783 * special case of mantissa a power of two. 784 */ 785 irv = dsign 786 ? STRTOG_Normal | STRTOG_Inexlo 787 : STRTOG_Normal | STRTOG_Inexhi; 788 if (dsign || bbbits > 1 || denorm || rve1 == emin) 789 break; 790 delta = lshift(delta,1); 791 if (cmp(delta, bs) > 0) { 792 irv = STRTOG_Normal | STRTOG_Inexlo; 793 goto drop_down; 794 } 795 break; 796 } 797 if (i == 0) { 798 /* exactly half-way between */ 799 if (dsign) { 800 if (denorm && all_on(rvb, rvbits)) { 801 /*boundary case -- increment exponent*/ 802 rvb->wds = 1; 803 rvb->x[0] = 1; 804 rve = emin + nbits - (rvbits = 1); 805 irv = STRTOG_Normal | STRTOG_Inexhi; 806 denorm = 0; 807 break; 808 } 809 irv = STRTOG_Normal | STRTOG_Inexlo; 810 } 811 else if (bbbits == 1) { 812 irv = STRTOG_Normal; 813 drop_down: 814 /* boundary case -- decrement exponent */ 815 if (rve1 == emin) { 816 irv = STRTOG_Normal | STRTOG_Inexhi; 817 if (rvb->wds == 1 && rvb->x[0] == 1) 818 sudden_underflow = 1; 819 break; 820 } 821 rve -= nbits; 822 rvb = set_ones(rvb, rvbits = nbits); 823 break; 824 } 825 else 826 irv = STRTOG_Normal | STRTOG_Inexhi; 827 if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1)) 828 break; 829 if (dsign) { 830 rvb = increment(rvb); 831 if ( (j = rvbits & kmask) !=0) 832 j = ULbits - j; 833 if (hi0bits(rvb->x[(unsigned int)(rvb->wds - 1) 834 >> kshift]) 835 != j) 836 rvbits++; 837 irv = STRTOG_Normal | STRTOG_Inexhi; 838 } 839 else { 840 if (bbbits == 1) 841 goto undfl; 842 decrement(rvb); 843 irv = STRTOG_Normal | STRTOG_Inexlo; 844 } 845 break; 846 } 847 if ((dval(adj) = ratio(delta, bs)) <= 2.) { 848 adj1: 849 inex = STRTOG_Inexlo; 850 if (dsign) { 851 asub = 0; 852 inex = STRTOG_Inexhi; 853 } 854 else if (denorm && bbbits <= 1) { 855 undfl: 856 rvb->wds = 0; 857 rve = emin; 858 irv = STRTOG_Underflow | STRTOG_Inexlo; 859 break; 860 } 861 adj0 = dval(adj) = 1.; 862 } 863 else { 864 adj0 = dval(adj) *= 0.5; 865 if (dsign) { 866 asub = 0; 867 inex = STRTOG_Inexlo; 868 } 869 if (dval(adj) < 2147483647.) { 870 L = adj0; 871 adj0 -= L; 872 switch(rd) { 873 case 0: 874 if (adj0 >= .5) 875 goto inc_L; 876 break; 877 case 1: 878 if (asub && adj0 > 0.) 879 goto inc_L; 880 break; 881 case 2: 882 if (!asub && adj0 > 0.) { 883 inc_L: 884 L++; 885 inex = STRTOG_Inexact - inex; 886 } 887 } 888 dval(adj) = L; 889 } 890 } 891 y = rve + rvbits; 892 893 /* adj *= ulp(dval(rv)); */ 894 /* if (asub) rv -= adj; else rv += adj; */ 895 896 if (!denorm && rvbits < nbits) { 897 rvb = lshift(rvb, j = nbits - rvbits); 898 rve -= j; 899 rvbits = nbits; 900 } 901 ab = d2b(dval(adj), &abe, &abits); 902 if (abe < 0) 903 rshift(ab, -abe); 904 else if (abe > 0) 905 ab = lshift(ab, abe); 906 rvb0 = rvb; 907 if (asub) { 908 /* rv -= adj; */ 909 j = hi0bits(rvb->x[rvb->wds-1]); 910 rvb = diff(rvb, ab); 911 k = rvb0->wds - 1; 912 if (denorm) 913 /* do nothing */; 914 else if (rvb->wds <= k 915 || hi0bits( rvb->x[k]) > 916 hi0bits(rvb0->x[k])) { 917 /* unlikely; can only have lost 1 high bit */ 918 if (rve1 == emin) { 919 --rvbits; 920 denorm = 1; 921 } 922 else { 923 rvb = lshift(rvb, 1); 924 --rve; 925 --rve1; 926 L = finished = 0; 927 } 928 } 929 } 930 else { 931 rvb = sum(rvb, ab); 932 k = rvb->wds - 1; 933 if (k >= rvb0->wds 934 || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) { 935 if (denorm) { 936 if (++rvbits == nbits) 937 denorm = 0; 938 } 939 else { 940 rshift(rvb, 1); 941 rve++; 942 rve1++; 943 L = 0; 944 } 945 } 946 } 947 Bfree(ab); 948 Bfree(rvb0); 949 if (finished) 950 break; 951 952 z = rve + rvbits; 953 if (y == z && L) { 954 /* Can we stop now? */ 955 tol = dval(adj) * 5e-16; /* > max rel error */ 956 dval(adj) = adj0 - .5; 957 if (dval(adj) < -tol) { 958 if (adj0 > tol) { 959 irv |= inex; 960 break; 961 } 962 } 963 else if (dval(adj) > tol && adj0 < 1. - tol) { 964 irv |= inex; 965 break; 966 } 967 } 968 bb0 = denorm ? 0 : trailz(rvb); 969 Bfree(bb); 970 Bfree(bd); 971 Bfree(bs); 972 Bfree(delta); 973 } 974 if (!denorm && (j = nbits - rvbits)) { 975 if (j > 0) 976 rvb = lshift(rvb, j); 977 else 978 rshift(rvb, -j); 979 rve -= j; 980 } 981 *expt = rve; 982 Bfree(bb); 983 Bfree(bd); 984 Bfree(bs); 985 Bfree(bd0); 986 Bfree(delta); 987 if (rve > fpi->emax) { 988 huge: 989 rvb->wds = 0; 990 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; 991 #ifndef NO_ERRNO 992 errno = ERANGE; 993 #endif 994 #ifdef INFNAN_CHECK 995 infnanexp: 996 #endif 997 *expt = fpi->emax + 1; 998 } 999 ret: 1000 if (denorm) { 1001 if (sudden_underflow) { 1002 rvb->wds = 0; 1003 irv = STRTOG_Underflow | STRTOG_Inexlo; 1004 } 1005 else { 1006 irv = (irv & ~STRTOG_Retmask) | 1007 (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); 1008 if (irv & STRTOG_Inexact) 1009 irv |= STRTOG_Underflow; 1010 } 1011 } 1012 if (se) 1013 *se = __UNCONST(s); 1014 if (sign) 1015 irv |= STRTOG_Neg; 1016 if (rvb) { 1017 copybits(bits, nbits, rvb); 1018 Bfree(rvb); 1019 } 1020 return irv; 1021 } 1022