1 #include "fconv.h" 2 3 static int quorem(Bigint *, Bigint *); 4 5 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 6 * 7 * Inspired by "How to Print Floating-Point Numbers Accurately" by 8 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. 9 * 10 * Modifications: 11 * 1. Rather than iterating, we use a simple numeric overestimate 12 * to determine k = floor(log10(d)). We scale relevant 13 * quantities using O(log2(k)) rather than O(k) multiplications. 14 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 15 * try to generate digits strictly left to right. Instead, we 16 * compute with fewer bits and propagate the carry if necessary 17 * when rounding the final digit up. This is often faster. 18 * 3. Under the assumption that input will be rounded nearest, 19 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 20 * That is, we allow equality in stopping tests when the 21 * round-nearest rule will give the same floating-point value 22 * as would satisfaction of the stopping test with strict 23 * inequality. 24 * 4. We remove common factors of powers of 2 from relevant 25 * quantities. 26 * 5. When converting floating-point integers less than 1e16, 27 * we use floating-point arithmetic rather than resorting 28 * to multiple-precision integers. 29 * 6. When asked to produce fewer than 15 digits, we first try 30 * to get by with floating-point arithmetic; we resort to 31 * multiple-precision integer arithmetic only if we cannot 32 * guarantee that the floating-point calculation has given 33 * the correctly rounded result. For k requested digits and 34 * "uniformly" distributed input, the probability is 35 * something like 10^(k-15) that we must resort to the long 36 * calculation. 37 */ 38 39 char * 40 _dtoa(double darg, int mode, int ndigits, int *decpt, int *sign, char **rve) 41 { 42 /* Arguments ndigits, decpt, sign are similar to those 43 of ecvt and fcvt; trailing zeros are suppressed from 44 the returned string. If not null, *rve is set to point 45 to the end of the return value. If d is +-Infinity or NaN, 46 then *decpt is set to 9999. 47 48 mode: 49 0 ==> shortest string that yields d when read in 50 and rounded to nearest. 51 1 ==> like 0, but with Steele & White stopping rule; 52 e.g. with IEEE P754 arithmetic , mode 0 gives 53 1e23 whereas mode 1 gives 9.999999999999999e22. 54 2 ==> max(1,ndigits) significant digits. This gives a 55 return value similar to that of ecvt, except 56 that trailing zeros are suppressed. 57 3 ==> through ndigits past the decimal point. This 58 gives a return value similar to that from fcvt, 59 except that trailing zeros are suppressed, and 60 ndigits can be negative. 61 4-9 should give the same return values as 2-3, i.e., 62 4 <= mode <= 9 ==> same return as mode 63 2 + (mode & 1). These modes are mainly for 64 debugging; often they run slower but sometimes 65 faster than modes 2-3. 66 4,5,8,9 ==> left-to-right digit generation. 67 6-9 ==> don't try fast floating-point estimate 68 (if applicable). 69 70 Values of mode other than 0-9 are treated as mode 0. 71 72 Sufficient space is allocated to the return value 73 to hold the suppressed trailing zeros. 74 */ 75 76 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 77 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 78 spec_case, try_quick; 79 long L; 80 #ifndef Sudden_Underflow 81 int denorm; 82 unsigned long x; 83 #endif 84 Bigint *b, *b1, *delta, *mlo, *mhi, *S; 85 double ds; 86 Dul d2, eps; 87 char *s, *s0; 88 static Bigint *result; 89 static int result_k; 90 Dul d; 91 92 d.d = darg; 93 if (result) { 94 result->k = result_k; 95 result->maxwds = 1 << result_k; 96 Bfree(result); 97 result = 0; 98 } 99 100 if (word0(d) & Sign_bit) { 101 /* set sign for everything, including 0's and NaNs */ 102 *sign = 1; 103 word0(d) &= ~Sign_bit; /* clear sign bit */ 104 } 105 else 106 *sign = 0; 107 108 #if defined(IEEE_Arith) + defined(VAX) 109 #ifdef IEEE_Arith 110 if ((word0(d) & Exp_mask) == Exp_mask) 111 #else 112 if (word0(d) == 0x8000) 113 #endif 114 { 115 /* Infinity or NaN */ 116 *decpt = 9999; 117 s = 118 #ifdef IEEE_Arith 119 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : 120 #endif 121 "NaN"; 122 if (rve) 123 *rve = 124 #ifdef IEEE_Arith 125 s[3] ? s + 8 : 126 #endif 127 s + 3; 128 return s; 129 } 130 #endif 131 #ifdef IBM 132 d.d += 0; /* normalize */ 133 #endif 134 if (!d.d) { 135 *decpt = 1; 136 s = "0"; 137 if (rve) 138 *rve = s + 1; 139 return s; 140 } 141 142 b = d2b(d.d, &be, &bbits); 143 #ifdef Sudden_Underflow 144 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); 145 #else 146 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) { 147 #endif 148 d2.d = d.d; 149 word0(d2) &= Frac_mask1; 150 word0(d2) |= Exp_11; 151 #ifdef IBM 152 if (j = 11 - hi0bits(word0(d2) & Frac_mask)) 153 d2.d /= 1 << j; 154 #endif 155 156 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 157 * log10(x) = log(x) / log(10) 158 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 159 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 160 * 161 * This suggests computing an approximation k to log10(d) by 162 * 163 * k = (i - Bias)*0.301029995663981 164 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 165 * 166 * We want k to be too large rather than too small. 167 * The error in the first-order Taylor series approximation 168 * is in our favor, so we just round up the constant enough 169 * to compensate for any error in the multiplication of 170 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 171 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 172 * adding 1e-13 to the constant term more than suffices. 173 * Hence we adjust the constant term to 0.1760912590558. 174 * (We could get a more accurate k by invoking log10, 175 * but this is probably not worthwhile.) 176 */ 177 178 i -= Bias; 179 #ifdef IBM 180 i <<= 2; 181 i += j; 182 #endif 183 #ifndef Sudden_Underflow 184 denorm = 0; 185 } 186 else { 187 /* d is denormalized */ 188 189 i = bbits + be + (Bias + (P-1) - 1); 190 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32 191 : word1(d) << 32 - i; 192 d2.d = x; 193 word0(d2) -= 31*Exp_msk1; /* adjust exponent */ 194 i -= (Bias + (P-1) - 1) + 1; 195 denorm = 1; 196 } 197 #endif 198 ds = (d2.d-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 199 k = floor(ds); 200 k_check = 1; 201 if (k >= 0 && k <= Ten_pmax) { 202 if (d.d < tens[k]) 203 k--; 204 k_check = 0; 205 } 206 j = bbits - i - 1; 207 if (j >= 0) { 208 b2 = 0; 209 s2 = j; 210 } 211 else { 212 b2 = -j; 213 s2 = 0; 214 } 215 if (k >= 0) { 216 b5 = 0; 217 s5 = k; 218 s2 += k; 219 } 220 else { 221 b2 -= k; 222 b5 = -k; 223 s5 = 0; 224 } 225 if (mode < 0 || mode > 9) 226 mode = 0; 227 try_quick = 1; 228 if (mode > 5) { 229 mode -= 4; 230 try_quick = 0; 231 } 232 leftright = 1; 233 switch(mode) { 234 case 0: 235 case 1: 236 ilim = ilim1 = -1; 237 i = 18; 238 ndigits = 0; 239 break; 240 case 2: 241 leftright = 0; 242 /* no break */ 243 case 4: 244 if (ndigits <= 0) 245 ndigits = 1; 246 ilim = ilim1 = i = ndigits; 247 break; 248 case 3: 249 leftright = 0; 250 /* no break */ 251 case 5: 252 i = ndigits + k + 1; 253 ilim = i; 254 ilim1 = i - 1; 255 if (i <= 0) 256 i = 1; 257 } 258 j = sizeof(unsigned long); 259 for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j <= i; 260 j <<= 1) result_k++; 261 result = Balloc(result_k); 262 s = s0 = (char *)result; 263 264 if (ilim >= 0 && ilim <= Quick_max && try_quick) { 265 266 /* Try to get by with floating-point arithmetic. */ 267 268 i = 0; 269 d2.d = d.d; 270 k0 = k; 271 ilim0 = ilim; 272 ieps = 2; /* conservative */ 273 if (k > 0) { 274 ds = tens[k&0xf]; 275 j = k >> 4; 276 if (j & Bletch) { 277 /* prevent overflows */ 278 j &= Bletch - 1; 279 d.d /= bigtens[n_bigtens-1]; 280 ieps++; 281 } 282 for(; j; j >>= 1, i++) 283 if (j & 1) { 284 ieps++; 285 ds *= bigtens[i]; 286 } 287 d.d /= ds; 288 } 289 else if (j1 = -k) { 290 d.d *= tens[j1 & 0xf]; 291 for(j = j1 >> 4; j; j >>= 1, i++) 292 if (j & 1) { 293 ieps++; 294 d.d *= bigtens[i]; 295 } 296 } 297 if (k_check && d.d < 1. && ilim > 0) { 298 if (ilim1 <= 0) 299 goto fast_failed; 300 ilim = ilim1; 301 k--; 302 d.d *= 10.; 303 ieps++; 304 } 305 eps.d = ieps*d.d + 7.; 306 word0(eps) -= (P-1)*Exp_msk1; 307 if (ilim == 0) { 308 S = mhi = 0; 309 d.d -= 5.; 310 if (d.d > eps.d) 311 goto one_digit; 312 if (d.d < -eps.d) 313 goto no_digits; 314 goto fast_failed; 315 } 316 #ifndef No_leftright 317 if (leftright) { 318 /* Use Steele & White method of only 319 * generating digits needed. 320 */ 321 eps.d = 0.5/tens[ilim-1] - eps.d; 322 for(i = 0;;) { 323 L = floor(d.d); 324 d.d -= L; 325 *s++ = '0' + (int)L; 326 if (d.d < eps.d) 327 goto ret1; 328 if (1. - d.d < eps.d) 329 goto bump_up; 330 if (++i >= ilim) 331 break; 332 eps.d *= 10.; 333 d.d *= 10.; 334 } 335 } 336 else { 337 #endif 338 /* Generate ilim digits, then fix them up. */ 339 eps.d *= tens[ilim-1]; 340 for(i = 1;; i++, d.d *= 10.) { 341 L = floor(d.d); 342 d.d -= L; 343 *s++ = '0' + (int)L; 344 if (i == ilim) { 345 if (d.d > 0.5 + eps.d) 346 goto bump_up; 347 else if (d.d < 0.5 - eps.d) { 348 while(*--s == '0'); 349 s++; 350 goto ret1; 351 } 352 break; 353 } 354 } 355 #ifndef No_leftright 356 } 357 #endif 358 fast_failed: 359 s = s0; 360 d.d = d2.d; 361 k = k0; 362 ilim = ilim0; 363 } 364 365 /* Do we have a "small" integer? */ 366 367 if (be >= 0 && k <= Int_max) { 368 /* Yes. */ 369 ds = tens[k]; 370 if (ndigits < 0 && ilim <= 0) { 371 S = mhi = 0; 372 if (ilim < 0 || d.d <= 5*ds) 373 goto no_digits; 374 goto one_digit; 375 } 376 for(i = 1;; i++) { 377 L = floor(d.d / ds); 378 d.d -= L*ds; 379 #ifdef Check_FLT_ROUNDS 380 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 381 if (d.d < 0) { 382 L--; 383 d.d += ds; 384 } 385 #endif 386 *s++ = '0' + (int)L; 387 if (i == ilim) { 388 d.d += d.d; 389 if (d.d > ds || d.d == ds && L & 1) { 390 bump_up: 391 while(*--s == '9') 392 if (s == s0) { 393 k++; 394 *s = '0'; 395 break; 396 } 397 ++*s++; 398 } 399 break; 400 } 401 d.d *= 10.; 402 if (d.d == 0.) 403 break; 404 } 405 goto ret1; 406 } 407 408 m2 = b2; 409 m5 = b5; 410 mhi = mlo = 0; 411 if (leftright) { 412 if (mode < 2) { 413 i = 414 #ifndef Sudden_Underflow 415 denorm ? be + (Bias + (P-1) - 1 + 1) : 416 #endif 417 #ifdef IBM 418 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); 419 #else 420 1 + P - bbits; 421 #endif 422 } 423 else { 424 j = ilim - 1; 425 if (m5 >= j) 426 m5 -= j; 427 else { 428 s5 += j -= m5; 429 b5 += j; 430 m5 = 0; 431 } 432 if ((i = ilim) < 0) { 433 m2 -= i; 434 i = 0; 435 } 436 } 437 b2 += i; 438 s2 += i; 439 mhi = i2b(1); 440 } 441 if (m2 > 0 && s2 > 0) { 442 i = m2 < s2 ? m2 : s2; 443 b2 -= i; 444 m2 -= i; 445 s2 -= i; 446 } 447 if (b5 > 0) { 448 if (leftright) { 449 if (m5 > 0) { 450 mhi = pow5mult(mhi, m5); 451 b1 = mult(mhi, b); 452 Bfree(b); 453 b = b1; 454 } 455 if (j = b5 - m5) 456 b = pow5mult(b, j); 457 } 458 else 459 b = pow5mult(b, b5); 460 } 461 S = i2b(1); 462 if (s5 > 0) 463 S = pow5mult(S, s5); 464 465 /* Check for special case that d is a normalized power of 2. */ 466 467 if (mode < 2) { 468 if (!word1(d) && !(word0(d) & Bndry_mask) 469 #ifndef Sudden_Underflow 470 && word0(d) & Exp_mask 471 #endif 472 ) { 473 /* The special case */ 474 b2 += Log2P; 475 s2 += Log2P; 476 spec_case = 1; 477 } 478 else 479 spec_case = 0; 480 } 481 482 /* Arrange for convenient computation of quotients: 483 * shift left if necessary so divisor has 4 leading 0 bits. 484 * 485 * Perhaps we should just compute leading 28 bits of S once 486 * and for all and pass them and a shift to quorem, so it 487 * can do shifts and ors to compute the numerator for q. 488 */ 489 #ifdef Pack_32 490 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) 491 i = 32 - i; 492 #else 493 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) 494 i = 16 - i; 495 #endif 496 if (i > 4) { 497 i -= 4; 498 b2 += i; 499 m2 += i; 500 s2 += i; 501 } 502 else if (i < 4) { 503 i += 28; 504 b2 += i; 505 m2 += i; 506 s2 += i; 507 } 508 if (b2 > 0) 509 b = lshift(b, b2); 510 if (s2 > 0) 511 S = lshift(S, s2); 512 if (k_check) { 513 if (cmp(b,S) < 0) { 514 k--; 515 b = multadd(b, 10, 0); /* we botched the k estimate */ 516 if (leftright) 517 mhi = multadd(mhi, 10, 0); 518 ilim = ilim1; 519 } 520 } 521 if (ilim <= 0 && mode > 2) { 522 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 523 /* no digits, fcvt style */ 524 no_digits: 525 k = -1 - ndigits; 526 goto ret; 527 } 528 one_digit: 529 *s++ = '1'; 530 k++; 531 goto ret; 532 } 533 if (leftright) { 534 if (m2 > 0) 535 mhi = lshift(mhi, m2); 536 537 /* Compute mlo -- check for special case 538 * that d is a normalized power of 2. 539 */ 540 541 mlo = mhi; 542 if (spec_case) { 543 mhi = Balloc(mhi->k); 544 Bcopy(mhi, mlo); 545 mhi = lshift(mhi, Log2P); 546 } 547 548 for(i = 1;;i++) { 549 dig = quorem(b,S) + '0'; 550 /* Do we yet have the shortest decimal string 551 * that will round to d? 552 */ 553 j = cmp(b, mlo); 554 delta = diff(S, mhi); 555 j1 = delta->sign ? 1 : cmp(b, delta); 556 Bfree(delta); 557 #ifndef ROUND_BIASED 558 if (j1 == 0 && !mode && !(word1(d) & 1)) { 559 if (dig == '9') 560 goto round_9_up; 561 if (j > 0) 562 dig++; 563 *s++ = dig; 564 goto ret; 565 } 566 #endif 567 if (j < 0 || j == 0 && !mode 568 #ifndef ROUND_BIASED 569 && !(word1(d) & 1) 570 #endif 571 ) { 572 if (j1 > 0) { 573 b = lshift(b, 1); 574 j1 = cmp(b, S); 575 if ((j1 > 0 || j1 == 0 && dig & 1) 576 && dig++ == '9') 577 goto round_9_up; 578 } 579 *s++ = dig; 580 goto ret; 581 } 582 if (j1 > 0) { 583 if (dig == '9') { /* possible if i == 1 */ 584 round_9_up: 585 *s++ = '9'; 586 goto roundoff; 587 } 588 *s++ = dig + 1; 589 goto ret; 590 } 591 *s++ = dig; 592 if (i == ilim) 593 break; 594 b = multadd(b, 10, 0); 595 if (mlo == mhi) 596 mlo = mhi = multadd(mhi, 10, 0); 597 else { 598 mlo = multadd(mlo, 10, 0); 599 mhi = multadd(mhi, 10, 0); 600 } 601 } 602 } 603 else 604 for(i = 1;; i++) { 605 *s++ = dig = quorem(b,S) + '0'; 606 if (i >= ilim) 607 break; 608 b = multadd(b, 10, 0); 609 } 610 611 /* Round off last digit */ 612 613 b = lshift(b, 1); 614 j = cmp(b, S); 615 if (j > 0 || j == 0 && dig & 1) { 616 roundoff: 617 while(*--s == '9') 618 if (s == s0) { 619 k++; 620 *s++ = '1'; 621 goto ret; 622 } 623 ++*s++; 624 } 625 else { 626 while(*--s == '0'); 627 s++; 628 } 629 ret: 630 Bfree(S); 631 if (mhi) { 632 if (mlo && mlo != mhi) 633 Bfree(mlo); 634 Bfree(mhi); 635 } 636 ret1: 637 Bfree(b); 638 *s = 0; 639 *decpt = k + 1; 640 if (rve) 641 *rve = s; 642 return s0; 643 } 644 645 static int 646 quorem(Bigint *b, Bigint *S) 647 { 648 int n; 649 long borrow, y; 650 unsigned long carry, q, ys; 651 unsigned long *bx, *bxe, *sx, *sxe; 652 #ifdef Pack_32 653 long z; 654 unsigned long si, zs; 655 #endif 656 657 n = S->wds; 658 #ifdef DEBUG 659 /*debug*/ if (b->wds > n) 660 /*debug*/ Bug("oversize b in quorem"); 661 #endif 662 if (b->wds < n) 663 return 0; 664 sx = S->x; 665 sxe = sx + --n; 666 bx = b->x; 667 bxe = bx + n; 668 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 669 #ifdef DEBUG 670 /*debug*/ if (q > 9) 671 /*debug*/ Bug("oversized quotient in quorem"); 672 #endif 673 if (q) { 674 borrow = 0; 675 carry = 0; 676 do { 677 #ifdef Pack_32 678 si = *sx++; 679 ys = (si & 0xffff) * q + carry; 680 zs = (si >> 16) * q + (ys >> 16); 681 carry = zs >> 16; 682 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 683 borrow = y >> 16; 684 Sign_Extend(borrow, y); 685 z = (*bx >> 16) - (zs & 0xffff) + borrow; 686 borrow = z >> 16; 687 Sign_Extend(borrow, z); 688 Storeinc(bx, z, y); 689 #else 690 ys = *sx++ * q + carry; 691 carry = ys >> 16; 692 y = *bx - (ys & 0xffff) + borrow; 693 borrow = y >> 16; 694 Sign_Extend(borrow, y); 695 *bx++ = y & 0xffff; 696 #endif 697 } 698 while(sx <= sxe); 699 if (!*bxe) { 700 bx = b->x; 701 while(--bxe > bx && !*bxe) 702 --n; 703 b->wds = n; 704 } 705 } 706 if (cmp(b, S) >= 0) { 707 q++; 708 borrow = 0; 709 carry = 0; 710 bx = b->x; 711 sx = S->x; 712 do { 713 #ifdef Pack_32 714 si = *sx++; 715 ys = (si & 0xffff) + carry; 716 zs = (si >> 16) + (ys >> 16); 717 carry = zs >> 16; 718 y = (*bx & 0xffff) - (ys & 0xffff) + borrow; 719 borrow = y >> 16; 720 Sign_Extend(borrow, y); 721 z = (*bx >> 16) - (zs & 0xffff) + borrow; 722 borrow = z >> 16; 723 Sign_Extend(borrow, z); 724 Storeinc(bx, z, y); 725 #else 726 ys = *sx++ + carry; 727 carry = ys >> 16; 728 y = *bx - (ys & 0xffff) + borrow; 729 borrow = y >> 16; 730 Sign_Extend(borrow, y); 731 *bx++ = y & 0xffff; 732 #endif 733 } 734 while(sx <= sxe); 735 bx = b->x; 736 bxe = bx + n; 737 if (!*bxe) { 738 while(--bxe > bx && !*bxe) 739 --n; 740 b->wds = n; 741 } 742 } 743 return q; 744 } 745