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