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