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