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