1 /* derived from /netlib/fp/dtoa.c assuming IEEE, Standard C */ 2 /* kudos to dmg@bell-labs.com, gripes to ehg@bell-labs.com */ 3 #include "lib9.h" 4 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ 5 #define FREE_DTOA_LOCK(n) /*nothing*/ 6 7 /* let's provide reasonable defaults for usual implementation of IEEE f.p. */ 8 #ifndef DBL_DIG 9 #define DBL_DIG 15 10 #endif 11 #ifndef DBL_MAX_10_EXP 12 #define DBL_MAX_10_EXP 308 13 #endif 14 #ifndef DBL_MAX_EXP 15 #define DBL_MAX_EXP 1024 16 #endif 17 #ifndef FLT_RADIX 18 #define FLT_RADIX 2 19 #endif 20 #ifndef FLT_ROUNDS 21 #define FLT_ROUNDS 1 22 #endif 23 #ifndef Storeinc 24 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 25 #endif 26 27 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000; 28 29 #ifdef __LITTLE_ENDIAN 30 #define word0(x) ((unsigned long *)&x)[1] 31 #define word1(x) ((unsigned long *)&x)[0] 32 #else 33 #define word0(x) ((unsigned long *)&x)[0] 34 #define word1(x) ((unsigned long *)&x)[1] 35 #endif 36 37 /* #define P DBL_MANT_DIG */ 38 /* Ten_pmax = floor(P*log(2)/log(5)) */ 39 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 40 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 41 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 42 43 #define Exp_shift 20 44 #define Exp_shift1 20 45 #define Exp_msk1 0x100000 46 #define Exp_msk11 0x100000 47 #define Exp_mask 0x7ff00000 48 #define P 53 49 #define Bias 1023 50 #define Emin (-1022) 51 #define Exp_1 0x3ff00000 52 #define Exp_11 0x3ff00000 53 #define Ebits 11 54 #define Frac_mask 0xfffff 55 #define Frac_mask1 0xfffff 56 #define Ten_pmax 22 57 #define Bletch 0x10 58 #define Bndry_mask 0xfffff 59 #define Bndry_mask1 0xfffff 60 #define LSB 1 61 #define Sign_bit 0x80000000 62 #define Log2P 1 63 #define Tiny0 0 64 #define Tiny1 1 65 #define Quick_max 14 66 #define Int_max 14 67 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */ 68 #define Avoid_Underflow 69 70 #define rounded_product(a,b) a *= b 71 #define rounded_quotient(a,b) a /= b 72 73 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 74 #define Big1 0xffffffff 75 76 #define Kmax 15 77 78 struct 79 Bigint { 80 struct Bigint *next; 81 int k, maxwds, sign, wds; 82 unsigned long x[1]; 83 }; 84 85 typedef struct Bigint Bigint; 86 87 static Bigint *freelist[Kmax+1]; 88 89 static Bigint * 90 Balloc(int k) 91 { 92 int x; 93 Bigint * rv; 94 95 ACQUIRE_DTOA_LOCK(0); 96 if (rv = freelist[k]) { 97 freelist[k] = rv->next; 98 } else { 99 x = 1 << k; 100 rv = (Bigint * )malloc(sizeof(Bigint) + (x - 1) * sizeof(unsigned long)); 101 if(rv == nil) 102 return nil; 103 rv->k = k; 104 rv->maxwds = x; 105 } 106 FREE_DTOA_LOCK(0); 107 rv->sign = rv->wds = 0; 108 return rv; 109 } 110 111 static void 112 Bfree(Bigint *v) 113 { 114 if (v) { 115 ACQUIRE_DTOA_LOCK(0); 116 v->next = freelist[v->k]; 117 freelist[v->k] = v; 118 FREE_DTOA_LOCK(0); 119 } 120 } 121 122 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 123 y->wds*sizeof(long) + 2*sizeof(int)) 124 125 static Bigint * 126 multadd(Bigint *b, int m, int a) /* multiply by m and add a */ 127 { 128 int i, wds; 129 unsigned long * x, y; 130 unsigned long xi, z; 131 Bigint * b1; 132 133 wds = b->wds; 134 x = b->x; 135 i = 0; 136 do { 137 xi = *x; 138 y = (xi & 0xffff) * m + a; 139 z = (xi >> 16) * m + (y >> 16); 140 a = (int)(z >> 16); 141 *x++ = (z << 16) + (y & 0xffff); 142 } while (++i < wds); 143 if (a) { 144 if (wds >= b->maxwds) { 145 b1 = Balloc(b->k + 1); 146 Bcopy(b1, b); 147 Bfree(b); 148 b = b1; 149 } 150 b->x[wds++] = a; 151 b->wds = wds; 152 } 153 return b; 154 } 155 156 static Bigint * 157 s2b(const char *s, int nd0, int nd, unsigned long y9) 158 { 159 Bigint * b; 160 int i, k; 161 long x, y; 162 163 x = (nd + 8) / 9; 164 for (k = 0, y = 1; x > y; y <<= 1, k++) 165 ; 166 b = Balloc(k); 167 b->x[0] = y9; 168 b->wds = 1; 169 170 i = 9; 171 if (9 < nd0) { 172 s += 9; 173 do 174 b = multadd(b, 10, *s++ - '0'); 175 while (++i < nd0); 176 s++; 177 } else 178 s += 10; 179 for (; i < nd; i++) 180 b = multadd(b, 10, *s++ - '0'); 181 return b; 182 } 183 184 static int 185 hi0bits(register unsigned long x) 186 { 187 register int k = 0; 188 189 if (!(x & 0xffff0000)) { 190 k = 16; 191 x <<= 16; 192 } 193 if (!(x & 0xff000000)) { 194 k += 8; 195 x <<= 8; 196 } 197 if (!(x & 0xf0000000)) { 198 k += 4; 199 x <<= 4; 200 } 201 if (!(x & 0xc0000000)) { 202 k += 2; 203 x <<= 2; 204 } 205 if (!(x & 0x80000000)) { 206 k++; 207 if (!(x & 0x40000000)) 208 return 32; 209 } 210 return k; 211 } 212 213 static int 214 lo0bits(unsigned long *y) 215 { 216 register int k; 217 register unsigned long x = *y; 218 219 if (x & 7) { 220 if (x & 1) 221 return 0; 222 if (x & 2) { 223 *y = x >> 1; 224 return 1; 225 } 226 *y = x >> 2; 227 return 2; 228 } 229 k = 0; 230 if (!(x & 0xffff)) { 231 k = 16; 232 x >>= 16; 233 } 234 if (!(x & 0xff)) { 235 k += 8; 236 x >>= 8; 237 } 238 if (!(x & 0xf)) { 239 k += 4; 240 x >>= 4; 241 } 242 if (!(x & 0x3)) { 243 k += 2; 244 x >>= 2; 245 } 246 if (!(x & 1)) { 247 k++; 248 x >>= 1; 249 if (!x & 1) 250 return 32; 251 } 252 *y = x; 253 return k; 254 } 255 256 static Bigint * 257 i2b(int i) 258 { 259 Bigint * b; 260 261 b = Balloc(1); 262 b->x[0] = i; 263 b->wds = 1; 264 return b; 265 } 266 267 static Bigint * 268 mult(Bigint *a, Bigint *b) 269 { 270 Bigint * c; 271 int k, wa, wb, wc; 272 unsigned long carry, y, z; 273 unsigned long * x, *xa, *xae, *xb, *xbe, *xc, *xc0; 274 unsigned long z2; 275 276 if (a->wds < b->wds) { 277 c = a; 278 a = b; 279 b = c; 280 } 281 k = a->k; 282 wa = a->wds; 283 wb = b->wds; 284 wc = wa + wb; 285 if (wc > a->maxwds) 286 k++; 287 c = Balloc(k); 288 for (x = c->x, xa = x + wc; x < xa; x++) 289 *x = 0; 290 xa = a->x; 291 xae = xa + wa; 292 xb = b->x; 293 xbe = xb + wb; 294 xc0 = c->x; 295 for (; xb < xbe; xb++, xc0++) { 296 if (y = *xb & 0xffff) { 297 x = xa; 298 xc = xc0; 299 carry = 0; 300 do { 301 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 302 carry = z >> 16; 303 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 304 carry = z2 >> 16; 305 Storeinc(xc, z2, z); 306 } while (x < xae); 307 *xc = carry; 308 } 309 if (y = *xb >> 16) { 310 x = xa; 311 xc = xc0; 312 carry = 0; 313 z2 = *xc; 314 do { 315 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 316 carry = z >> 16; 317 Storeinc(xc, z, z2); 318 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 319 carry = z2 >> 16; 320 } while (x < xae); 321 *xc = z2; 322 } 323 } 324 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) 325 ; 326 c->wds = wc; 327 return c; 328 } 329 330 static Bigint *p5s; 331 332 static Bigint * 333 pow5mult(Bigint *b, int k) 334 { 335 Bigint * b1, *p5, *p51; 336 int i; 337 static int p05[3] = { 338 5, 25, 125 }; 339 340 if (i = k & 3) 341 b = multadd(b, p05[i-1], 0); 342 343 if (!(k >>= 2)) 344 return b; 345 if (!(p5 = p5s)) { 346 /* first time */ 347 ACQUIRE_DTOA_LOCK(1); 348 if (!(p5 = p5s)) { 349 p5 = p5s = i2b(625); 350 p5->next = 0; 351 } 352 FREE_DTOA_LOCK(1); 353 } 354 for (; ; ) { 355 if (k & 1) { 356 b1 = mult(b, p5); 357 Bfree(b); 358 b = b1; 359 } 360 if (!(k >>= 1)) 361 break; 362 if (!(p51 = p5->next)) { 363 ACQUIRE_DTOA_LOCK(1); 364 if (!(p51 = p5->next)) { 365 p51 = p5->next = mult(p5, p5); 366 p51->next = 0; 367 } 368 FREE_DTOA_LOCK(1); 369 } 370 p5 = p51; 371 } 372 return b; 373 } 374 375 static Bigint * 376 lshift(Bigint *b, int k) 377 { 378 int i, k1, n, n1; 379 Bigint * b1; 380 unsigned long * x, *x1, *xe, z; 381 382 n = k >> 5; 383 k1 = b->k; 384 n1 = n + b->wds + 1; 385 for (i = b->maxwds; n1 > i; i <<= 1) 386 k1++; 387 b1 = Balloc(k1); 388 x1 = b1->x; 389 for (i = 0; i < n; i++) 390 *x1++ = 0; 391 x = b->x; 392 xe = x + b->wds; 393 if (k &= 0x1f) { 394 k1 = 32 - k; 395 z = 0; 396 do { 397 *x1++ = *x << k | z; 398 z = *x++ >> k1; 399 } while (x < xe); 400 if (*x1 = z) 401 ++n1; 402 } else 403 do 404 *x1++ = *x++; 405 while (x < xe); 406 b1->wds = n1 - 1; 407 Bfree(b); 408 return b1; 409 } 410 411 static int 412 cmp(Bigint *a, Bigint *b) 413 { 414 unsigned long * xa, *xa0, *xb, *xb0; 415 int i, j; 416 417 i = a->wds; 418 j = b->wds; 419 if (i -= j) 420 return i; 421 xa0 = a->x; 422 xa = xa0 + j; 423 xb0 = b->x; 424 xb = xb0 + j; 425 for (; ; ) { 426 if (*--xa != *--xb) 427 return * xa < *xb ? -1 : 1; 428 if (xa <= xa0) 429 break; 430 } 431 return 0; 432 } 433 434 static Bigint * 435 diff(Bigint *a, Bigint *b) 436 { 437 Bigint * c; 438 int i, wa, wb; 439 long borrow, y; /* We need signed shifts here. */ 440 unsigned long * xa, *xae, *xb, *xbe, *xc; 441 long z; 442 443 i = cmp(a, b); 444 if (!i) { 445 c = Balloc(0); 446 c->wds = 1; 447 c->x[0] = 0; 448 return c; 449 } 450 if (i < 0) { 451 c = a; 452 a = b; 453 b = c; 454 i = 1; 455 } else 456 i = 0; 457 c = Balloc(a->k); 458 c->sign = i; 459 wa = a->wds; 460 xa = a->x; 461 xae = xa + wa; 462 wb = b->wds; 463 xb = b->x; 464 xbe = xb + wb; 465 xc = c->x; 466 borrow = 0; 467 do { 468 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 469 borrow = y >> 16; 470 Sign_Extend(borrow, y); 471 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 472 borrow = z >> 16; 473 Sign_Extend(borrow, z); 474 Storeinc(xc, z, y); 475 } while (xb < xbe); 476 while (xa < xae) { 477 y = (*xa & 0xffff) + borrow; 478 borrow = y >> 16; 479 Sign_Extend(borrow, y); 480 z = (*xa++ >> 16) + borrow; 481 borrow = z >> 16; 482 Sign_Extend(borrow, z); 483 Storeinc(xc, z, y); 484 } 485 while (!*--xc) 486 wa--; 487 c->wds = wa; 488 return c; 489 } 490 491 static double 492 ulp(double x) 493 { 494 register long L; 495 double a; 496 497 L = (word0(x) & Exp_mask) - (P - 1) * Exp_msk1; 498 #ifndef Sudden_Underflow 499 if (L > 0) { 500 #endif 501 word0(a) = L; 502 word1(a) = 0; 503 #ifndef Sudden_Underflow 504 } else { 505 L = -L >> Exp_shift; 506 if (L < Exp_shift) { 507 word0(a) = 0x80000 >> L; 508 word1(a) = 0; 509 } else { 510 word0(a) = 0; 511 L -= Exp_shift; 512 word1(a) = L >= 31 ? 1 : 1 << 31 - L; 513 } 514 } 515 #endif 516 return a; 517 } 518 519 static double 520 b2d(Bigint *a, int *e) 521 { 522 unsigned long * xa, *xa0, w, y, z; 523 int k; 524 double d; 525 #define d0 word0(d) 526 #define d1 word1(d) 527 528 xa0 = a->x; 529 xa = xa0 + a->wds; 530 y = *--xa; 531 k = hi0bits(y); 532 *e = 32 - k; 533 if (k < Ebits) { 534 d0 = Exp_1 | y >> Ebits - k; 535 w = xa > xa0 ? *--xa : 0; 536 d1 = y << (32 - Ebits) + k | w >> Ebits - k; 537 goto ret_d; 538 } 539 z = xa > xa0 ? *--xa : 0; 540 if (k -= Ebits) { 541 d0 = Exp_1 | y << k | z >> 32 - k; 542 y = xa > xa0 ? *--xa : 0; 543 d1 = z << k | y >> 32 - k; 544 } else { 545 d0 = Exp_1 | y; 546 d1 = z; 547 } 548 ret_d: 549 #undef d0 550 #undef d1 551 return d; 552 } 553 554 static Bigint * 555 d2b(double d, int *e, int *bits) 556 { 557 Bigint * b; 558 int de, i, k; 559 unsigned long * x, y, z; 560 #define d0 word0(d) 561 #define d1 word1(d) 562 563 b = Balloc(1); 564 x = b->x; 565 566 z = d0 & Frac_mask; 567 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 568 #ifdef Sudden_Underflow 569 de = (int)(d0 >> Exp_shift); 570 z |= Exp_msk11; 571 #else 572 if (de = (int)(d0 >> Exp_shift)) 573 z |= Exp_msk1; 574 #endif 575 if (y = d1) { 576 if (k = lo0bits(&y)) { 577 x[0] = y | z << 32 - k; 578 z >>= k; 579 } else 580 x[0] = y; 581 i = b->wds = (x[1] = z) ? 2 : 1; 582 } else { 583 k = lo0bits(&z); 584 x[0] = z; 585 i = b->wds = 1; 586 k += 32; 587 } 588 #ifndef Sudden_Underflow 589 if (de) { 590 #endif 591 *e = de - Bias - (P - 1) + k; 592 *bits = P - k; 593 #ifndef Sudden_Underflow 594 } else { 595 *e = de - Bias - (P - 1) + 1 + k; 596 *bits = 32 * i - hi0bits(x[i-1]); 597 } 598 #endif 599 return b; 600 } 601 602 #undef d0 603 #undef d1 604 605 static double 606 ratio(Bigint *a, Bigint *b) 607 { 608 double da, db; 609 int k, ka, kb; 610 611 da = b2d(a, &ka); 612 db = b2d(b, &kb); 613 k = ka - kb + 32 * (a->wds - b->wds); 614 if (k > 0) 615 word0(da) += k * Exp_msk1; 616 else { 617 k = -k; 618 word0(db) += k * Exp_msk1; 619 } 620 return da / db; 621 } 622 623 static const double 624 tens[] = { 625 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 626 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 627 1e20, 1e21, 1e22 628 }; 629 630 static const double 631 bigtens[] = { 632 1e16, 1e32, 1e64, 1e128, 1e256 }; 633 634 static const double tinytens[] = { 635 1e-16, 1e-32, 1e-64, 1e-128, 636 9007199254740992.e-256 637 }; 638 639 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 640 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 641 #define Scale_Bit 0x10 642 #define n_bigtens 5 643 644 #define NAN_WORD0 0x7ff80000 645 646 #define NAN_WORD1 0 647 648 static int 649 match(const char **sp, char *t) 650 { 651 int c, d; 652 const char * s = *sp; 653 654 while (d = *t++) { 655 if ((c = *++s) >= 'A' && c <= 'Z') 656 c += 'a' - 'A'; 657 if (c != d) 658 return 0; 659 } 660 *sp = s + 1; 661 return 1; 662 } 663 664 double 665 strtod(const char *s00, char **se) 666 { 667 int scale; 668 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, 669 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 670 const char * s, *s0, *s1; 671 double aadj, aadj1, adj, rv, rv0; 672 long L; 673 unsigned long y, z; 674 Bigint * bb, *bb1, *bd, *bd0, *bs, *delta; 675 sign = nz0 = nz = 0; 676 rv = 0.; 677 for (s = s00; ; s++) 678 switch (*s) { 679 case '-': 680 sign = 1; 681 /* no break */ 682 case '+': 683 if (*++s) 684 goto break2; 685 /* no break */ 686 case 0: 687 s = s00; 688 goto ret; 689 case '\t': 690 case '\n': 691 case '\v': 692 case '\f': 693 case '\r': 694 case ' ': 695 continue; 696 default: 697 goto break2; 698 } 699 break2: 700 if (*s == '0') { 701 nz0 = 1; 702 while (*++s == '0') 703 ; 704 if (!*s) 705 goto ret; 706 } 707 s0 = s; 708 y = z = 0; 709 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 710 if (nd < 9) 711 y = 10 * y + c - '0'; 712 else if (nd < 16) 713 z = 10 * z + c - '0'; 714 nd0 = nd; 715 if (c == '.') { 716 c = *++s; 717 if (!nd) { 718 for (; c == '0'; c = *++s) 719 nz++; 720 if (c > '0' && c <= '9') { 721 s0 = s; 722 nf += nz; 723 nz = 0; 724 goto have_dig; 725 } 726 goto dig_done; 727 } 728 for (; c >= '0' && c <= '9'; c = *++s) { 729 have_dig: 730 nz++; 731 if (c -= '0') { 732 nf += nz; 733 for (i = 1; i < nz; i++) 734 if (nd++ < 9) 735 y *= 10; 736 else if (nd <= DBL_DIG + 1) 737 z *= 10; 738 if (nd++ < 9) 739 y = 10 * y + c; 740 else if (nd <= DBL_DIG + 1) 741 z = 10 * z + c; 742 nz = 0; 743 } 744 } 745 } 746 dig_done: 747 e = 0; 748 if (c == 'e' || c == 'E') { 749 if (!nd && !nz && !nz0) { 750 s = s00; 751 goto ret; 752 } 753 s00 = s; 754 esign = 0; 755 switch (c = *++s) { 756 case '-': 757 esign = 1; 758 case '+': 759 c = *++s; 760 } 761 if (c >= '0' && c <= '9') { 762 while (c == '0') 763 c = *++s; 764 if (c > '0' && c <= '9') { 765 L = c - '0'; 766 s1 = s; 767 while ((c = *++s) >= '0' && c <= '9') 768 L = 10 * L + c - '0'; 769 if (s - s1 > 8 || L > 19999) 770 /* Avoid confusion from exponents 771 * so large that e might overflow. 772 */ 773 e = 19999; /* safe for 16 bit ints */ 774 else 775 e = (int)L; 776 if (esign) 777 e = -e; 778 } else 779 e = 0; 780 } else 781 s = s00; 782 } 783 if (!nd) { 784 if (!nz && !nz0) { 785 /* Check for Nan and Infinity */ 786 switch (c) { 787 case 'i': 788 case 'I': 789 if (match(&s, "nfinity")) { 790 word0(rv) = 0x7ff00000; 791 word1(rv) = 0; 792 goto ret; 793 } 794 break; 795 case 'n': 796 case 'N': 797 if (match(&s, "an")) { 798 word0(rv) = NAN_WORD0; 799 word1(rv) = NAN_WORD1; 800 goto ret; 801 } 802 } 803 s = s00; 804 } 805 goto ret; 806 } 807 e1 = e -= nf; 808 809 /* Now we have nd0 digits, starting at s0, followed by a 810 * decimal point, followed by nd-nd0 digits. The number we're 811 * after is the integer represented by those digits times 812 * 10**e */ 813 814 if (!nd0) 815 nd0 = nd; 816 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 817 rv = y; 818 if (k > 9) 819 rv = tens[k - 9] * rv + z; 820 bd0 = 0; 821 if (nd <= DBL_DIG 822 && FLT_ROUNDS == 1 823 ) { 824 if (!e) 825 goto ret; 826 if (e > 0) { 827 if (e <= Ten_pmax) { 828 /* rv = */ rounded_product(rv, tens[e]); 829 goto ret; 830 } 831 i = DBL_DIG - nd; 832 if (e <= Ten_pmax + i) { 833 /* A fancier test would sometimes let us do 834 * this for larger i values. 835 */ 836 e -= i; 837 rv *= tens[i]; 838 /* rv = */ rounded_product(rv, tens[e]); 839 goto ret; 840 } 841 } else if (e >= -Ten_pmax) { 842 /* rv = */ rounded_quotient(rv, tens[-e]); 843 goto ret; 844 } 845 } 846 e1 += nd - k; 847 848 scale = 0; 849 850 /* Get starting approximation = rv * 10**e1 */ 851 852 if (e1 > 0) { 853 if (i = e1 & 15) 854 rv *= tens[i]; 855 if (e1 &= ~15) { 856 if (e1 > DBL_MAX_10_EXP) { 857 ovfl: 858 /* Can't trust HUGE_VAL */ 859 word0(rv) = Exp_mask; 860 word1(rv) = 0; 861 if (bd0) 862 goto retfree; 863 goto ret; 864 } 865 if (e1 >>= 4) { 866 for (j = 0; e1 > 1; j++, e1 >>= 1) 867 if (e1 & 1) 868 rv *= bigtens[j]; 869 /* The last multiplication could overflow. */ 870 word0(rv) -= P * Exp_msk1; 871 rv *= bigtens[j]; 872 if ((z = word0(rv) & Exp_mask) 873 > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) 874 goto ovfl; 875 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) { 876 /* set to largest number */ 877 /* (Can't trust DBL_MAX) */ 878 word0(rv) = Big0; 879 word1(rv) = Big1; 880 } else 881 word0(rv) += P * Exp_msk1; 882 } 883 884 } 885 } else if (e1 < 0) { 886 e1 = -e1; 887 if (i = e1 & 15) 888 rv /= tens[i]; 889 if (e1 &= ~15) { 890 e1 >>= 4; 891 if (e1 >= 1 << n_bigtens) 892 goto undfl; 893 if (e1 & Scale_Bit) 894 scale = P; 895 for (j = 0; e1 > 0; j++, e1 >>= 1) 896 if (e1 & 1) 897 rv *= tinytens[j]; 898 if (!rv) { 899 undfl: 900 rv = 0.; 901 if (bd0) 902 goto retfree; 903 goto ret; 904 } 905 } 906 } 907 908 /* Now the hard part -- adjusting rv to the correct value.*/ 909 910 /* Put digits into bd: true value = bd * 10^e */ 911 912 bd0 = s2b(s0, nd0, nd, y); 913 914 for (; ; ) { 915 bd = Balloc(bd0->k); 916 Bcopy(bd, bd0); 917 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 918 bs = i2b(1); 919 920 if (e >= 0) { 921 bb2 = bb5 = 0; 922 bd2 = bd5 = e; 923 } else { 924 bb2 = bb5 = -e; 925 bd2 = bd5 = 0; 926 } 927 if (bbe >= 0) 928 bb2 += bbe; 929 else 930 bd2 -= bbe; 931 bs2 = bb2; 932 #ifdef Sudden_Underflow 933 j = P + 1 - bbbits; 934 #else 935 i = bbe + bbbits - 1; /* logb(rv) */ 936 if (i < Emin) /* denormal */ 937 j = bbe + (P - Emin); 938 else 939 j = P + 1 - bbbits; 940 #endif 941 bb2 += j; 942 bd2 += j; 943 bd2 += scale; 944 i = bb2 < bd2 ? bb2 : bd2; 945 if (i > bs2) 946 i = bs2; 947 if (i > 0) { 948 bb2 -= i; 949 bd2 -= i; 950 bs2 -= i; 951 } 952 if (bb5 > 0) { 953 bs = pow5mult(bs, bb5); 954 bb1 = mult(bs, bb); 955 Bfree(bb); 956 bb = bb1; 957 } 958 if (bb2 > 0) 959 bb = lshift(bb, bb2); 960 if (bd5 > 0) 961 bd = pow5mult(bd, bd5); 962 if (bd2 > 0) 963 bd = lshift(bd, bd2); 964 if (bs2 > 0) 965 bs = lshift(bs, bs2); 966 delta = diff(bb, bd); 967 dsign = delta->sign; 968 delta->sign = 0; 969 i = cmp(delta, bs); 970 if (i < 0) { 971 /* Error is less than half an ulp -- check for 972 * special case of mantissa a power of two. 973 */ 974 if (dsign || word1(rv) || word0(rv) & Bndry_mask 975 || (word0(rv) & Exp_mask) <= Exp_msk1 976 ) { 977 if (!delta->x[0] && delta->wds == 1) 978 dsign = 2; 979 break; 980 } 981 delta = lshift(delta, Log2P); 982 if (cmp(delta, bs) > 0) 983 goto drop_down; 984 break; 985 } 986 if (i == 0) { 987 /* exactly half-way between */ 988 if (dsign) { 989 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 990 && word1(rv) == 0xffffffff) { 991 /*boundary case -- increment exponent*/ 992 word0(rv) = (word0(rv) & Exp_mask) 993 + Exp_msk1 994 ; 995 word1(rv) = 0; 996 dsign = 0; 997 break; 998 } 999 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { 1000 dsign = 2; 1001 drop_down: 1002 /* boundary case -- decrement exponent */ 1003 #ifdef Sudden_Underflow 1004 L = word0(rv) & Exp_mask; 1005 if (L <= Exp_msk1) 1006 goto undfl; 1007 L -= Exp_msk1; 1008 #else 1009 L = (word0(rv) & Exp_mask) - Exp_msk1; 1010 #endif 1011 word0(rv) = L | Bndry_mask1; 1012 word1(rv) = 0xffffffff; 1013 break; 1014 } 1015 if (!(word1(rv) & LSB)) 1016 break; 1017 if (dsign) 1018 rv += ulp(rv); 1019 else { 1020 rv -= ulp(rv); 1021 #ifndef Sudden_Underflow 1022 if (!rv) 1023 goto undfl; 1024 #endif 1025 } 1026 dsign = 1 - dsign; 1027 break; 1028 } 1029 if ((aadj = ratio(delta, bs)) <= 2.) { 1030 if (dsign) 1031 aadj = aadj1 = 1.; 1032 else if (word1(rv) || word0(rv) & Bndry_mask) { 1033 #ifndef Sudden_Underflow 1034 if (word1(rv) == Tiny1 && !word0(rv)) 1035 goto undfl; 1036 #endif 1037 aadj = 1.; 1038 aadj1 = -1.; 1039 } else { 1040 /* special case -- power of FLT_RADIX to be */ 1041 /* rounded down... */ 1042 1043 if (aadj < 2. / FLT_RADIX) 1044 aadj = 1. / FLT_RADIX; 1045 else 1046 aadj *= 0.5; 1047 aadj1 = -aadj; 1048 } 1049 } else { 1050 aadj *= 0.5; 1051 aadj1 = dsign ? aadj : -aadj; 1052 if (FLT_ROUNDS == 0) 1053 aadj1 += 0.5; 1054 } 1055 y = word0(rv) & Exp_mask; 1056 1057 /* Check for overflow */ 1058 1059 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) { 1060 rv0 = rv; 1061 word0(rv) -= P * Exp_msk1; 1062 adj = aadj1 * ulp(rv); 1063 rv += adj; 1064 if ((word0(rv) & Exp_mask) >= 1065 Exp_msk1 * (DBL_MAX_EXP + Bias - P)) { 1066 if (word0(rv0) == Big0 && word1(rv0) == Big1) 1067 goto ovfl; 1068 word0(rv) = Big0; 1069 word1(rv) = Big1; 1070 goto cont; 1071 } else 1072 word0(rv) += P * Exp_msk1; 1073 } else { 1074 #ifdef Sudden_Underflow 1075 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) { 1076 rv0 = rv; 1077 word0(rv) += P * Exp_msk1; 1078 adj = aadj1 * ulp(rv); 1079 rv += adj; 1080 if ((word0(rv) & Exp_mask) <= P * Exp_msk1) { 1081 if (word0(rv0) == Tiny0 1082 && word1(rv0) == Tiny1) 1083 goto undfl; 1084 word0(rv) = Tiny0; 1085 word1(rv) = Tiny1; 1086 goto cont; 1087 } else 1088 word0(rv) -= P * Exp_msk1; 1089 } else { 1090 adj = aadj1 * ulp(rv); 1091 rv += adj; 1092 } 1093 #else 1094 /* Compute adj so that the IEEE rounding rules will 1095 * correctly round rv + adj in some half-way cases. 1096 * If rv * ulp(rv) is denormalized (i.e., 1097 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 1098 * trouble from bits lost to denormalization; 1099 * example: 1.2e-307 . 1100 */ 1101 if (y <= (P - 1) * Exp_msk1 && aadj >= 1.) { 1102 aadj1 = (double)(int)(aadj + 0.5); 1103 if (!dsign) 1104 aadj1 = -aadj1; 1105 } 1106 adj = aadj1 * ulp(rv); 1107 rv += adj; 1108 #endif 1109 } 1110 z = word0(rv) & Exp_mask; 1111 if (!scale) 1112 if (y == z) { 1113 /* Can we stop now? */ 1114 L = aadj; 1115 aadj -= L; 1116 /* The tolerances below are conservative. */ 1117 if (dsign || word1(rv) || word0(rv) & Bndry_mask) { 1118 if (aadj < .4999999 || aadj > .5000001) 1119 break; 1120 } else if (aadj < .4999999 / FLT_RADIX) 1121 break; 1122 } 1123 cont: 1124 Bfree(bb); 1125 Bfree(bd); 1126 Bfree(bs); 1127 Bfree(delta); 1128 } 1129 if (scale) { 1130 if ((word0(rv) & Exp_mask) <= P * Exp_msk1 1131 && word1(rv) & 1 1132 && dsign != 2) 1133 if (dsign) 1134 rv += ulp(rv); 1135 else 1136 word1(rv) &= ~1; 1137 word0(rv0) = Exp_1 - P * Exp_msk1; 1138 word1(rv0) = 0; 1139 rv *= rv0; 1140 } 1141 retfree: 1142 Bfree(bb); 1143 Bfree(bd); 1144 Bfree(bs); 1145 Bfree(bd0); 1146 Bfree(delta); 1147 ret: 1148 if (se) 1149 *se = (char *)s; 1150 return sign ? -rv : rv; 1151 } 1152 1153 static int 1154 quorem(Bigint *b, Bigint *S) 1155 { 1156 int n; 1157 long borrow, y; 1158 unsigned long carry, q, ys; 1159 unsigned long * bx, *bxe, *sx, *sxe; 1160 long z; 1161 unsigned long si, zs; 1162 1163 n = S->wds; 1164 if (b->wds < n) 1165 return 0; 1166 sx = S->x; 1167 sxe = sx + --n; 1168 bx = b->x; 1169 bxe = bx + n; 1170 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1171 if (q) { 1172 borrow = 0; 1173 carry = 0; 1174 do { 1175 si = *sx++; 1176 ys = (si & 0xffff) * q + carry; 1177 zs = (si >> 16) * q + (ys >> 16); 1178 carry = zs >> 16; 1179 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1180 borrow = y >> 16; 1181 Sign_Extend(borrow, y); 1182 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1183 borrow = z >> 16; 1184 Sign_Extend(borrow, z); 1185 Storeinc(bx, z, y); 1186 } while (sx <= sxe); 1187 if (!*bxe) { 1188 bx = b->x; 1189 while (--bxe > bx && !*bxe) 1190 --n; 1191 b->wds = n; 1192 } 1193 } 1194 if (cmp(b, S) >= 0) { 1195 q++; 1196 borrow = 0; 1197 carry = 0; 1198 bx = b->x; 1199 sx = S->x; 1200 do { 1201 si = *sx++; 1202 ys = (si & 0xffff) + carry; 1203 zs = (si >> 16) + (ys >> 16); 1204 carry = zs >> 16; 1205 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 1206 borrow = y >> 16; 1207 Sign_Extend(borrow, y); 1208 z = (*bx >> 16) - (zs & 0xffff) + borrow; 1209 borrow = z >> 16; 1210 Sign_Extend(borrow, z); 1211 Storeinc(bx, z, y); 1212 } while (sx <= sxe); 1213 bx = b->x; 1214 bxe = bx + n; 1215 if (!*bxe) { 1216 while (--bxe > bx && !*bxe) 1217 --n; 1218 b->wds = n; 1219 } 1220 } 1221 return q; 1222 } 1223 1224 static char * 1225 rv_alloc(int i) 1226 { 1227 int j, k, *r; 1228 1229 j = sizeof(unsigned long); 1230 for (k = 0; 1231 sizeof(Bigint) - sizeof(unsigned long) - sizeof(int) + j <= i; 1232 j <<= 1) 1233 k++; 1234 r = (int * )Balloc(k); 1235 *r = k; 1236 return 1237 (char *)(r + 1); 1238 } 1239 1240 static char * 1241 nrv_alloc(char *s, char **rve, int n) 1242 { 1243 char *rv, *t; 1244 1245 t = rv = rv_alloc(n); 1246 while (*t = *s++) 1247 t++; 1248 if (rve) 1249 *rve = t; 1250 return rv; 1251 } 1252 1253 /* freedtoa(s) must be used to free values s returned by dtoa 1254 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 1255 * but for consistency with earlier versions of dtoa, it is optional 1256 * when MULTIPLE_THREADS is not defined. 1257 */ 1258 1259 void 1260 freedtoa(char *s) 1261 { 1262 Bigint * b = (Bigint * )((int *)s - 1); 1263 b->maxwds = 1 << (b->k = *(int * )b); 1264 Bfree(b); 1265 } 1266 1267 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 1268 * 1269 * Inspired by "How to Print Floating-Point Numbers Accurately" by 1270 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 1271 * 1272 * Modifications: 1273 * 1. Rather than iterating, we use a simple numeric overestimate 1274 * to determine k = floor(log10(d)). We scale relevant 1275 * quantities using O(log2(k)) rather than O(k) multiplications. 1276 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 1277 * try to generate digits strictly left to right. Instead, we 1278 * compute with fewer bits and propagate the carry if necessary 1279 * when rounding the final digit up. This is often faster. 1280 * 3. Under the assumption that input will be rounded nearest, 1281 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 1282 * That is, we allow equality in stopping tests when the 1283 * round-nearest rule will give the same floating-point value 1284 * as would satisfaction of the stopping test with strict 1285 * inequality. 1286 * 4. We remove common factors of powers of 2 from relevant 1287 * quantities. 1288 * 5. When converting floating-point integers less than 1e16, 1289 * we use floating-point arithmetic rather than resorting 1290 * to multiple-precision integers. 1291 * 6. When asked to produce fewer than 15 digits, we first try 1292 * to get by with floating-point arithmetic; we resort to 1293 * multiple-precision integer arithmetic only if we cannot 1294 * guarantee that the floating-point calculation has given 1295 * the correctly rounded result. For k requested digits and 1296 * "uniformly" distributed input, the probability is 1297 * something like 10^(k-15) that we must resort to the long 1298 * calculation. 1299 */ 1300 1301 char * 1302 dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve) 1303 { 1304 /* Arguments ndigits, decpt, sign are similar to those 1305 of ecvt and fcvt; trailing zeros are suppressed from 1306 the returned string. If not null, *rve is set to point 1307 to the end of the return value. If d is +-Infinity or NaN, 1308 then *decpt is set to 9999. 1309 1310 mode: 1311 0 ==> shortest string that yields d when read in 1312 and rounded to nearest. 1313 1 ==> like 0, but with Steele & White stopping rule; 1314 e.g. with IEEE P754 arithmetic , mode 0 gives 1315 1e23 whereas mode 1 gives 9.999999999999999e22. 1316 2 ==> max(1,ndigits) significant digits. This gives a 1317 return value similar to that of ecvt, except 1318 that trailing zeros are suppressed. 1319 3 ==> through ndigits past the decimal point. This 1320 gives a return value similar to that from fcvt, 1321 except that trailing zeros are suppressed, and 1322 ndigits can be negative. 1323 4-9 should give the same return values as 2-3, i.e., 1324 4 <= mode <= 9 ==> same return as mode 1325 2 + (mode & 1). These modes are mainly for 1326 debugging; often they run slower but sometimes 1327 faster than modes 2-3. 1328 4,5,8,9 ==> left-to-right digit generation. 1329 6-9 ==> don't try fast floating-point estimate 1330 (if applicable). 1331 1332 Values of mode other than 0-9 are treated as mode 0. 1333 1334 Sufficient space is allocated to the return value 1335 to hold the suppressed trailing zeros. 1336 */ 1337 1338 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 1339 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 1340 spec_case, try_quick; 1341 long L; 1342 #ifndef Sudden_Underflow 1343 int denorm; 1344 unsigned long x; 1345 #endif 1346 Bigint * b, *b1, *delta, *mlo, *mhi, *S; 1347 double d2, ds, eps; 1348 char *s, *s0; 1349 1350 if (word0(d) & Sign_bit) { 1351 /* set sign for everything, including 0's and NaNs */ 1352 *sign = 1; 1353 word0(d) &= ~Sign_bit; /* clear sign bit */ 1354 } else 1355 *sign = 0; 1356 1357 if ((word0(d) & Exp_mask) == Exp_mask) { 1358 /* Infinity or NaN */ 1359 *decpt = 9999; 1360 if (!word1(d) && !(word0(d) & 0xfffff)) 1361 return nrv_alloc("Infinity", rve, 8); 1362 return nrv_alloc("NaN", rve, 3); 1363 } 1364 if (!d) { 1365 *decpt = 1; 1366 return nrv_alloc("0", rve, 1); 1367 } 1368 1369 b = d2b(d, &be, &bbits); 1370 #ifdef Sudden_Underflow 1371 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1)); 1372 #else 1373 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))) { 1374 #endif 1375 word0(d2) = (word0(d) & Frac_mask1) | Exp_11; 1376 word1(d2) = word1(d); 1377 1378 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 1379 * log10(x) = log(x) / log(10) 1380 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 1381 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 1382 * 1383 * This suggests computing an approximation k to log10(d) by 1384 * 1385 * k = (i - Bias)*0.301029995663981 1386 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 1387 * 1388 * We want k to be too large rather than too small. 1389 * The error in the first-order Taylor series approximation 1390 * is in our favor, so we just round up the constant enough 1391 * to compensate for any error in the multiplication of 1392 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 1393 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 1394 * adding 1e-13 to the constant term more than suffices. 1395 * Hence we adjust the constant term to 0.1760912590558. 1396 * (We could get a more accurate k by invoking log10, 1397 * but this is probably not worthwhile.) 1398 */ 1399 1400 i -= Bias; 1401 #ifndef Sudden_Underflow 1402 denorm = 0; 1403 } else { 1404 /* d is denormalized */ 1405 1406 i = bbits + be + (Bias + (P - 1) - 1); 1407 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32 1408 : word1(d) << 32 - i; 1409 d2 = x; 1410 word0(d2) -= 31 * Exp_msk1; /* adjust exponent */ 1411 i -= (Bias + (P - 1) - 1) + 1; 1412 denorm = 1; 1413 } 1414 #endif 1415 ds = (d2 - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981; 1416 k = (int)ds; 1417 if (ds < 0. && ds != k) 1418 k--; /* want k = floor(ds) */ 1419 k_check = 1; 1420 if (k >= 0 && k <= Ten_pmax) { 1421 if (d < tens[k]) 1422 k--; 1423 k_check = 0; 1424 } 1425 j = bbits - i - 1; 1426 if (j >= 0) { 1427 b2 = 0; 1428 s2 = j; 1429 } else { 1430 b2 = -j; 1431 s2 = 0; 1432 } 1433 if (k >= 0) { 1434 b5 = 0; 1435 s5 = k; 1436 s2 += k; 1437 } else { 1438 b2 -= k; 1439 b5 = -k; 1440 s5 = 0; 1441 } 1442 if (mode < 0 || mode > 9) 1443 mode = 0; 1444 try_quick = 1; 1445 if (mode > 5) { 1446 mode -= 4; 1447 try_quick = 0; 1448 } 1449 leftright = 1; 1450 switch (mode) { 1451 case 0: 1452 case 1: 1453 ilim = ilim1 = -1; 1454 i = 18; 1455 ndigits = 0; 1456 break; 1457 case 2: 1458 leftright = 0; 1459 /* no break */ 1460 case 4: 1461 if (ndigits <= 0) 1462 ndigits = 1; 1463 ilim = ilim1 = i = ndigits; 1464 break; 1465 case 3: 1466 leftright = 0; 1467 /* no break */ 1468 case 5: 1469 i = ndigits + k + 1; 1470 ilim = i; 1471 ilim1 = i - 1; 1472 if (i <= 0) 1473 i = 1; 1474 } 1475 s = s0 = rv_alloc(i); 1476 1477 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 1478 1479 /* Try to get by with floating-point arithmetic. */ 1480 1481 i = 0; 1482 d2 = d; 1483 k0 = k; 1484 ilim0 = ilim; 1485 ieps = 2; /* conservative */ 1486 if (k > 0) { 1487 ds = tens[k&0xf]; 1488 j = k >> 4; 1489 if (j & Bletch) { 1490 /* prevent overflows */ 1491 j &= Bletch - 1; 1492 d /= bigtens[n_bigtens-1]; 1493 ieps++; 1494 } 1495 for (; j; j >>= 1, i++) 1496 if (j & 1) { 1497 ieps++; 1498 ds *= bigtens[i]; 1499 } 1500 d /= ds; 1501 } else if (j1 = -k) { 1502 d *= tens[j1 & 0xf]; 1503 for (j = j1 >> 4; j; j >>= 1, i++) 1504 if (j & 1) { 1505 ieps++; 1506 d *= bigtens[i]; 1507 } 1508 } 1509 if (k_check && d < 1. && ilim > 0) { 1510 if (ilim1 <= 0) 1511 goto fast_failed; 1512 ilim = ilim1; 1513 k--; 1514 d *= 10.; 1515 ieps++; 1516 } 1517 eps = ieps * d + 7.; 1518 word0(eps) -= (P - 1) * Exp_msk1; 1519 if (ilim == 0) { 1520 S = mhi = 0; 1521 d -= 5.; 1522 if (d > eps) 1523 goto one_digit; 1524 if (d < -eps) 1525 goto no_digits; 1526 goto fast_failed; 1527 } 1528 /* Generate ilim digits, then fix them up. */ 1529 eps *= tens[ilim-1]; 1530 for (i = 1; ; i++, d *= 10.) { 1531 L = d; 1532 d -= L; 1533 *s++ = '0' + (int)L; 1534 if (i == ilim) { 1535 if (d > 0.5 + eps) 1536 goto bump_up; 1537 else if (d < 0.5 - eps) { 1538 while (*--s == '0') 1539 ; 1540 s++; 1541 goto ret1; 1542 } 1543 break; 1544 } 1545 } 1546 fast_failed: 1547 s = s0; 1548 d = d2; 1549 k = k0; 1550 ilim = ilim0; 1551 } 1552 1553 /* Do we have a "small" integer? */ 1554 1555 if (be >= 0 && k <= Int_max) { 1556 /* Yes. */ 1557 ds = tens[k]; 1558 if (ndigits < 0 && ilim <= 0) { 1559 S = mhi = 0; 1560 if (ilim < 0 || d <= 5 * ds) 1561 goto no_digits; 1562 goto one_digit; 1563 } 1564 for (i = 1; ; i++) { 1565 L = d / ds; 1566 d -= L * ds; 1567 *s++ = '0' + (int)L; 1568 if (i == ilim) { 1569 d += d; 1570 if (d > ds || d == ds && L & 1) { 1571 bump_up: 1572 while (*--s == '9') 1573 if (s == s0) { 1574 k++; 1575 *s = '0'; 1576 break; 1577 } 1578 ++ * s++; 1579 } 1580 break; 1581 } 1582 if (!(d *= 10.)) 1583 break; 1584 } 1585 goto ret1; 1586 } 1587 1588 m2 = b2; 1589 m5 = b5; 1590 mhi = mlo = 0; 1591 if (leftright) { 1592 if (mode < 2) { 1593 i = 1594 #ifndef Sudden_Underflow 1595 denorm ? be + (Bias + (P - 1) - 1 + 1) : 1596 #endif 1597 1 + P - bbits; 1598 } else { 1599 j = ilim - 1; 1600 if (m5 >= j) 1601 m5 -= j; 1602 else { 1603 s5 += j -= m5; 1604 b5 += j; 1605 m5 = 0; 1606 } 1607 if ((i = ilim) < 0) { 1608 m2 -= i; 1609 i = 0; 1610 } 1611 } 1612 b2 += i; 1613 s2 += i; 1614 mhi = i2b(1); 1615 } 1616 if (m2 > 0 && s2 > 0) { 1617 i = m2 < s2 ? m2 : s2; 1618 b2 -= i; 1619 m2 -= i; 1620 s2 -= i; 1621 } 1622 if (b5 > 0) { 1623 if (leftright) { 1624 if (m5 > 0) { 1625 mhi = pow5mult(mhi, m5); 1626 b1 = mult(mhi, b); 1627 Bfree(b); 1628 b = b1; 1629 } 1630 if (j = b5 - m5) 1631 b = pow5mult(b, j); 1632 } else 1633 b = pow5mult(b, b5); 1634 } 1635 S = i2b(1); 1636 if (s5 > 0) 1637 S = pow5mult(S, s5); 1638 1639 /* Check for special case that d is a normalized power of 2. */ 1640 1641 spec_case = 0; 1642 if (mode < 2) { 1643 if (!word1(d) && !(word0(d) & Bndry_mask) 1644 #ifndef Sudden_Underflow 1645 && word0(d) & Exp_mask 1646 #endif 1647 ) { 1648 /* The special case */ 1649 b2 += Log2P; 1650 s2 += Log2P; 1651 spec_case = 1; 1652 } 1653 } 1654 1655 /* Arrange for convenient computation of quotients: 1656 * shift left if necessary so divisor has 4 leading 0 bits. 1657 * 1658 * Perhaps we should just compute leading 28 bits of S once 1659 * and for all and pass them and a shift to quorem, so it 1660 * can do shifts and ors to compute the numerator for q. 1661 */ 1662 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) 1663 i = 32 - i; 1664 if (i > 4) { 1665 i -= 4; 1666 b2 += i; 1667 m2 += i; 1668 s2 += i; 1669 } else if (i < 4) { 1670 i += 28; 1671 b2 += i; 1672 m2 += i; 1673 s2 += i; 1674 } 1675 if (b2 > 0) 1676 b = lshift(b, b2); 1677 if (s2 > 0) 1678 S = lshift(S, s2); 1679 if (k_check) { 1680 if (cmp(b, S) < 0) { 1681 k--; 1682 b = multadd(b, 10, 0); /* we botched the k estimate */ 1683 if (leftright) 1684 mhi = multadd(mhi, 10, 0); 1685 ilim = ilim1; 1686 } 1687 } 1688 if (ilim <= 0 && mode > 2) { 1689 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) { 1690 /* no digits, fcvt style */ 1691 no_digits: 1692 k = -1 - ndigits; 1693 goto ret; 1694 } 1695 one_digit: 1696 *s++ = '1'; 1697 k++; 1698 goto ret; 1699 } 1700 if (leftright) { 1701 if (m2 > 0) 1702 mhi = lshift(mhi, m2); 1703 1704 /* Compute mlo -- check for special case 1705 * that d is a normalized power of 2. 1706 */ 1707 1708 mlo = mhi; 1709 if (spec_case) { 1710 mhi = Balloc(mhi->k); 1711 Bcopy(mhi, mlo); 1712 mhi = lshift(mhi, Log2P); 1713 } 1714 1715 for (i = 1; ; i++) { 1716 dig = quorem(b, S) + '0'; 1717 /* Do we yet have the shortest decimal string 1718 * that will round to d? 1719 */ 1720 j = cmp(b, mlo); 1721 delta = diff(S, mhi); 1722 j1 = delta->sign ? 1 : cmp(b, delta); 1723 Bfree(delta); 1724 if (j1 == 0 && !mode && !(word1(d) & 1)) { 1725 if (dig == '9') 1726 goto round_9_up; 1727 if (j > 0) 1728 dig++; 1729 *s++ = dig; 1730 goto ret; 1731 } 1732 if (j < 0 || j == 0 && !mode 1733 && !(word1(d) & 1) 1734 ) { 1735 if (j1 > 0) { 1736 b = lshift(b, 1); 1737 j1 = cmp(b, S); 1738 if ((j1 > 0 || j1 == 0 && dig & 1) 1739 && dig++ == '9') 1740 goto round_9_up; 1741 } 1742 *s++ = dig; 1743 goto ret; 1744 } 1745 if (j1 > 0) { 1746 if (dig == '9') { /* possible if i == 1 */ 1747 round_9_up: 1748 *s++ = '9'; 1749 goto roundoff; 1750 } 1751 *s++ = dig + 1; 1752 goto ret; 1753 } 1754 *s++ = dig; 1755 if (i == ilim) 1756 break; 1757 b = multadd(b, 10, 0); 1758 if (mlo == mhi) 1759 mlo = mhi = multadd(mhi, 10, 0); 1760 else { 1761 mlo = multadd(mlo, 10, 0); 1762 mhi = multadd(mhi, 10, 0); 1763 } 1764 } 1765 } else 1766 for (i = 1; ; i++) { 1767 *s++ = dig = quorem(b, S) + '0'; 1768 if (i >= ilim) 1769 break; 1770 b = multadd(b, 10, 0); 1771 } 1772 1773 /* Round off last digit */ 1774 1775 b = lshift(b, 1); 1776 j = cmp(b, S); 1777 if (j > 0 || j == 0 && dig & 1) { 1778 roundoff: 1779 while (*--s == '9') 1780 if (s == s0) { 1781 k++; 1782 *s++ = '1'; 1783 goto ret; 1784 } 1785 ++ * s++; 1786 } else { 1787 while (*--s == '0') 1788 ; 1789 s++; 1790 } 1791 ret: 1792 Bfree(S); 1793 if (mhi) { 1794 if (mlo && mlo != mhi) 1795 Bfree(mlo); 1796 Bfree(mhi); 1797 } 1798 ret1: 1799 Bfree(b); 1800 *s = 0; 1801 *decpt = k + 1; 1802 if (rve) 1803 *rve = s; 1804 return s0; 1805 } 1806 1807