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                                                 while(*--s == '0'){}
 421                                                 s++;
 422                                                 if (dval(d))
 423                                                         inex = STRTOG_Inexlo;
 424                                                 goto ret1;
 425                                                 }
 426                                         break;
 427                                         }
 428                                 }
 429 #ifndef No_leftright
 430                         }
 431 #endif
 432  fast_failed:
 433                 s = s0;
 434                 dval(d) = d2;
 435                 k = k0;
 436                 ilim = ilim0;
 437                 }
 438
 439         /* Do we have a "small" integer? */
 440
 441         if (be >= 0 && k <= Int_max) {
 442                 /* Yes. */
 443                 ds = tens[k];
 444                 if (ndigits < 0 && ilim <= 0) {
 445                         S = mhi = 0;
 446                         if (ilim < 0 || dval(d) <= 5*ds)
 447                                 goto no_digits;
 448                         goto one_digit;
 449                         }
 450                 for(i = 1;; i++, dval(d) *= 10.) {
 451                         L = dval(d) / ds;
 452                         dval(d) -= L*ds;
 453 #ifdef Check_FLT_ROUNDS
 454                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 455                         if (dval(d) < 0) {
 456                                 L--;
 457                                 dval(d) += ds;
 458                                 }
 459 #endif
 460                         *s++ = '0' + (int)L;
 461                         if (dval(d) == 0.)
 462                                 break;
 463                         if (i == ilim) {
 464                                 if (rdir) {
 465                                         if (rdir == 1)
 466                                                 goto bump_up;
 467                                         inex = STRTOG_Inexlo;
 468                                         goto ret1;
 469                                         }
 470                                 dval(d) += dval(d);
 471                                 if (dval(d) > ds || dval(d) == ds && L & 1) {
 472  bump_up:
 473                                         inex = STRTOG_Inexhi;
 474                                         while(*--s == '9')
 475                                                 if (s == s0) {
 476                                                         k++;
 477                                                         *s = '0';
 478                                                         break;
 479                                                         }
 480                                         ++*s++;
 481                                         }
 482                                 else
 483                                         inex = STRTOG_Inexlo;
 484                                 break;
 485                                 }
 486                         }
 487                 goto ret1;
 488                 }
 489
 490         m2 = b2;
 491         m5 = b5;
 492         mhi = mlo = 0;
 493         if (leftright) {
 494                 if (mode < 2) {
 495                         i = nbits - bbits;
 496                         if (be - i++ < fpi->emin)
 497                                 /* denormal */
 498                                 i = be - fpi->emin + 1;
 499                         }
 500                 else {
 501                         j = ilim - 1;
 502                         if (m5 >= j)
 503                                 m5 -= j;
 504                         else {
 505                                 s5 += j -= m5;
 506                                 b5 += j;
 507                                 m5 = 0;
 508                                 }
 509                         if ((i = ilim) < 0) {
 510                                 m2 -= i;
 511                                 i = 0;
 512                                 }
 513                         }
 514                 b2 += i;
 515                 s2 += i;
 516                 mhi = i2b(1);
 517                 }
 518         if (m2 > 0 && s2 > 0) {
 519                 i = m2 < s2 ? m2 : s2;
 520                 b2 -= i;
 521                 m2 -= i;
 522                 s2 -= i;
 523                 }
 524         if (b5 > 0) {
 525                 if (leftright) {
 526                         if (m5 > 0) {
 527                                 mhi = pow5mult(mhi, m5);
 528                                 b1 = mult(mhi, b);
 529                                 Bfree(b);
 530                                 b = b1;
 531                                 }
 532                         if ( (j = b5 - m5) !=0)
 533                                 b = pow5mult(b, j);
 534                         }
 535                 else
 536                         b = pow5mult(b, b5);
 537                 }
 538         S = i2b(1);
 539         if (s5 > 0)
 540                 S = pow5mult(S, s5);
 541
 542         /* Check for special case that d is a normalized power of 2. */
 543
 544         spec_case = 0;
 545         if (mode < 2) {
 546                 if (bbits == 1 && be0 > fpi->emin + 1) {
 547                         /* The special case */
 548                         b2++;
 549                         s2++;
 550                         spec_case = 1;
 551                         }
 552                 }
 553
 554         /* Arrange for convenient computation of quotients:
 555          * shift left if necessary so divisor has 4 leading 0 bits.
 556          *
 557          * Perhaps we should just compute leading 28 bits of S once
 558          * and for all and pass them and a shift to quorem, so it
 559          * can do shifts and ors to compute the numerator for q.
 560          */
 561 #ifdef Pack_32
 562         if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
 563                 i = 32 - i;
 564 #else
 565         if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
 566                 i = 16 - i;
 567 #endif
 568         if (i > 4) {
 569                 i -= 4;
 570                 b2 += i;
 571                 m2 += i;
 572                 s2 += i;
 573                 }
 574         else if (i < 4) {
 575                 i += 28;
 576                 b2 += i;
 577                 m2 += i;
 578                 s2 += i;
 579                 }
 580         if (b2 > 0)
 581                 b = lshift(b, b2);
 582         if (s2 > 0)
 583                 S = lshift(S, s2);
 584         if (k_check) {
 585                 if (cmp(b,S) < 0) {
 586                         k--;
 587                         b = multadd(b, 10, 0);  /* we botched the k estimate */
 588                         if (leftright)
 589                                 mhi = multadd(mhi, 10, 0);
 590                         ilim = ilim1;
 591                         }
 592                 }
 593         if (ilim <= 0 && mode > 2) {
 594                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 595                         /* no digits, fcvt style */
 596  no_digits:
 597                         k = -1 - ndigits;
 598                         inex = STRTOG_Inexlo;
 599                         goto ret;
 600                         }
 601  one_digit:
 602                 inex = STRTOG_Inexhi;
 603                 *s++ = '1';
 604                 k++;
 605                 goto ret;
 606                 }
 607         if (leftright) {
 608                 if (m2 > 0)
 609                         mhi = lshift(mhi, m2);
 610
 611                 /* Compute mlo -- check for special case
 612                  * that d is a normalized power of 2.
 613                  */
 614
 615                 mlo = mhi;
 616                 if (spec_case) {
 617                         mhi = Balloc(mhi->k);
 618                         Bcopy(mhi, mlo);
 619                         mhi = lshift(mhi, 1);
 620                         }
 621
 622                 for(i = 1;;i++) {
 623                         dig = quorem(b,S) + '0';
 624                         /* Do we yet have the shortest decimal string
 625                          * that will round to d?
 626                          */
 627                         j = cmp(b, mlo);
 628                         delta = diff(S, mhi);
 629                         j1 = delta->sign ? 1 : cmp(b, delta);
 630                         Bfree(delta);
 631 #ifndef ROUND_BIASED
 632                         if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
 633                                 if (dig == '9')
 634                                         goto round_9_up;
 635                                 if (j <= 0) {
 636                                         if (b->wds > 1 || b->x[0])
 637                                                 inex = STRTOG_Inexlo;
 638                                         }
 639                                 else {
 640                                         dig++;
 641                                         inex = STRTOG_Inexhi;
 642                                         }
 643                                 *s++ = dig;
 644                                 goto ret;
 645                                 }
 646 #endif
 647                         if (j < 0 || j == 0 && !mode
 648 #ifndef ROUND_BIASED
 649                                                         && !(bits[0] & 1)
 650 #endif
 651                                         ) {
 652                                 if (rdir && (b->wds > 1 || b->x[0])) {
 653                                         if (rdir == 2) {
 654                                                 inex = STRTOG_Inexlo;
 655                                                 goto accept;
 656                                                 }
 657                                         while (cmp(S,mhi) > 0) {
 658                                                 *s++ = dig;
 659                                                 mhi1 = multadd(mhi, 10, 0);
 660                                                 if (mlo == mhi)
 661                                                         mlo = mhi1;
 662                                                 mhi = mhi1;
 663                                                 b = multadd(b, 10, 0);
 664                                                 dig = quorem(b,S) + '0';
 665                                                 }
 666                                         if (dig++ == '9')
 667                                                 goto round_9_up;
 668                                         inex = STRTOG_Inexhi;
 669                                         goto accept;
 670                                         }
 671                                 if (j1 > 0) {
 672                                         b = lshift(b, 1);
 673                                         j1 = cmp(b, S);
 674                                         if ((j1 > 0 || j1 == 0 && dig & 1)
 675                                         && dig++ == '9')
 676                                                 goto round_9_up;
 677                                         inex = STRTOG_Inexhi;
 678                                         }
 679                                 if (b->wds > 1 || b->x[0])
 680                                         inex = STRTOG_Inexlo;
 681  accept:
 682                                 *s++ = dig;
 683                                 goto ret;
 684                                 }
 685                         if (j1 > 0 && rdir != 2) {
 686                                 if (dig == '9') { /* possible if i == 1 */
 687  round_9_up:
 688                                         *s++ = '9';
 689                                         inex = STRTOG_Inexhi;
 690                                         goto roundoff;
 691                                         }
 692                                 inex = STRTOG_Inexhi;
 693                                 *s++ = dig + 1;
 694                                 goto ret;
 695                                 }
 696                         *s++ = dig;
 697                         if (i == ilim)
 698                                 break;
 699                         b = multadd(b, 10, 0);
 700                         if (mlo == mhi)
 701                                 mlo = mhi = multadd(mhi, 10, 0);
 702                         else {
 703                                 mlo = multadd(mlo, 10, 0);
 704                                 mhi = multadd(mhi, 10, 0);
 705                                 }
 706                         }
 707                 }
 708         else
 709                 for(i = 1;; i++) {
 710                         *s++ = dig = quorem(b,S) + '0';
 711                         if (i >= ilim)
 712                                 break;
 713                         b = multadd(b, 10, 0);
 714                         }
 715
 716         /* Round off last digit */
 717
 718         if (rdir) {
 719                 if (rdir == 2 || b->wds <= 1 && !b->x[0])
 720                         goto chopzeros;
 721                 goto roundoff;
 722                 }
 723         b = lshift(b, 1);
 724         j = cmp(b, S);
 725         if (j > 0 || j == 0 && dig & 1) {
 726  roundoff:
 727                 inex = STRTOG_Inexhi;
 728                 while(*--s == '9')
 729                         if (s == s0) {
 730                                 k++;
 731                                 *s++ = '1';
 732                                 goto ret;
 733                                 }
 734                 ++*s++;
 735                 }
 736         else {
 737  chopzeros:
 738                 if (b->wds > 1 || b->x[0])
 739                         inex = STRTOG_Inexlo;
 740                 while(*--s == '0'){}
 741                 s++;
 742                 }
 743  ret:
 744         Bfree(S);
 745         if (mhi) {
 746                 if (mlo && mlo != mhi)
 747                         Bfree(mlo);
 748                 Bfree(mhi);
 749                 }
 750  ret1:
 751         Bfree(b);
 752         *s = 0;
 753         *decpt = k + 1;
 754         if (rve)
 755                 *rve = s;
 756         *kindp |= inex;
 757         return s0;
 758         }