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