1 /* $NetBSD: dtoa.c,v 1.12 2024/01/20 14:52:46 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 = (unsigned int)jj1 >> 4; j; j = (unsigned int)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 } 545 else { 546 b = pow5mult(b, b5); 547 if (b == NULL) 548 return NULL; 549 } 550 } 551 S = i2b(1); 552 if (S == NULL) 553 return NULL; 554 if (s5 > 0) { 555 S = pow5mult(S, s5); 556 if (S == NULL) 557 return NULL; 558 } 559 560 /* Check for special case that d is a normalized power of 2. */ 561 562 spec_case = 0; 563 if ((mode < 2 || leftright) 564 #ifdef Honor_FLT_ROUNDS 565 && Rounding == 1 566 #endif 567 ) { 568 if (!word1(&d) && !(word0(&d) & Bndry_mask) 569 #ifndef Sudden_Underflow 570 && word0(&d) & (Exp_mask & ~Exp_msk1) 571 #endif 572 ) { 573 /* The special case */ 574 b2 += Log2P; 575 s2 += Log2P; 576 spec_case = 1; 577 } 578 } 579 580 /* Arrange for convenient computation of quotients: 581 * shift left if necessary so divisor has 4 leading 0 bits. 582 * 583 * Perhaps we should just compute leading 28 bits of S once 584 * and for all and pass them and a shift to quorem, so it 585 * can do shifts and ors to compute the numerator for q. 586 */ 587 #ifdef Pack_32 588 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0) 589 i = 32 - i; 590 #else 591 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0) 592 i = 16 - i; 593 #endif 594 if (i > 4) { 595 i -= 4; 596 b2 += i; 597 m2 += i; 598 s2 += i; 599 } 600 else if (i < 4) { 601 i += 28; 602 b2 += i; 603 m2 += i; 604 s2 += i; 605 } 606 if (b2 > 0) { 607 b = lshift(b, b2); 608 if (b == NULL) 609 return NULL; 610 } 611 if (s2 > 0) { 612 S = lshift(S, s2); 613 if (S == NULL) 614 return NULL; 615 } 616 if (k_check) { 617 if (cmp(b,S) < 0) { 618 k--; 619 b = multadd(b, 10, 0); /* we botched the k estimate */ 620 if (b == NULL) 621 return NULL; 622 if (leftright) { 623 mhi = multadd(mhi, 10, 0); 624 if (mhi == NULL) 625 return NULL; 626 } 627 ilim = ilim1; 628 } 629 } 630 if (ilim <= 0 && (mode == 3 || mode == 5)) { 631 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 632 /* no digits, fcvt style */ 633 no_digits: 634 k = -1 - ndigits; 635 goto ret; 636 } 637 one_digit: 638 *s++ = '1'; 639 k++; 640 goto ret; 641 } 642 if (leftright) { 643 if (m2 > 0) { 644 mhi = lshift(mhi, m2); 645 if (mhi == NULL) 646 return NULL; 647 } 648 649 /* Compute mlo -- check for special case 650 * that d is a normalized power of 2. 651 */ 652 653 mlo = mhi; 654 if (spec_case) { 655 mhi = Balloc(mhi->k); 656 if (mhi == NULL) 657 return NULL; 658 Bcopy(mhi, mlo); 659 mhi = lshift(mhi, Log2P); 660 if (mhi == NULL) 661 return NULL; 662 } 663 664 for(i = 1;;i++) { 665 dig = quorem(b,S) + '0'; 666 /* Do we yet have the shortest decimal string 667 * that will round to d? 668 */ 669 j = cmp(b, mlo); 670 delta = diff(S, mhi); 671 if (delta == NULL) 672 return NULL; 673 jj1 = delta->sign ? 1 : cmp(b, delta); 674 Bfree(delta); 675 #ifndef ROUND_BIASED 676 if (jj1 == 0 && mode != 1 && !(word1(&d) & 1) 677 #ifdef Honor_FLT_ROUNDS 678 && Rounding >= 1 679 #endif 680 ) { 681 if (dig == '9') 682 goto round_9_up; 683 if (j > 0) 684 dig++; 685 #ifdef SET_INEXACT 686 else if (!b->x[0] && b->wds <= 1) 687 inexact = 0; 688 #endif 689 *s++ = dig; 690 goto ret; 691 } 692 #endif 693 if (j < 0 || (j == 0 && mode != 1 694 #ifndef ROUND_BIASED 695 && !(word1(&d) & 1) 696 #endif 697 )) { 698 if (!b->x[0] && b->wds <= 1) { 699 #ifdef SET_INEXACT 700 inexact = 0; 701 #endif 702 goto accept_dig; 703 } 704 #ifdef Honor_FLT_ROUNDS 705 if (mode > 1) 706 switch(Rounding) { 707 case 0: goto accept_dig; 708 case 2: goto keep_dig; 709 } 710 #endif /*Honor_FLT_ROUNDS*/ 711 if (jj1 > 0) { 712 b = lshift(b, 1); 713 if (b == NULL) 714 return NULL; 715 jj1 = cmp(b, S); 716 #ifdef ROUND_BIASED 717 if (jj1 >= 0 /*)*/ 718 #else 719 if ((jj1 > 0 || (jj1 == 0 && dig & 1)) 720 #endif 721 && dig++ == '9') 722 goto round_9_up; 723 } 724 accept_dig: 725 *s++ = dig; 726 goto ret; 727 } 728 if (jj1 > 0) { 729 #ifdef Honor_FLT_ROUNDS 730 if (!Rounding) 731 goto accept_dig; 732 #endif 733 if (dig == '9') { /* possible if i == 1 */ 734 round_9_up: 735 *s++ = '9'; 736 goto roundoff; 737 } 738 *s++ = dig + 1; 739 goto ret; 740 } 741 #ifdef Honor_FLT_ROUNDS 742 keep_dig: 743 #endif 744 *s++ = dig; 745 if (i == ilim) 746 break; 747 b = multadd(b, 10, 0); 748 if (b == NULL) 749 return NULL; 750 if (mlo == mhi) { 751 mlo = mhi = multadd(mhi, 10, 0); 752 if (mlo == NULL) 753 return NULL; 754 } 755 else { 756 mlo = multadd(mlo, 10, 0); 757 if (mlo == NULL) 758 return NULL; 759 mhi = multadd(mhi, 10, 0); 760 if (mhi == NULL) 761 return NULL; 762 } 763 } 764 } 765 else 766 for(i = 1;; i++) { 767 *s++ = dig = quorem(b,S) + '0'; 768 if (!b->x[0] && b->wds <= 1) { 769 #ifdef SET_INEXACT 770 inexact = 0; 771 #endif 772 goto ret; 773 } 774 if (i >= ilim) 775 break; 776 b = multadd(b, 10, 0); 777 if (b == NULL) 778 return NULL; 779 } 780 781 /* Round off last digit */ 782 783 #ifdef Honor_FLT_ROUNDS 784 switch(Rounding) { 785 case 0: goto trimzeros; 786 case 2: goto roundoff; 787 } 788 #endif 789 b = lshift(b, 1); 790 if (b == NULL) 791 return NULL; 792 j = cmp(b, S); 793 #ifdef ROUND_BIASED 794 if (j >= 0) 795 #else 796 if (j > 0 || (j == 0 && dig & 1)) 797 #endif 798 { 799 roundoff: 800 while(*--s == '9') 801 if (s == s0) { 802 k++; 803 *s++ = '1'; 804 goto ret; 805 } 806 ++*s++; 807 } 808 else { 809 #ifdef Honor_FLT_ROUNDS 810 trimzeros: 811 #endif 812 while(*--s == '0'); 813 s++; 814 } 815 ret: 816 Bfree(S); 817 if (mhi) { 818 if (mlo && mlo != mhi) 819 Bfree(mlo); 820 Bfree(mhi); 821 } 822 ret1: 823 #ifdef SET_INEXACT 824 if (inexact) { 825 if (!oldinexact) { 826 word0(&d) = Exp_1 + (70 << Exp_shift); 827 word1(&d) = 0; 828 dval(&d) += 1.; 829 } 830 } 831 else if (!oldinexact) 832 clear_inexact(); 833 #endif 834 Bfree(b); 835 if (s == s0) { /* don't return empty string */ 836 *s++ = '0'; 837 k = 0; 838 } 839 *s = 0; 840 *decpt = k + 1; 841 if (rve) 842 *rve = s; 843 return s0; 844 } 845
Indexes created Sun Sep 21 20:09:37 GMT 2025