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