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