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 4 /* Let x be the exact mathematical number defined by a decimal 5 * string s. Then atof(s) is the round-nearest-even IEEE 6 * floating point value. 7 * Let y be an IEEE floating point value and let s be the string 8 * printed as %.17g. Then atof(s) is exactly y. 9 */ 10 #include <u.h> 11 #include <libc.h> 12 13 static Lock _dtoalk[2]; 14 #define ACQUIRE_DTOA_LOCK(n) lock(&_dtoalk[n]) 15 #define FREE_DTOA_LOCK(n) unlock(&_dtoalk[n]) 16 17 #define PRIVATE_mem ((2000+sizeof(double)-1)/sizeof(double)) 18 static double private_mem[PRIVATE_mem], *pmem_next = private_mem; 19 20 #define FLT_ROUNDS 1 21 #define DBL_DIG 15 22 #define DBL_MAX_10_EXP 308 23 #define DBL_MAX_EXP 1024 24 #define FLT_RADIX 2 25 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 26 27 /* Ten_pmax = floor(P*log(2)/log(5)) */ 28 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 29 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 30 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 31 32 #define Exp_shift 20 33 #define Exp_shift1 20 34 #define Exp_msk1 0x100000 35 #define Exp_msk11 0x100000 36 #define Exp_mask 0x7ff00000 37 #define P 53 38 #define Bias 1023 39 #define Emin (-1022) 40 #define Exp_1 0x3ff00000 41 #define Exp_11 0x3ff00000 42 #define Ebits 11 43 #define Frac_mask 0xfffff 44 #define Frac_mask1 0xfffff 45 #define Ten_pmax 22 46 #define Bletch 0x10 47 #define Bndry_mask 0xfffff 48 #define Bndry_mask1 0xfffff 49 #define LSB 1 50 #define Sign_bit 0x80000000 51 #define Log2P 1 52 #define Tiny0 0 53 #define Tiny1 1 54 #define Quick_max 14 55 #define Int_max 14 56 #define Avoid_Underflow 57 58 #define rounded_product(a,b) a *= b 59 #define rounded_quotient(a,b) a /= b 60 61 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 62 #define Big1 0xffffffff 63 64 #define FFFFFFFF 0xffffffffUL 65 66 #define Kmax 15 67 68 typedef struct Bigint Bigint; 69 typedef struct Ulongs Ulongs; 70 71 struct Bigint { 72 Bigint *next; 73 int k, maxwds, sign, wds; 74 unsigned x[1]; 75 }; 76 77 struct Ulongs { 78 ulong hi; 79 ulong lo; 80 }; 81 82 static Bigint *freelist[Kmax+1]; 83 84 Ulongs 85 double2ulongs(double d) 86 { 87 FPdbleword dw; 88 Ulongs uls; 89 90 dw.x = d; 91 uls.hi = dw.hi; 92 uls.lo = dw.lo; 93 return uls; 94 } 95 96 double 97 ulongs2double(Ulongs uls) 98 { 99 FPdbleword dw; 100 101 dw.hi = uls.hi; 102 dw.lo = uls.lo; 103 return dw.x; 104 } 105 106 static Bigint * 107 Balloc(int k) 108 { 109 int x; 110 Bigint * rv; 111 unsigned int len; 112 113 ACQUIRE_DTOA_LOCK(0); 114 if (rv = freelist[k]) { 115 freelist[k] = rv->next; 116 } else { 117 x = 1 << k; 118 len = (sizeof(Bigint) + (x - 1) * sizeof(unsigned int) + sizeof(double) -1) 119 / sizeof(double); 120 if (pmem_next - private_mem + len <= PRIVATE_mem) { 121 rv = (Bigint * )pmem_next; 122 pmem_next += len; 123 } else 124 rv = (Bigint * )malloc(len * sizeof(double)); 125 rv->k = k; 126 rv->maxwds = x; 127 } 128 FREE_DTOA_LOCK(0); 129 rv->sign = rv->wds = 0; 130 return rv; 131 } 132 133 static void 134 Bfree(Bigint *v) 135 { 136 if (v) { 137 ACQUIRE_DTOA_LOCK(0); 138 v->next = freelist[v->k]; 139 freelist[v->k] = v; 140 FREE_DTOA_LOCK(0); 141 } 142 } 143 144 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 145 y->wds*sizeof(int) + 2*sizeof(int)) 146 147 static Bigint * 148 multadd(Bigint *b, int m, int a) /* multiply by m and add a */ 149 { 150 int i, wds; 151 unsigned int carry, *x, y; 152 unsigned int xi, z; 153 Bigint * b1; 154 155 wds = b->wds; 156 x = b->x; 157 i = 0; 158 carry = a; 159 do { 160 xi = *x; 161 y = (xi & 0xffff) * m + carry; 162 z = (xi >> 16) * m + (y >> 16); 163 carry = z >> 16; 164 *x++ = (z << 16) + (y & 0xffff); 165 } while (++i < wds); 166 if (carry) { 167 if (wds >= b->maxwds) { 168 b1 = Balloc(b->k + 1); 169 Bcopy(b1, b); 170 Bfree(b); 171 b = b1; 172 } 173 b->x[wds++] = carry; 174 b->wds = wds; 175 } 176 return b; 177 } 178 179 static Bigint * 180 s2b(const char *s, int nd0, int nd, unsigned int y9) 181 { 182 Bigint * b; 183 int i, k; 184 int x, y; 185 186 x = (nd + 8) / 9; 187 for (k = 0, y = 1; x > y; y <<= 1, k++) 188 ; 189 b = Balloc(k); 190 b->x[0] = y9; 191 b->wds = 1; 192 193 i = 9; 194 if (9 < nd0) { 195 s += 9; 196 do 197 b = multadd(b, 10, *s++ - '0'); 198 while (++i < nd0); 199 s++; 200 } else 201 s += 10; 202 for (; i < nd; i++) 203 b = multadd(b, 10, *s++ - '0'); 204 return b; 205 } 206 207 static int 208 hi0bits(register unsigned int x) 209 { 210 register int k = 0; 211 212 if (!(x & 0xffff0000)) { 213 k = 16; 214 x <<= 16; 215 } 216 if (!(x & 0xff000000)) { 217 k += 8; 218 x <<= 8; 219 } 220 if (!(x & 0xf0000000)) { 221 k += 4; 222 x <<= 4; 223 } 224 if (!(x & 0xc0000000)) { 225 k += 2; 226 x <<= 2; 227 } 228 if (!(x & 0x80000000)) { 229 k++; 230 if (!(x & 0x40000000)) 231 return 32; 232 } 233 return k; 234 } 235 236 static int 237 lo0bits(unsigned int *y) 238 { 239 register int k; 240 register unsigned int x = *y; 241 242 if (x & 7) { 243 if (x & 1) 244 return 0; 245 if (x & 2) { 246 *y = x >> 1; 247 return 1; 248 } 249 *y = x >> 2; 250 return 2; 251 } 252 k = 0; 253 if (!(x & 0xffff)) { 254 k = 16; 255 x >>= 16; 256 } 257 if (!(x & 0xff)) { 258 k += 8; 259 x >>= 8; 260 } 261 if (!(x & 0xf)) { 262 k += 4; 263 x >>= 4; 264 } 265 if (!(x & 0x3)) { 266 k += 2; 267 x >>= 2; 268 } 269 if (!(x & 1)) { 270 k++; 271 x >>= 1; 272 if (!x & 1) 273 return 32; 274 } 275 *y = x; 276 return k; 277 } 278 279 static Bigint * 280 i2b(int i) 281 { 282 Bigint * b; 283 284 b = Balloc(1); 285 b->x[0] = i; 286 b->wds = 1; 287 return b; 288 } 289 290 static Bigint * 291 mult(Bigint *a, Bigint *b) 292 { 293 Bigint * c; 294 int k, wa, wb, wc; 295 unsigned int * x, *xa, *xae, *xb, *xbe, *xc, *xc0; 296 unsigned int y; 297 unsigned int carry, z; 298 unsigned int z2; 299 300 if (a->wds < b->wds) { 301 c = a; 302 a = b; 303 b = c; 304 } 305 k = a->k; 306 wa = a->wds; 307 wb = b->wds; 308 wc = wa + wb; 309 if (wc > a->maxwds) 310 k++; 311 c = Balloc(k); 312 for (x = c->x, xa = x + wc; x < xa; x++) 313 *x = 0; 314 xa = a->x; 315 xae = xa + wa; 316 xb = b->x; 317 xbe = xb + wb; 318 xc0 = c->x; 319 for (; xb < xbe; xb++, xc0++) { 320 if (y = *xb & 0xffff) { 321 x = xa; 322 xc = xc0; 323 carry = 0; 324 do { 325 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 326 carry = z >> 16; 327 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 328 carry = z2 >> 16; 329 Storeinc(xc, z2, z); 330 } while (x < xae); 331 *xc = carry; 332 } 333 if (y = *xb >> 16) { 334 x = xa; 335 xc = xc0; 336 carry = 0; 337 z2 = *xc; 338 do { 339 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 340 carry = z >> 16; 341 Storeinc(xc, z, z2); 342 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 343 carry = z2 >> 16; 344 } while (x < xae); 345 *xc = z2; 346 } 347 } 348 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) 349 ; 350 c->wds = wc; 351 return c; 352 } 353 354 static Bigint *p5s; 355 356 static Bigint * 357 pow5mult(Bigint *b, int k) 358 { 359 Bigint * b1, *p5, *p51; 360 int i; 361 static int p05[3] = { 362 5, 25, 125 }; 363 364 if (i = k & 3) 365 b = multadd(b, p05[i-1], 0); 366 367 if (!(k >>= 2)) 368 return b; 369 if (!(p5 = p5s)) { 370 /* first time */ 371 ACQUIRE_DTOA_LOCK(1); 372 if (!(p5 = p5s)) { 373 p5 = p5s = i2b(625); 374 p5->next = 0; 375 } 376 FREE_DTOA_LOCK(1); 377 } 378 for (; ; ) { 379 if (k & 1) { 380 b1 = mult(b, p5); 381 Bfree(b); 382 b = b1; 383 } 384 if (!(k >>= 1)) 385 break; 386 if (!(p51 = p5->next)) { 387 ACQUIRE_DTOA_LOCK(1); 388 if (!(p51 = p5->next)) { 389 p51 = p5->next = mult(p5, p5); 390 p51->next = 0; 391 } 392 FREE_DTOA_LOCK(1); 393 } 394 p5 = p51; 395 } 396 return b; 397 } 398 399 static Bigint * 400 lshift(Bigint *b, int k) 401 { 402 int i, k1, n, n1; 403 Bigint * b1; 404 unsigned int * x, *x1, *xe, z; 405 406 n = k >> 5; 407 k1 = b->k; 408 n1 = n + b->wds + 1; 409 for (i = b->maxwds; n1 > i; i <<= 1) 410 k1++; 411 b1 = Balloc(k1); 412 x1 = b1->x; 413 for (i = 0; i < n; i++) 414 *x1++ = 0; 415 x = b->x; 416 xe = x + b->wds; 417 if (k &= 0x1f) { 418 k1 = 32 - k; 419 z = 0; 420 do { 421 *x1++ = *x << k | z; 422 z = *x++ >> k1; 423 } while (x < xe); 424 if (*x1 = z) 425 ++n1; 426 } else 427 do 428 *x1++ = *x++; 429 while (x < xe); 430 b1->wds = n1 - 1; 431 Bfree(b); 432 return b1; 433 } 434 435 static int 436 cmp(Bigint *a, Bigint *b) 437 { 438 unsigned int * xa, *xa0, *xb, *xb0; 439 int i, j; 440 441 i = a->wds; 442 j = b->wds; 443 if (i -= j) 444 return i; 445 xa0 = a->x; 446 xa = xa0 + j; 447 xb0 = b->x; 448 xb = xb0 + j; 449 for (; ; ) { 450 if (*--xa != *--xb) 451 return * xa < *xb ? -1 : 1; 452 if (xa <= xa0) 453 break; 454 } 455 return 0; 456 } 457 458 static Bigint * 459 diff(Bigint *a, Bigint *b) 460 { 461 Bigint * c; 462 int i, wa, wb; 463 unsigned int * xa, *xae, *xb, *xbe, *xc; 464 unsigned int borrow, y; 465 unsigned int z; 466 467 i = cmp(a, b); 468 if (!i) { 469 c = Balloc(0); 470 c->wds = 1; 471 c->x[0] = 0; 472 return c; 473 } 474 if (i < 0) { 475 c = a; 476 a = b; 477 b = c; 478 i = 1; 479 } else 480 i = 0; 481 c = Balloc(a->k); 482 c->sign = i; 483 wa = a->wds; 484 xa = a->x; 485 xae = xa + wa; 486 wb = b->wds; 487 xb = b->x; 488 xbe = xb + wb; 489 xc = c->x; 490 borrow = 0; 491 do { 492 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 493 borrow = (y & 0x10000) >> 16; 494 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 495 borrow = (z & 0x10000) >> 16; 496 Storeinc(xc, z, y); 497 } while (xb < xbe); 498 while (xa < xae) { 499 y = (*xa & 0xffff) - borrow; 500 borrow = (y & 0x10000) >> 16; 501 z = (*xa++ >> 16) - borrow; 502 borrow = (z & 0x10000) >> 16; 503 Storeinc(xc, z, y); 504 } 505 while (!*--xc) 506 wa--; 507 c->wds = wa; 508 return c; 509 } 510 511 static double 512 ulp(double x) 513 { 514 ulong L; 515 Ulongs uls; 516 517 uls = double2ulongs(x); 518 L = (uls.hi & Exp_mask) - (P - 1) * Exp_msk1; 519 return ulongs2double((Ulongs){L, 0}); 520 } 521 522 static double 523 b2d(Bigint *a, int *e) 524 { 525 unsigned *xa, *xa0, w, y, z; 526 int k; 527 ulong d0, d1; 528 529 xa0 = a->x; 530 xa = xa0 + a->wds; 531 y = *--xa; 532 k = hi0bits(y); 533 *e = 32 - k; 534 if (k < Ebits) { 535 w = xa > xa0 ? *--xa : 0; 536 d1 = y << (32 - Ebits) + k | w >> Ebits - k; 537 return ulongs2double((Ulongs){Exp_1 | y >> Ebits - k, d1}); 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 return ulongs2double((Ulongs){d0, d1}); 549 } 550 551 static Bigint * 552 d2b(double d, int *e, int *bits) 553 { 554 Bigint * b; 555 int de, i, k; 556 unsigned *x, y, z; 557 Ulongs uls; 558 559 b = Balloc(1); 560 x = b->x; 561 562 uls = double2ulongs(d); 563 z = uls.hi & Frac_mask; 564 uls.hi &= 0x7fffffff; /* clear sign bit, which we ignore */ 565 de = (int)(uls.hi >> Exp_shift); 566 z |= Exp_msk11; 567 if (y = uls.lo) { /* assignment = */ 568 if (k = lo0bits(&y)) { /* assignment = */ 569 x[0] = y | z << 32 - k; 570 z >>= k; 571 } else 572 x[0] = y; 573 i = b->wds = (x[1] = z) ? 2 : 1; 574 } else { 575 k = lo0bits(&z); 576 x[0] = z; 577 i = b->wds = 1; 578 k += 32; 579 } 580 USED(i); 581 *e = de - Bias - (P - 1) + k; 582 *bits = P - k; 583 return b; 584 } 585 586 static double 587 ratio(Bigint *a, Bigint *b) 588 { 589 double da, db; 590 int k, ka, kb; 591 Ulongs uls; 592 593 da = b2d(a, &ka); 594 db = b2d(b, &kb); 595 k = ka - kb + 32 * (a->wds - b->wds); 596 if (k > 0) { 597 uls = double2ulongs(da); 598 uls.hi += k * Exp_msk1; 599 da = ulongs2double(uls); 600 } else { 601 k = -k; 602 uls = double2ulongs(db); 603 uls.hi += k * Exp_msk1; 604 db = ulongs2double(uls); 605 } 606 return da / db; 607 } 608 609 static const double 610 tens[] = { 611 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 612 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 613 1e20, 1e21, 1e22 614 }; 615 616 static const double 617 bigtens[] = { 618 1e16, 1e32, 1e64, 1e128, 1e256 }; 619 620 static const double tinytens[] = { 621 1e-16, 1e-32, 1e-64, 1e-128, 622 9007199254740992.e-256 623 }; 624 625 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 626 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 627 #define Scale_Bit 0x10 628 #define n_bigtens 5 629 630 #define NAN_WORD0 0x7ff80000 631 632 #define NAN_WORD1 0 633 634 static int 635 match(const char **sp, char *t) 636 { 637 int c, d; 638 const char * s = *sp; 639 640 while (d = *t++) { 641 if ((c = *++s) >= 'A' && c <= 'Z') 642 c += 'a' - 'A'; 643 if (c != d) 644 return 0; 645 } 646 *sp = s + 1; 647 return 1; 648 } 649 650 static void 651 gethex(double *rvp, const char **sp) 652 { 653 unsigned int c, x[2]; 654 const char * s; 655 int havedig, udx0, xshift; 656 657 x[0] = x[1] = 0; 658 havedig = xshift = 0; 659 udx0 = 1; 660 s = *sp; 661 while (c = *(const unsigned char * )++s) { 662 if (c >= '0' && c <= '9') 663 c -= '0'; 664 else if (c >= 'a' && c <= 'f') 665 c += 10 - 'a'; 666 else if (c >= 'A' && c <= 'F') 667 c += 10 - 'A'; 668 else if (c <= ' ') { 669 if (udx0 && havedig) { 670 udx0 = 0; 671 xshift = 1; 672 } 673 continue; 674 } else if (/*(*/ c == ')') { 675 *sp = s + 1; 676 break; 677 } else 678 return; /* invalid form: don't change *sp */ 679 havedig = 1; 680 if (xshift) { 681 xshift = 0; 682 x[0] = x[1]; 683 x[1] = 0; 684 } 685 if (udx0) 686 x[0] = (x[0] << 4) | (x[1] >> 28); 687 x[1] = (x[1] << 4) | c; 688 } 689 if ((x[0] &= 0xfffff) || x[1]) 690 *rvp = ulongs2double((Ulongs){Exp_mask | x[0], x[1]}); 691 } 692 693 static int 694 quorem(Bigint *b, Bigint *S) 695 { 696 int n; 697 unsigned int * bx, *bxe, q, *sx, *sxe; 698 unsigned int borrow, carry, y, ys; 699 unsigned int si, z, zs; 700 701 n = S->wds; 702 if (b->wds < n) 703 return 0; 704 sx = S->x; 705 sxe = sx + --n; 706 bx = b->x; 707 bxe = bx + n; 708 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 709 if (q) { 710 borrow = 0; 711 carry = 0; 712 do { 713 si = *sx++; 714 ys = (si & 0xffff) * q + carry; 715 zs = (si >> 16) * q + (ys >> 16); 716 carry = zs >> 16; 717 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 718 borrow = (y & 0x10000) >> 16; 719 z = (*bx >> 16) - (zs & 0xffff) - borrow; 720 borrow = (z & 0x10000) >> 16; 721 Storeinc(bx, z, y); 722 } while (sx <= sxe); 723 if (!*bxe) { 724 bx = b->x; 725 while (--bxe > bx && !*bxe) 726 --n; 727 b->wds = n; 728 } 729 } 730 if (cmp(b, S) >= 0) { 731 q++; 732 borrow = 0; 733 carry = 0; 734 bx = b->x; 735 sx = S->x; 736 do { 737 si = *sx++; 738 ys = (si & 0xffff) + carry; 739 zs = (si >> 16) + (ys >> 16); 740 carry = zs >> 16; 741 y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 742 borrow = (y & 0x10000) >> 16; 743 z = (*bx >> 16) - (zs & 0xffff) - borrow; 744 borrow = (z & 0x10000) >> 16; 745 Storeinc(bx, z, y); 746 } while (sx <= sxe); 747 bx = b->x; 748 bxe = bx + n; 749 if (!*bxe) { 750 while (--bxe > bx && !*bxe) 751 --n; 752 b->wds = n; 753 } 754 } 755 return q; 756 } 757 758 static char * 759 rv_alloc(int i) 760 { 761 int j, k, *r; 762 763 j = sizeof(unsigned int); 764 for (k = 0; 765 sizeof(Bigint) - sizeof(unsigned int) - sizeof(int) + j <= i; 766 j <<= 1) 767 k++; 768 r = (int * )Balloc(k); 769 *r = k; 770 return 771 (char *)(r + 1); 772 } 773 774 static char * 775 nrv_alloc(char *s, char **rve, int n) 776 { 777 char *rv, *t; 778 779 t = rv = rv_alloc(n); 780 while (*t = *s++) 781 t++; 782 if (rve) 783 *rve = t; 784 return rv; 785 } 786 787 /* freedtoa(s) must be used to free values s returned by dtoa 788 * when MULTIPLE_THREADS is #defined. It should be used in all cases, 789 * but for consistency with earlier versions of dtoa, it is optional 790 * when MULTIPLE_THREADS is not defined. 791 */ 792 793 void 794 freedtoa(char *s) 795 { 796 Bigint * b = (Bigint * )((int *)s - 1); 797 b->maxwds = 1 << (b->k = *(int * )b); 798 Bfree(b); 799 } 800 801 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 802 * 803 * Inspired by "How to Print Floating-Point Numbers Accurately" by 804 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 805 * 806 * Modifications: 807 * 1. Rather than iterating, we use a simple numeric overestimate 808 * to determine k = floor(log10(d)). We scale relevant 809 * quantities using O(log2(k)) rather than O(k) multiplications. 810 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 811 * try to generate digits strictly left to right. Instead, we 812 * compute with fewer bits and propagate the carry if necessary 813 * when rounding the final digit up. This is often faster. 814 * 3. Under the assumption that input will be rounded nearest, 815 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 816 * That is, we allow equality in stopping tests when the 817 * round-nearest rule will give the same floating-point value 818 * as would satisfaction of the stopping test with strict 819 * inequality. 820 * 4. We remove common factors of powers of 2 from relevant 821 * quantities. 822 * 5. When converting floating-point integers less than 1e16, 823 * we use floating-point arithmetic rather than resorting 824 * to multiple-precision integers. 825 * 6. When asked to produce fewer than 15 digits, we first try 826 * to get by with floating-point arithmetic; we resort to 827 * multiple-precision integer arithmetic only if we cannot 828 * guarantee that the floating-point calculation has given 829 * the correctly rounded result. For k requested digits and 830 * "uniformly" distributed input, the probability is 831 * something like 10^(k-15) that we must resort to the int 832 * calculation. 833 */ 834 835 char * 836 dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve) 837 { 838 /* Arguments ndigits, decpt, sign are similar to those 839 of ecvt and fcvt; trailing zeros are suppressed from 840 the returned string. If not null, *rve is set to point 841 to the end of the return value. If d is +-Infinity or NaN, 842 then *decpt is set to 9999. 843 844 mode: 845 0 ==> shortest string that yields d when read in 846 and rounded to nearest. 847 1 ==> like 0, but with Steele & White stopping rule; 848 e.g. with IEEE P754 arithmetic , mode 0 gives 849 1e23 whereas mode 1 gives 9.999999999999999e22. 850 2 ==> max(1,ndigits) significant digits. This gives a 851 return value similar to that of ecvt, except 852 that trailing zeros are suppressed. 853 3 ==> through ndigits past the decimal point. This 854 gives a return value similar to that from fcvt, 855 except that trailing zeros are suppressed, and 856 ndigits can be negative. 857 4-9 should give the same return values as 2-3, i.e., 858 4 <= mode <= 9 ==> same return as mode 859 2 + (mode & 1). These modes are mainly for 860 debugging; often they run slower but sometimes 861 faster than modes 2-3. 862 4,5,8,9 ==> left-to-right digit generation. 863 6-9 ==> don't try fast floating-point estimate 864 (if applicable). 865 866 Values of mode other than 0-9 are treated as mode 0. 867 868 Sufficient space is allocated to the return value 869 to hold the suppressed trailing zeros. 870 */ 871 872 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 873 j, j1, k, k0, k_check, L, leftright, m2, m5, s2, s5, 874 spec_case, try_quick; 875 Bigint * b, *b1, *delta, *mlo=nil, *mhi, *S; 876 double d2, ds, eps; 877 char *s, *s0; 878 Ulongs ulsd, ulsd2; 879 880 ulsd = double2ulongs(d); 881 if (ulsd.hi & Sign_bit) { 882 /* set sign for everything, including 0's and NaNs */ 883 *sign = 1; 884 ulsd.hi &= ~Sign_bit; /* clear sign bit */ 885 } else 886 *sign = 0; 887 888 if ((ulsd.hi & Exp_mask) == Exp_mask) { 889 /* Infinity or NaN */ 890 *decpt = 9999; 891 if (!ulsd.lo && !(ulsd.hi & 0xfffff)) 892 return nrv_alloc("Infinity", rve, 8); 893 return nrv_alloc("NaN", rve, 3); 894 } 895 d = ulongs2double(ulsd); 896 897 if (!d) { 898 *decpt = 1; 899 return nrv_alloc("0", rve, 1); 900 } 901 902 b = d2b(d, &be, &bbits); 903 i = (int)(ulsd.hi >> Exp_shift1 & (Exp_mask >> Exp_shift1)); 904 905 ulsd2 = ulsd; 906 ulsd2.hi &= Frac_mask1; 907 ulsd2.hi |= Exp_11; 908 d2 = ulongs2double(ulsd2); 909 910 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 911 * log10(x) = log(x) / log(10) 912 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 913 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 914 * 915 * This suggests computing an approximation k to log10(d) by 916 * 917 * k = (i - Bias)*0.301029995663981 918 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 919 * 920 * We want k to be too large rather than too small. 921 * The error in the first-order Taylor series approximation 922 * is in our favor, so we just round up the constant enough 923 * to compensate for any error in the multiplication of 924 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 925 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 926 * adding 1e-13 to the constant term more than suffices. 927 * Hence we adjust the constant term to 0.1760912590558. 928 * (We could get a more accurate k by invoking log10, 929 * but this is probably not worthwhile.) 930 */ 931 932 i -= Bias; 933 ds = (d2 - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981; 934 k = (int)ds; 935 if (ds < 0. && ds != k) 936 k--; /* want k = floor(ds) */ 937 k_check = 1; 938 if (k >= 0 && k <= Ten_pmax) { 939 if (d < tens[k]) 940 k--; 941 k_check = 0; 942 } 943 j = bbits - i - 1; 944 if (j >= 0) { 945 b2 = 0; 946 s2 = j; 947 } else { 948 b2 = -j; 949 s2 = 0; 950 } 951 if (k >= 0) { 952 b5 = 0; 953 s5 = k; 954 s2 += k; 955 } else { 956 b2 -= k; 957 b5 = -k; 958 s5 = 0; 959 } 960 if (mode < 0 || mode > 9) 961 mode = 0; 962 try_quick = 1; 963 if (mode > 5) { 964 mode -= 4; 965 try_quick = 0; 966 } 967 leftright = 1; 968 switch (mode) { 969 case 0: 970 case 1: 971 default: 972 ilim = ilim1 = -1; 973 i = 18; 974 ndigits = 0; 975 break; 976 case 2: 977 leftright = 0; 978 /* no break */ 979 case 4: 980 if (ndigits <= 0) 981 ndigits = 1; 982 ilim = ilim1 = i = ndigits; 983 break; 984 case 3: 985 leftright = 0; 986 /* no break */ 987 case 5: 988 i = ndigits + k + 1; 989 ilim = i; 990 ilim1 = i - 1; 991 if (i <= 0) 992 i = 1; 993 } 994 s = s0 = rv_alloc(i); 995 996 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 997 998 /* Try to get by with floating-point arithmetic. */ 999 1000 i = 0; 1001 d2 = d; 1002 k0 = k; 1003 ilim0 = ilim; 1004 ieps = 2; /* conservative */ 1005 if (k > 0) { 1006 ds = tens[k&0xf]; 1007 j = k >> 4; 1008 if (j & Bletch) { 1009 /* prevent overflows */ 1010 j &= Bletch - 1; 1011 d /= bigtens[n_bigtens-1]; 1012 ieps++; 1013 } 1014 for (; j; j >>= 1, i++) 1015 if (j & 1) { 1016 ieps++; 1017 ds *= bigtens[i]; 1018 } 1019 d /= ds; 1020 } else if (j1 = -k) { 1021 d *= tens[j1 & 0xf]; 1022 for (j = j1 >> 4; j; j >>= 1, i++) 1023 if (j & 1) { 1024 ieps++; 1025 d *= bigtens[i]; 1026 } 1027 } 1028 if (k_check && d < 1. && ilim > 0) { 1029 if (ilim1 <= 0) 1030 goto fast_failed; 1031 ilim = ilim1; 1032 k--; 1033 d *= 10.; 1034 ieps++; 1035 } 1036 eps = ieps * d + 7.; 1037 1038 ulsd = double2ulongs(eps); 1039 ulsd.hi -= (P - 1) * Exp_msk1; 1040 eps = ulongs2double(ulsd); 1041 1042 if (ilim == 0) { 1043 S = mhi = 0; 1044 d -= 5.; 1045 if (d > eps) 1046 goto one_digit; 1047 if (d < -eps) 1048 goto no_digits; 1049 goto fast_failed; 1050 } 1051 /* Generate ilim digits, then fix them up. */ 1052 eps *= tens[ilim-1]; 1053 for (i = 1; ; i++, d *= 10.) { 1054 L = d; 1055 // assert(L < 10); 1056 d -= L; 1057 *s++ = '0' + (int)L; 1058 if (i == ilim) { 1059 if (d > 0.5 + eps) 1060 goto bump_up; 1061 else if (d < 0.5 - eps) { 1062 while (*--s == '0') 1063 ; 1064 s++; 1065 goto ret1; 1066 } 1067 break; 1068 } 1069 } 1070 fast_failed: 1071 s = s0; 1072 d = d2; 1073 k = k0; 1074 ilim = ilim0; 1075 } 1076 1077 /* Do we have a "small" integer? */ 1078 1079 if (be >= 0 && k <= Int_max) { 1080 /* Yes. */ 1081 ds = tens[k]; 1082 if (ndigits < 0 && ilim <= 0) { 1083 S = mhi = 0; 1084 if (ilim < 0 || d <= 5 * ds) 1085 goto no_digits; 1086 goto one_digit; 1087 } 1088 for (i = 1; ; i++) { 1089 L = d / ds; 1090 d -= L * ds; 1091 *s++ = '0' + (int)L; 1092 if (i == ilim) { 1093 d += d; 1094 if (d > ds || d == ds && L & 1) { 1095 bump_up: 1096 while (*--s == '9') 1097 if (s == s0) { 1098 k++; 1099 *s = '0'; 1100 break; 1101 } 1102 ++ * s++; 1103 } 1104 break; 1105 } 1106 if (!(d *= 10.)) 1107 break; 1108 } 1109 goto ret1; 1110 } 1111 1112 m2 = b2; 1113 m5 = b5; 1114 mhi = mlo = 0; 1115 if (leftright) { 1116 if (mode < 2) { 1117 i = 1118 1 + P - bbits; 1119 } else { 1120 j = ilim - 1; 1121 if (m5 >= j) 1122 m5 -= j; 1123 else { 1124 s5 += j -= m5; 1125 b5 += j; 1126 m5 = 0; 1127 } 1128 if ((i = ilim) < 0) { 1129 m2 -= i; 1130 i = 0; 1131 } 1132 } 1133 b2 += i; 1134 s2 += i; 1135 mhi = i2b(1); 1136 } 1137 if (m2 > 0 && s2 > 0) { 1138 i = m2 < s2 ? m2 : s2; 1139 b2 -= i; 1140 m2 -= i; 1141 s2 -= i; 1142 } 1143 if (b5 > 0) { 1144 if (leftright) { 1145 if (m5 > 0) { 1146 mhi = pow5mult(mhi, m5); 1147 b1 = mult(mhi, b); 1148 Bfree(b); 1149 b = b1; 1150 } 1151 if (j = b5 - m5) 1152 b = pow5mult(b, j); 1153 } else 1154 b = pow5mult(b, b5); 1155 } 1156 S = i2b(1); 1157 if (s5 > 0) 1158 S = pow5mult(S, s5); 1159 1160 /* Check for special case that d is a normalized power of 2. */ 1161 1162 spec_case = 0; 1163 if (mode < 2) { 1164 ulsd = double2ulongs(d); 1165 if (!ulsd.lo && !(ulsd.hi & Bndry_mask)) { 1166 /* The special case */ 1167 b2 += Log2P; 1168 s2 += Log2P; 1169 spec_case = 1; 1170 } 1171 } 1172 1173 /* Arrange for convenient computation of quotients: 1174 * shift left if necessary so divisor has 4 leading 0 bits. 1175 * 1176 * Perhaps we should just compute leading 28 bits of S once 1177 * and for all and pass them and a shift to quorem, so it 1178 * can do shifts and ors to compute the numerator for q. 1179 */ 1180 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) 1181 i = 32 - i; 1182 if (i > 4) { 1183 i -= 4; 1184 b2 += i; 1185 m2 += i; 1186 s2 += i; 1187 } else if (i < 4) { 1188 i += 28; 1189 b2 += i; 1190 m2 += i; 1191 s2 += i; 1192 } 1193 if (b2 > 0) 1194 b = lshift(b, b2); 1195 if (s2 > 0) 1196 S = lshift(S, s2); 1197 if (k_check) { 1198 if (cmp(b, S) < 0) { 1199 k--; 1200 b = multadd(b, 10, 0); /* we botched the k estimate */ 1201 if (leftright) 1202 mhi = multadd(mhi, 10, 0); 1203 ilim = ilim1; 1204 } 1205 } 1206 if (ilim <= 0 && mode > 2) { 1207 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) { 1208 /* no digits, fcvt style */ 1209 no_digits: 1210 k = -1 - ndigits; 1211 goto ret; 1212 } 1213 one_digit: 1214 *s++ = '1'; 1215 k++; 1216 goto ret; 1217 } 1218 if (leftright) { 1219 if (m2 > 0) 1220 mhi = lshift(mhi, m2); 1221 1222 /* Compute mlo -- check for special case 1223 * that d is a normalized power of 2. 1224 */ 1225 1226 mlo = mhi; 1227 if (spec_case) { 1228 mhi = Balloc(mhi->k); 1229 Bcopy(mhi, mlo); 1230 mhi = lshift(mhi, Log2P); 1231 } 1232 1233 for (i = 1; ; i++) { 1234 dig = quorem(b, S) + '0'; 1235 /* Do we yet have the shortest decimal string 1236 * that will round to d? 1237 */ 1238 j = cmp(b, mlo); 1239 delta = diff(S, mhi); 1240 j1 = delta->sign ? 1 : cmp(b, delta); 1241 Bfree(delta); 1242 ulsd = double2ulongs(d); 1243 if (j1 == 0 && !mode && !(ulsd.lo & 1)) { 1244 if (dig == '9') 1245 goto round_9_up; 1246 if (j > 0) 1247 dig++; 1248 *s++ = dig; 1249 goto ret; 1250 } 1251 if (j < 0 || j == 0 && !mode && !(ulsd.lo & 1)) { 1252 if (j1 > 0) { 1253 b = lshift(b, 1); 1254 j1 = cmp(b, S); 1255 if ((j1 > 0 || j1 == 0 && dig & 1) 1256 && dig++ == '9') 1257 goto round_9_up; 1258 } 1259 *s++ = dig; 1260 goto ret; 1261 } 1262 if (j1 > 0) { 1263 if (dig == '9') { /* possible if i == 1 */ 1264 round_9_up: 1265 *s++ = '9'; 1266 goto roundoff; 1267 } 1268 *s++ = dig + 1; 1269 goto ret; 1270 } 1271 *s++ = dig; 1272 if (i == ilim) 1273 break; 1274 b = multadd(b, 10, 0); 1275 if (mlo == mhi) 1276 mlo = mhi = multadd(mhi, 10, 0); 1277 else { 1278 mlo = multadd(mlo, 10, 0); 1279 mhi = multadd(mhi, 10, 0); 1280 } 1281 } 1282 } else 1283 for (i = 1; ; i++) { 1284 *s++ = dig = quorem(b, S) + '0'; 1285 if (i >= ilim) 1286 break; 1287 b = multadd(b, 10, 0); 1288 } 1289 1290 /* Round off last digit */ 1291 1292 b = lshift(b, 1); 1293 j = cmp(b, S); 1294 if (j > 0 || j == 0 && dig & 1) { 1295 roundoff: 1296 while (*--s == '9') 1297 if (s == s0) { 1298 k++; 1299 *s++ = '1'; 1300 goto ret; 1301 } 1302 ++ * s++; 1303 } else { 1304 while (*--s == '0') 1305 ; 1306 s++; 1307 } 1308 ret: 1309 Bfree(S); 1310 if (mhi) { 1311 if (mlo && mlo != mhi) 1312 Bfree(mlo); 1313 Bfree(mhi); 1314 } 1315 ret1: 1316 Bfree(b); 1317 *s = 0; 1318 *decpt = k + 1; 1319 if (rve) 1320 *rve = s; 1321 return s0; 1322 } 1323