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