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