gdtoa.c revision 1.7
1/* $NetBSD: gdtoa.c,v 1.7 2019/08/01 02:27:43 riastradh 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 9999. 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 ((s2 += i) > 0) 607 S = lshift(S, s2); 608 if (k_check) { 609 if (cmp(b,S) < 0) { 610 k--; 611 b = multadd(b, 10, 0); /* we botched the k estimate */ 612 if (b == NULL) 613 return NULL; 614 if (leftright) { 615 mhi = multadd(mhi, 10, 0); 616 if (mhi == NULL) 617 return NULL; 618 } 619 ilim = ilim1; 620 } 621 } 622 if (ilim <= 0 && mode > 2) { 623 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { 624 /* no digits, fcvt style */ 625 no_digits: 626 k = -1 - ndigits; 627 inex = STRTOG_Inexlo; 628 goto ret; 629 } 630 one_digit: 631 inex = STRTOG_Inexhi; 632 *s++ = '1'; 633 k++; 634 goto ret; 635 } 636 if (leftright) { 637 if (m2 > 0) { 638 mhi = lshift(mhi, m2); 639 if (mhi == NULL) 640 return NULL; 641 } 642 643 /* Compute mlo -- check for special case 644 * that d is a normalized power of 2. 645 */ 646 647 mlo = mhi; 648 if (spec_case) { 649 mhi = Balloc(mhi->k); 650 if (mhi == NULL) 651 return NULL; 652 Bcopy(mhi, mlo); 653 mhi = lshift(mhi, 1); 654 if (mhi == NULL) 655 return NULL; 656 } 657 658 for(i = 1;;i++) { 659 dig = quorem(b,S) + '0'; 660 /* Do we yet have the shortest decimal string 661 * that will round to d? 662 */ 663 j = cmp(b, mlo); 664 delta = diff(S, mhi); 665 if (delta == NULL) 666 return NULL; 667 jj1 = delta->sign ? 1 : cmp(b, delta); 668 Bfree(delta); 669#ifndef ROUND_BIASED 670 if (jj1 == 0 && !mode && !(bits[0] & 1) && !rdir) { 671 if (dig == '9') 672 goto round_9_up; 673 if (j <= 0) { 674 if (b->wds > 1 || b->x[0]) 675 inex = STRTOG_Inexlo; 676 } 677 else { 678 dig++; 679 inex = STRTOG_Inexhi; 680 } 681 *s++ = dig; 682 goto ret; 683 } 684#endif 685 if (j < 0 || (j == 0 && !mode 686#ifndef ROUND_BIASED 687 && !(bits[0] & 1) 688#endif 689 )) { 690 if (rdir && (b->wds > 1 || b->x[0])) { 691 if (rdir == 2) { 692 inex = STRTOG_Inexlo; 693 goto accept; 694 } 695 while (cmp(S,mhi) > 0) { 696 *s++ = dig; 697 mhi1 = multadd(mhi, 10, 0); 698 if (mhi1 == NULL) 699 return NULL; 700 if (mlo == mhi) 701 mlo = mhi1; 702 mhi = mhi1; 703 b = multadd(b, 10, 0); 704 if (b == NULL) 705 return NULL; 706 dig = quorem(b,S) + '0'; 707 } 708 if (dig++ == '9') 709 goto round_9_up; 710 inex = STRTOG_Inexhi; 711 goto accept; 712 } 713 if (jj1 > 0) { 714 b = lshift(b, 1); 715 if (b == NULL) 716 return NULL; 717 jj1 = cmp(b, S); 718#ifdef ROUND_BIASED 719 if (jj1 >= 0 /*)*/ 720#else 721 if ((jj1 > 0 || (jj1 == 0 && dig & 1)) 722#endif 723 && dig++ == '9') 724 goto round_9_up; 725 inex = STRTOG_Inexhi; 726 } 727 if (b->wds > 1 || b->x[0]) 728 inex = STRTOG_Inexlo; 729 accept: 730 *s++ = dig; 731 goto ret; 732 } 733 if (jj1 > 0 && rdir != 2) { 734 if (dig == '9') { /* possible if i == 1 */ 735 round_9_up: 736 *s++ = '9'; 737 inex = STRTOG_Inexhi; 738 goto roundoff; 739 } 740 inex = STRTOG_Inexhi; 741 *s++ = dig + 1; 742 goto ret; 743 } 744 *s++ = dig; 745 if (i == ilim) 746 break; 747 b = multadd(b, 10, 0); 748 if (b == NULL) 749 return NULL; 750 if (mlo == mhi) { 751 mlo = mhi = multadd(mhi, 10, 0); 752 if (mlo == NULL) 753 return NULL; 754 } 755 else { 756 mlo = multadd(mlo, 10, 0); 757 if (mlo == NULL) 758 return NULL; 759 mhi = multadd(mhi, 10, 0); 760 if (mhi == NULL) 761 return NULL; 762 } 763 } 764 } 765 else 766 for(i = 1;; i++) { 767 *s++ = dig = quorem(b,S) + '0'; 768 if (i >= ilim) 769 break; 770 b = multadd(b, 10, 0); 771 if (b == NULL) 772 return NULL; 773 } 774 775 /* Round off last digit */ 776 777 if (rdir) { 778 if (rdir == 2 || (b->wds <= 1 && !b->x[0])) 779 goto chopzeros; 780 goto roundoff; 781 } 782 b = lshift(b, 1); 783 if (b == NULL) 784 return NULL; 785 j = cmp(b, S); 786#ifdef ROUND_BIASED 787 if (j >= 0) 788#else 789 if (j > 0 || (j == 0 && dig & 1)) 790#endif 791 { 792 roundoff: 793 inex = STRTOG_Inexhi; 794 while(*--s == '9') 795 if (s == s0) { 796 k++; 797 *s++ = '1'; 798 goto ret; 799 } 800 ++*s++; 801 } 802 else { 803 chopzeros: 804 if (b->wds > 1 || b->x[0]) 805 inex = STRTOG_Inexlo; 806 while(*--s == '0'){} 807 ++s; 808 } 809 ret: 810 Bfree(S); 811 if (mhi) { 812 if (mlo && mlo != mhi) 813 Bfree(mlo); 814 Bfree(mhi); 815 } 816 ret1: 817 Bfree(b); 818 *s = 0; 819 *decpt = k + 1; 820 if (rve) 821 *rve = s; 822 *kindp |= inex; 823 return s0; 824 } 825