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