1 /**************************************************************** 2 3 The author of this software is David M. Gay. 4 5 Copyright (C) 1998, 1999 by Lucent Technologies 6 All Rights Reserved 7 8 Permission to use, copy, modify, and distribute this software and 9 its documentation for any purpose and without fee is hereby 10 granted, provided that the above copyright notice appear in all 11 copies and that both that the copyright notice and this 12 permission notice and warranty disclaimer appear in supporting 13 documentation, and that the name of Lucent or any of its entities 14 not be used in advertising or publicity pertaining to 15 distribution of the software without specific, written prior 16 permission. 17 18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25 THIS SOFTWARE. 26 27 ****************************************************************/ 28 29 /* Please send bug reports to David M. Gay (dmg at acm dot org, 30 * with " at " changed at "@" and " dot " changed to "."). */ 31 32 #include "gdtoaimp.h" 33 34 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 35 * 36 * Inspired by "How to Print Floating-Point Numbers Accurately" by 37 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 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 * 77 dtoa 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 } 754