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