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