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