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