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