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