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