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