contrib/gdtoa/gdtoa.c

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