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 static Bigint * 35 #ifdef KR_headers 36 bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits; 37 #else 38 bitstob(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 * 112 gdtoa 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 if (dval(d)) 421 inex = STRTOG_Inexlo; 422 goto clear_trailing0; 423 } 424 break; 425 } 426 } 427 #ifndef No_leftright 428 } 429 #endif 430 fast_failed: 431 s = s0; 432 dval(d) = d2; 433 k = k0; 434 ilim = ilim0; 435 } 436 437 /* Do we have a "small" integer? */ 438 439 if (be >= 0 && k <= Int_max) { 440 /* Yes. */ 441 ds = tens[k]; 442 if (ndigits < 0 && ilim <= 0) { 443 S = mhi = 0; 444 if (ilim < 0 || dval(d) <= 5*ds) 445 goto no_digits; 446 goto one_digit; 447 } 448 for(i = 1;; i++, dval(d) *= 10.) { 449 L = dval(d) / ds; 450 dval(d) -= L*ds; 451 #ifdef Check_FLT_ROUNDS 452 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ 453 if (dval(d) < 0) { 454 L--; 455 dval(d) += ds; 456 } 457 #endif 458 *s++ = '0' + (int)L; 459 if (dval(d) == 0.) 460 break; 461 if (i == ilim) { 462 if (rdir) { 463 if (rdir == 1) 464 goto bump_up; 465 inex = STRTOG_Inexlo; 466 goto ret1; 467 } 468 dval(d) += dval(d); 469 if (dval(d) > ds || dval(d) == ds && L & 1) { 470 bump_up: 471 inex = STRTOG_Inexhi; 472 while(*--s == '9') 473 if (s == s0) { 474 k++; 475 *s = '0'; 476 break; 477 } 478 ++*s++; 479 } 480 else { 481 inex = STRTOG_Inexlo; 482 clear_trailing0: 483 while(*--s == '0'){} 484 ++s; 485 } 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 j1 = delta->sign ? 1 : cmp(b, delta); 632 Bfree(delta); 633 #ifndef ROUND_BIASED 634 if (j1 == 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 (j1 > 0) { 674 b = lshift(b, 1); 675 j1 = cmp(b, S); 676 if ((j1 > 0 || j1 == 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 (j1 > 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