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 d2, ds; 161 char *s, *s0; 162 U d, eps; 163 164#ifndef MULTIPLE_THREADS 165 if (dtoa_result) { 166 freedtoa(dtoa_result); 167 dtoa_result = 0; 168 } 169#endif 170 inex = 0; 171 kind = *kindp &= ~STRTOG_Inexact; 172 switch(kind & STRTOG_Retmask) { 173 case STRTOG_Zero: 174 goto ret_zero; 175 case STRTOG_Normal: 176 case STRTOG_Denormal: 177 break; 178 case STRTOG_Infinite: 179 *decpt = -32768; 180 return nrv_alloc("Infinity", rve, 8); 181 case STRTOG_NaN: 182 *decpt = -32768; 183 return nrv_alloc("NaN", rve, 3); 184 default: 185 return 0; 186 } 187 b = bitstob(bits, nbits = fpi->nbits, &bbits); 188 be0 = be; 189 if ( (i = trailz(b)) !=0) { 190 rshift(b, i); 191 be += i; 192 bbits -= i; 193 } 194 if (!b->wds) { 195 Bfree(b); 196 ret_zero: 197 *decpt = 1; 198 return nrv_alloc("0", rve, 1); 199 } 200 201 dval(&d) = b2d(b, &i); 202 i = be + bbits - 1; 203 word0(&d) &= Frac_mask1; 204 word0(&d) |= Exp_11; 205#ifdef IBM 206 if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0) 207 dval(&d) /= 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#ifdef IBM 232 i <<= 2; 233 i += j; 234#endif 235 ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; 236 237 /* correct assumption about exponent range */ 238 if ((j = i) < 0) 239 j = -j; 240 if ((j -= 1077) > 0) 241 ds += j * 7e-17; 242 243 k = (int)ds; 244 if (ds < 0. && ds != k) 245 k--; /* want k = floor(ds) */ 246 k_check = 1; 247#ifdef IBM 248 j = be + bbits - 1; 249 if ( (j1 = j & 3) !=0) 250 dval(&d) *= 1 << j1; 251 word0(&d) += j << Exp_shift - 2 & Exp_mask; 252#else 253 word0(&d) += (be + bbits - 1) << Exp_shift; 254#endif 255 if (k >= 0 && k <= Ten_pmax) { 256 if (dval(&d) < tens[k]) 257 k--; 258 k_check = 0; 259 } 260 j = bbits - i - 1; 261 if (j >= 0) { 262 b2 = 0; 263 s2 = j; 264 } 265 else { 266 b2 = -j; 267 s2 = 0; 268 } 269 if (k >= 0) { 270 b5 = 0; 271 s5 = k; 272 s2 += k; 273 } 274 else { 275 b2 -= k; 276 b5 = -k; 277 s5 = 0; 278 } 279 if (mode < 0 || mode > 9) 280 mode = 0; 281 try_quick = 1; 282 if (mode > 5) { 283 mode -= 4; 284 try_quick = 0; 285 } 286 leftright = 1; 287 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 288 /* silence erroneous "gcc -Wall" warning. */ 289 switch(mode) { 290 case 0: 291 case 1: 292 i = (int)(nbits * .30103) + 3; 293 ndigits = 0; 294 break; 295 case 2: 296 leftright = 0; 297 /* no break */ 298 case 4: 299 if (ndigits <= 0) 300 ndigits = 1; 301 ilim = ilim1 = i = ndigits; 302 break; 303 case 3: 304 leftright = 0; 305 /* no break */ 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 = 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 >>= 1, i++) 351 if (j & 1) { 352 ieps++; 353 ds *= bigtens[i]; 354 } 355 } 356 else { 357 ds = 1.; 358 if ( (j1 = -k) !=0) { 359 dval(&d) *= tens[j1 & 0xf]; 360 for(j = j1 >> 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 if (dval(&d)) 423 inex = STRTOG_Inexlo; 424 goto clear_trailing0; 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 clear_trailing0: 485 while(*--s == '0'){} 486 ++s; 487 } 488 break; 489 } 490 } 491 goto ret1; 492 } 493 494 m2 = b2; 495 m5 = b5; 496 mhi = mlo = 0; 497 if (leftright) { 498 if (mode < 2) { 499 i = nbits - bbits; 500 if (be - i++ < fpi->emin) 501 /* denormal */ 502 i = be - fpi->emin + 1; 503 } 504 else { 505 j = ilim - 1; 506 if (m5 >= j) 507 m5 -= j; 508 else { 509 s5 += j -= m5; 510 b5 += j; 511 m5 = 0; 512 } 513 if ((i = ilim) < 0) { 514 m2 -= i; 515 i = 0; 516 } 517 } 518 b2 += i; 519 s2 += i; 520 mhi = i2b(1); 521 } 522 if (m2 > 0 && s2 > 0) { 523 i = m2 < s2 ? m2 : s2; 524 b2 -= i; 525 m2 -= i; 526 s2 -= i; 527 } 528 if (b5 > 0) { 529 if (leftright) { 530 if (m5 > 0) { 531 mhi = pow5mult(mhi, m5); 532 b1 = mult(mhi, b); 533 Bfree(b); 534 b = b1; 535 } 536 if ( (j = b5 - m5) !=0) 537 b = pow5mult(b, j); 538 } 539 else 540 b = pow5mult(b, b5); 541 } 542 S = i2b(1); 543 if (s5 > 0) 544 S = pow5mult(S, s5); 545 546 /* Check for special case that d is a normalized power of 2. */ 547 548 spec_case = 0; 549 if (mode < 2) { 550 if (bbits == 1 && be0 > fpi->emin + 1) { 551 /* The special case */ 552 b2++; 553 s2++; 554 spec_case = 1; 555 } 556 } 557 558 /* Arrange for convenient computation of quotients: 559 * shift left if necessary so divisor has 4 leading 0 bits. 560 * 561 * Perhaps we should just compute leading 28 bits of S once 562 * and for all and pass them and a shift to quorem, so it 563 * can do shifts and ors to compute the numerator for q. 564 */ 565 i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask; 566 m2 += i; 567 if ((b2 += i) > 0) 568 b = lshift(b, b2); 569 if ((s2 += i) > 0) 570 S = lshift(S, s2); 571 if (k_check) { 572 if (cmp(b,S) < 0) { 573 k--; 574 b = multadd(b, 10, 0); /* we botched the k estimate */ 575 if (leftright) 576 mhi = multadd(mhi, 10, 0); 577 ilim = ilim1; 578 } 579 } 580 if (ilim <= 0 && mode > 2) { 581 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 582 /* no digits, fcvt style */ 583 no_digits: 584 k = -1 - ndigits; 585 inex = STRTOG_Inexlo; 586 goto ret; 587 } 588 one_digit: 589 inex = STRTOG_Inexhi; 590 *s++ = '1'; 591 k++; 592 goto ret; 593 } 594 if (leftright) { 595 if (m2 > 0) 596 mhi = lshift(mhi, m2); 597 598 /* Compute mlo -- check for special case 599 * that d is a normalized power of 2. 600 */ 601 602 mlo = mhi; 603 if (spec_case) { 604 mhi = Balloc(mhi->k); 605 Bcopy(mhi, mlo); 606 mhi = lshift(mhi, 1); 607 } 608 609 for(i = 1;;i++) { 610 dig = quorem(b,S) + '0'; 611 /* Do we yet have the shortest decimal string 612 * that will round to d? 613 */ 614 j = cmp(b, mlo); 615 delta = diff(S, mhi); 616 j1 = delta->sign ? 1 : cmp(b, delta); 617 Bfree(delta); 618#ifndef ROUND_BIASED 619 if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) { 620 if (dig == '9') 621 goto round_9_up; 622 if (j <= 0) { 623 if (b->wds > 1 || b->x[0]) 624 inex = STRTOG_Inexlo; 625 } 626 else { 627 dig++; 628 inex = STRTOG_Inexhi; 629 } 630 *s++ = dig; 631 goto ret; 632 } 633#endif 634 if (j < 0 || (j == 0 && !mode 635#ifndef ROUND_BIASED 636 && !(bits[0] & 1) 637#endif 638 )) { 639 if (rdir && (b->wds > 1 || b->x[0])) { 640 if (rdir == 2) { 641 inex = STRTOG_Inexlo; 642 goto accept; 643 } 644 while (cmp(S,mhi) > 0) { 645 *s++ = dig; 646 mhi1 = multadd(mhi, 10, 0); 647 if (mlo == mhi) 648 mlo = mhi1; 649 mhi = mhi1; 650 b = multadd(b, 10, 0); 651 dig = quorem(b,S) + '0'; 652 } 653 if (dig++ == '9') 654 goto round_9_up; 655 inex = STRTOG_Inexhi; 656 goto accept; 657 } 658 if (j1 > 0) { 659 b = lshift(b, 1); 660 j1 = cmp(b, S); 661 if ((j1 > 0 || (j1 == 0 && dig & 1)) 662 && dig++ == '9') 663 goto round_9_up; 664 inex = STRTOG_Inexhi; 665 } 666 if (b->wds > 1 || b->x[0]) 667 inex = STRTOG_Inexlo; 668 accept: 669 *s++ = dig; 670 goto ret; 671 } 672 if (j1 > 0 && rdir != 2) { 673 if (dig == '9') { /* possible if i == 1 */ 674 round_9_up: 675 *s++ = '9'; 676 inex = STRTOG_Inexhi; 677 goto roundoff; 678 } 679 inex = STRTOG_Inexhi; 680 *s++ = dig + 1; 681 goto ret; 682 } 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 (i >= ilim) 699 break; 700 b = multadd(b, 10, 0); 701 } 702 703 /* Round off last digit */ 704 705 if (rdir) { 706 if (rdir == 2 || (b->wds <= 1 && !b->x[0])) 707 goto chopzeros; 708 goto roundoff; 709 } 710 b = lshift(b, 1); 711 j = cmp(b, S); 712 if (j > 0 || (j == 0 && dig & 1)) { 713 roundoff: 714 inex = STRTOG_Inexhi; 715 while(*--s == '9') 716 if (s == s0) { 717 k++; 718 *s++ = '1'; 719 goto ret; 720 } 721 ++*s++; 722 } 723 else { 724 chopzeros: 725 if (b->wds > 1 || b->x[0]) 726 inex = STRTOG_Inexlo; 727 while(*--s == '0'){} 728 ++s; 729 } 730 ret: 731 Bfree(S); 732 if (mhi) { 733 if (mlo && mlo != mhi) 734 Bfree(mlo); 735 Bfree(mhi); 736 } 737 ret1: 738 Bfree(b); 739 *s = 0; 740 *decpt = k + 1; 741 if (rve) 742 *rve = s; 743 *kindp |= inex; 744 return s0; 745 } 746