StdLib/LibC/gdtoa/gdtoa.c

   1 /* $NetBSD: gdtoa.c,v 1.1.1.1.4.1.4.1 2008/04/08 21:10:55 jdc Exp $ */
   2
   3 /****************************************************************
   4
   5 The author of this software is David M. Gay.
   6
   7 Copyright (C) 1998, 1999 by Lucent Technologies
   8 All Rights Reserved
   9
  10 Permission to use, copy, modify, and distribute this software and
  11 its documentation for any purpose and without fee is hereby
  12 granted, provided that the above copyright notice appear in all
  13 copies and that both that the copyright notice and this
  14 permission notice and warranty disclaimer appear in supporting
  15 documentation, and that the name of Lucent or any of its entities
  16 not be used in advertising or publicity pertaining to
  17 distribution of the software without specific, written prior
  18 permission.
  19
  20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  27 THIS SOFTWARE.
  28
  29 ****************************************************************/
  30
  31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  32  * with " at " changed at "@" and " dot " changed to ".").  */
  33 #include  <LibConfig.h>
  34
  35 #include "gdtoaimp.h"
  36
  37 #if defined(_MSC_VER)
  38   /* Disable warnings about conversions to narrower data types. */
  39   #pragma warning ( disable : 4244 )
  40   // Squelch bogus warnings about uninitialized variable use.
  41   #pragma warning ( disable : 4701 )
  42 #endif
  43
  44 static Bigint *
  45 bitstob(ULong *bits, int nbits, int *bbits)
  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   if (b == NULL)
  63     return NULL;
  64   be = bits + (((unsigned int)nbits - 1) >> kshift);
  65   x = x0 = b->x;
  66   do {
  67     *x++ = *bits & ALL_ON;
  68 #ifdef Pack_16
  69     *x++ = (*bits >> 16) & ALL_ON;
  70 #endif
  71     } while(++bits <= be);
  72   i = x - x0;
  73   while(!x0[--i])
  74     if (!i) {
  75       b->wds = 0;
  76       *bbits = 0;
  77       goto ret;
  78       }
  79   b->wds = i + 1;
  80   *bbits = i*ULbits + 32 - hi0bits(b->x[i]);
  81  ret:
  82   return b;
  83   }
  84
  85 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  86  *
  87  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  88  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
  89  *
  90  * Modifications:
  91  *  1. Rather than iterating, we use a simple numeric overestimate
  92  *     to determine k = floor(log10(d)).  We scale relevant
  93  *     quantities using O(log2(k)) rather than O(k) multiplications.
  94  *  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  95  *     try to generate digits strictly left to right.  Instead, we
  96  *     compute with fewer bits and propagate the carry if necessary
  97  *     when rounding the final digit up.  This is often faster.
  98  *  3. Under the assumption that input will be rounded nearest,
  99  *     mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
 100  *     That is, we allow equality in stopping tests when the
 101  *     round-nearest rule will give the same floating-point value
 102  *     as would satisfaction of the stopping test with strict
 103  *     inequality.
 104  *  4. We remove common factors of powers of 2 from relevant
 105  *     quantities.
 106  *  5. When converting floating-point integers less than 1e16,
 107  *     we use floating-point arithmetic rather than resorting
 108  *     to multiple-precision integers.
 109  *  6. When asked to produce fewer than 15 digits, we first try
 110  *     to get by with floating-point arithmetic; we resort to
 111  *     multiple-precision integer arithmetic only if we cannot
 112  *     guarantee that the floating-point calculation has given
 113  *     the correctly rounded result.  For k requested digits and
 114  *     "uniformly" distributed input, the probability is
 115  *     something like 10^(k-15) that we must resort to the Long
 116  *     calculation.
 117  */
 118
 119  char *
 120 gdtoa
 121   (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
 122 {
 123  /* Arguments ndigits and decpt are similar to the second and third
 124   arguments of ecvt and fcvt; trailing zeros are suppressed from
 125   the returned string.  If not null, *rve is set to point
 126   to the end of the return value.  If d is +-Infinity or NaN,
 127   then *decpt is set to 9999.
 128
 129   mode:
 130     0 ==> shortest string that yields d when read in
 131       and rounded to nearest.
 132     1 ==> like 0, but with Steele & White stopping rule;
 133       e.g. with IEEE P754 arithmetic , mode 0 gives
 134       1e23 whereas mode 1 gives 9.999999999999999e22.
 135     2 ==> max(1,ndigits) significant digits.  This gives a
 136       return value similar to that of ecvt, except
 137       that trailing zeros are suppressed.
 138     3 ==> through ndigits past the decimal point.  This
 139       gives a return value similar to that from fcvt,
 140       except that trailing zeros are suppressed, and
 141       ndigits can be negative.
 142     4-9 should give the same return values as 2-3, i.e.,
 143       4 <= mode <= 9 ==> same return as mode
 144       2 + (mode & 1).  These modes are mainly for
 145       debugging; often they run slower but sometimes
 146       faster than modes 2-3.
 147     4,5,8,9 ==> left-to-right digit generation.
 148     6-9 ==> don't try fast floating-point estimate
 149       (if applicable).
 150
 151     Values of mode other than 0-9 are treated as mode 0.
 152
 153     Sufficient space is allocated to the return value
 154     to hold the suppressed trailing zeros.
 155   */
 156
 157   int bbits, b2, b5, be0, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, inex;
 158   int j, jj1, k, k0, k_check, kind, leftright, m2, m5, nbits;
 159   int rdir, s2, s5, spec_case, try_quick;
 160   Long L;
 161   Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
 162   double d, d2, ds, eps;
 163   char *s, *s0;
 164
 165 #ifndef MULTIPLE_THREADS
 166   if (dtoa_result) {
 167     freedtoa(dtoa_result);
 168     dtoa_result = 0;
 169     }
 170 #endif
 171   inex = 0;
 172   if (*kindp & STRTOG_NoMemory)
 173     return NULL;
 174   kind = *kindp &= ~STRTOG_Inexact;
 175   switch(kind & STRTOG_Retmask) {
 176     case STRTOG_Zero:
 177     goto ret_zero;
 178     case STRTOG_Normal:
 179     case STRTOG_Denormal:
 180     break;
 181     case STRTOG_Infinite:
 182     *decpt = -32768;
 183     return nrv_alloc("Infinity", rve, 8);
 184     case STRTOG_NaN:
 185     *decpt = -32768;
 186     return nrv_alloc("NaN", rve, 3);
 187     default:
 188     return 0;
 189     }
 190   b = bitstob(bits, nbits = fpi->nbits, &bbits);
 191   if (b == NULL)
 192     return NULL;
 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 ( (jj1 = j & 3) !=0)
 255     dval(d) *= 1 << jj1;
 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       /*FALLTHROUGH*/
 302     case 4:
 303       if (ndigits <= 0)
 304         ndigits = 1;
 305       ilim = ilim1 = i = ndigits;
 306       break;
 307     case 3:
 308       leftright = 0;
 309       /*FALLTHROUGH*/
 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((size_t)i);
 318   if (s == NULL)
 319     return NULL;
 320
 321   if ( (rdir = fpi->rounding - 1) !=0) {
 322     if (rdir < 0)
 323       rdir = 2;
 324     if (kind & STRTOG_Neg)
 325       rdir = 3 - rdir;
 326     }
 327
 328   /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
 329
 330   if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
 331 #ifndef IMPRECISE_INEXACT
 332     && k == 0
 333 #endif
 334                 ) {
 335
 336     /* Try to get by with floating-point arithmetic. */
 337
 338     i = 0;
 339     d2 = dval(d);
 340 #ifdef IBM
 341     if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
 342       dval(d) /= 1 << j;
 343 #endif
 344     k0 = k;
 345     ilim0 = ilim;
 346     ieps = 2; /* conservative */
 347     if (k > 0) {
 348       ds = tens[k&0xf];
 349       j = (unsigned int)k >> 4;
 350       if (j & Bletch) {
 351         /* prevent overflows */
 352         j &= Bletch - 1;
 353         dval(d) /= bigtens[n_bigtens-1];
 354         ieps++;
 355         }
 356       for(; j; j /= 2, i++)
 357         if (j & 1) {
 358           ieps++;
 359           ds *= bigtens[i];
 360           }
 361       }
 362     else  {
 363       ds = 1.;
 364       if ( (jj1 = -k) !=0) {
 365         dval(d) *= tens[jj1 & 0xf];
 366         for(j = jj1 >> 4; j; j >>= 1, i++)
 367           if (j & 1) {
 368             ieps++;
 369             dval(d) *= bigtens[i];
 370             }
 371         }
 372       }
 373     if (k_check && dval(d) < 1. && ilim > 0) {
 374       if (ilim1 <= 0)
 375         goto fast_failed;
 376       ilim = ilim1;
 377       k--;
 378       dval(d) *= 10.;
 379       ieps++;
 380       }
 381     dval(eps) = ieps*dval(d) + 7.;
 382     word0(eps) -= (P-1)*Exp_msk1;
 383     if (ilim == 0) {
 384       S = mhi = 0;
 385       dval(d) -= 5.;
 386       if (dval(d) > dval(eps))
 387         goto one_digit;
 388       if (dval(d) < -dval(eps))
 389         goto no_digits;
 390       goto fast_failed;
 391       }
 392 #ifndef No_leftright
 393     if (leftright) {
 394       /* Use Steele & White method of only
 395        * generating digits needed.
 396        */
 397       dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
 398       for(i = 0;;) {
 399         L = (Long)(dval(d)/ds);
 400         dval(d) -= L*ds;
 401         *s++ = '0' + (int)L;
 402         if (dval(d) < dval(eps)) {
 403           if (dval(d))
 404             inex = STRTOG_Inexlo;
 405           goto ret1;
 406           }
 407         if (ds - dval(d) < dval(eps))
 408           goto bump_up;
 409         if (++i >= ilim)
 410           break;
 411         dval(eps) *= 10.;
 412         dval(d) *= 10.;
 413         }
 414       }
 415     else {
 416 #endif
 417       /* Generate ilim digits, then fix them up. */
 418       dval(eps) *= tens[ilim-1];
 419       for(i = 1;; i++, dval(d) *= 10.) {
 420         if ( (L = (Long)(dval(d)/ds)) !=0)
 421           dval(d) -= L*ds;
 422         *s++ = '0' + (int)L;
 423         if (i == ilim) {
 424           ds *= 0.5;
 425           if (dval(d) > ds + dval(eps))
 426             goto bump_up;
 427           else if (dval(d) < ds - dval(eps)) {
 428             while(*--s == '0'){}
 429             s++;
 430             if (dval(d))
 431               inex = STRTOG_Inexlo;
 432             goto ret1;
 433             }
 434           break;
 435           }
 436         }
 437 #ifndef No_leftright
 438       }
 439 #endif
 440  fast_failed:
 441     s = s0;
 442     dval(d) = d2;
 443     k = k0;
 444     ilim = ilim0;
 445     }
 446
 447   /* Do we have a "small" integer? */
 448
 449   if (be >= 0 && k <= Int_max) {
 450     /* Yes. */
 451     ds = tens[k];
 452     if (ndigits < 0 && ilim <= 0) {
 453       S = mhi = 0;
 454       if (ilim < 0 || dval(d) <= 5*ds)
 455         goto no_digits;
 456       goto one_digit;
 457       }
 458     for(i = 1;; i++, dval(d) *= 10.) {
 459       L = dval(d) / ds;
 460       dval(d) -= L*ds;
 461 #ifdef Check_FLT_ROUNDS
 462       /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 463       if (dval(d) < 0) {
 464         L--;
 465         dval(d) += ds;
 466         }
 467 #endif
 468       *s++ = '0' + (int)L;
 469       if (dval(d) == 0.)
 470         break;
 471       if (i == ilim) {
 472         if (rdir) {
 473           if (rdir == 1)
 474             goto bump_up;
 475           inex = STRTOG_Inexlo;
 476           goto ret1;
 477           }
 478         dval(d) += dval(d);
 479         if (dval(d) > ds || (dval(d) == ds && L & 1)) {
 480  bump_up:
 481           inex = STRTOG_Inexhi;
 482           while(*--s == '9')
 483             if (s == s0) {
 484               k++;
 485               *s = '0';
 486               break;
 487               }
 488           ++*s++;
 489           }
 490         else
 491           inex = STRTOG_Inexlo;
 492         break;
 493         }
 494       }
 495     goto ret1;
 496     }
 497
 498   m2 = b2;
 499   m5 = b5;
 500   mhi = mlo = 0;
 501   if (leftright) {
 502     if (mode < 2) {
 503       i = nbits - bbits;
 504       if (be - i++ < fpi->emin)
 505         /* denormal */
 506         i = be - fpi->emin + 1;
 507       }
 508     else {
 509       j = ilim - 1;
 510       if (m5 >= j)
 511         m5 -= j;
 512       else {
 513         s5 += j -= m5;
 514         b5 += j;
 515         m5 = 0;
 516         }
 517       if ((i = ilim) < 0) {
 518         m2 -= i;
 519         i = 0;
 520         }
 521       }
 522     b2 += i;
 523     s2 += i;
 524     mhi = i2b(1);
 525     }
 526   if (m2 > 0 && s2 > 0) {
 527     i = m2 < s2 ? m2 : s2;
 528     b2 -= i;
 529     m2 -= i;
 530     s2 -= i;
 531     }
 532   if (b5 > 0) {
 533     if (leftright) {
 534       if (m5 > 0) {
 535         mhi = pow5mult(mhi, m5);
 536         if (mhi == NULL)
 537           return NULL;
 538         b1 = mult(mhi, b);
 539         if (b1 == NULL)
 540           return NULL;
 541         Bfree(b);
 542         b = b1;
 543         }
 544       if ( (j = b5 - m5) !=0) {
 545         b = pow5mult(b, j);
 546         if (b == NULL)
 547           return NULL;
 548         }
 549       }
 550     else {
 551       b = pow5mult(b, b5);
 552       if (b == NULL)
 553         return NULL;
 554       }
 555     }
 556   S = i2b(1);
 557   if (S == NULL)
 558     return NULL;
 559   if (s5 > 0) {
 560     S = pow5mult(S, s5);
 561     if (S == NULL)
 562       return NULL;
 563     }
 564
 565   /* Check for special case that d is a normalized power of 2. */
 566
 567   spec_case = 0;
 568   if (mode < 2) {
 569     if (bbits == 1 && be0 > fpi->emin + 1) {
 570       /* The special case */
 571       b2++;
 572       s2++;
 573       spec_case = 1;
 574       }
 575     }
 576
 577   /* Arrange for convenient computation of quotients:
 578    * shift left if necessary so divisor has 4 leading 0 bits.
 579    *
 580    * Perhaps we should just compute leading 28 bits of S once
 581    * and for all and pass them and a shift to quorem, so it
 582    * can do shifts and ors to compute the numerator for q.
 583    */
 584 #ifdef Pack_32
 585   if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
 586     i = 32 - i;
 587 #else
 588   if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
 589     i = 16 - i;
 590 #endif
 591   if (i > 4) {
 592     i -= 4;
 593     b2 += i;
 594     m2 += i;
 595     s2 += i;
 596     }
 597   else if (i < 4) {
 598     i += 28;
 599     b2 += i;
 600     m2 += i;
 601     s2 += i;
 602     }
 603   if (b2 > 0)
 604     b = lshift(b, b2);
 605   if (s2 > 0)
 606     S = lshift(S, s2);
 607   if (k_check) {
 608     if (cmp(b,S) < 0) {
 609       k--;
 610       b = multadd(b, 10, 0);  /* we botched the k estimate */
 611       if (b == NULL)
 612         return NULL;
 613       if (leftright) {
 614         mhi = multadd(mhi, 10, 0);
 615         if (mhi == NULL)
 616           return NULL;
 617         }
 618       ilim = ilim1;
 619       }
 620     }
 621   if (ilim <= 0 && mode > 2) {
 622     if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 623       /* no digits, fcvt style */
 624  no_digits:
 625       k = -1 - ndigits;
 626       inex = STRTOG_Inexlo;
 627       goto ret;
 628       }
 629  one_digit:
 630     inex = STRTOG_Inexhi;
 631     *s++ = '1';
 632     k++;
 633     goto ret;
 634     }
 635   if (leftright) {
 636     if (m2 > 0) {
 637       mhi = lshift(mhi, m2);
 638       if (mhi == NULL)
 639         return NULL;
 640       }
 641
 642     /* Compute mlo -- check for special case
 643      * that d is a normalized power of 2.
 644      */
 645
 646     mlo = mhi;
 647     if (spec_case) {
 648       mhi = Balloc(mhi->k);
 649       if (mhi == NULL)
 650         return NULL;
 651       Bcopy(mhi, mlo);
 652       mhi = lshift(mhi, 1);
 653       if (mhi == NULL)
 654         return NULL;
 655       }
 656
 657     for(i = 1;;i++) {
 658       dig = quorem(b,S) + '0';
 659       /* Do we yet have the shortest decimal string
 660        * that will round to d?
 661        */
 662       j = cmp(b, mlo);
 663       delta = diff(S, mhi);
 664       if (delta == NULL)
 665         return NULL;
 666       jj1 = delta->sign ? 1 : cmp(b, delta);
 667       Bfree(delta);
 668 #ifndef ROUND_BIASED
 669       if (jj1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
 670         if (dig == '9')
 671           goto round_9_up;
 672         if (j <= 0) {
 673           if (b->wds > 1 || b->x[0])
 674             inex = STRTOG_Inexlo;
 675           }
 676         else {
 677           dig++;
 678           inex = STRTOG_Inexhi;
 679           }
 680         *s++ = dig;
 681         goto ret;
 682         }
 683 #endif
 684       if (j < 0 || (j == 0 && !mode
 685 #ifndef ROUND_BIASED
 686               && !(bits[0] & 1)
 687 #endif
 688           )) {
 689         if (rdir && (b->wds > 1 || b->x[0])) {
 690           if (rdir == 2) {
 691             inex = STRTOG_Inexlo;
 692             goto accept;
 693             }
 694           while (cmp(S,mhi) > 0) {
 695             *s++ = dig;
 696             mhi1 = multadd(mhi, 10, 0);
 697             if (mhi1 == NULL)
 698               return NULL;
 699             if (mlo == mhi)
 700               mlo = mhi1;
 701             mhi = mhi1;
 702             b = multadd(b, 10, 0);
 703             if (b == NULL)
 704               return NULL;
 705             dig = quorem(b,S) + '0';
 706             }
 707           if (dig++ == '9')
 708             goto round_9_up;
 709           inex = STRTOG_Inexhi;
 710           goto accept;
 711           }
 712         if (jj1 > 0) {
 713           b = lshift(b, 1);
 714           if (b == NULL)
 715             return NULL;
 716           jj1 = cmp(b, S);
 717           if ((jj1 > 0 || (jj1 == 0 && dig & 1))
 718           && dig++ == '9')
 719             goto round_9_up;
 720           inex = STRTOG_Inexhi;
 721           }
 722         if (b->wds > 1 || b->x[0])
 723           inex = STRTOG_Inexlo;
 724  accept:
 725         *s++ = dig;
 726         goto ret;
 727         }
 728       if (jj1 > 0 && rdir != 2) {
 729         if (dig == '9') { /* possible if i == 1 */
 730  round_9_up:
 731           *s++ = '9';
 732           inex = STRTOG_Inexhi;
 733           goto roundoff;
 734           }
 735         inex = STRTOG_Inexhi;
 736         *s++ = dig + 1;
 737         goto ret;
 738         }
 739       *s++ = dig;
 740       if (i == ilim)
 741         break;
 742       b = multadd(b, 10, 0);
 743       if (b == NULL)
 744         return NULL;
 745       if (mlo == mhi) {
 746         mlo = mhi = multadd(mhi, 10, 0);
 747         if (mlo == NULL)
 748           return NULL;
 749         }
 750       else {
 751         mlo = multadd(mlo, 10, 0);
 752         if (mlo == NULL)
 753           return NULL;
 754         mhi = multadd(mhi, 10, 0);
 755         if (mhi == NULL)
 756           return NULL;
 757         }
 758       }
 759     }
 760   else
 761     for(i = 1;; i++) {
 762       *s++ = dig = quorem(b,S) + '0';
 763       if (i >= ilim)
 764         break;
 765       b = multadd(b, 10, 0);
 766       if (b == NULL)
 767         return NULL;
 768       }
 769
 770   /* Round off last digit */
 771
 772   if (rdir) {
 773     if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
 774       goto chopzeros;
 775     goto roundoff;
 776     }
 777   b = lshift(b, 1);
 778   if (b == NULL)
 779     return NULL;
 780   j = cmp(b, S);
 781   if (j > 0 || (j == 0 && dig & 1)) {
 782  roundoff:
 783     inex = STRTOG_Inexhi;
 784     while(*--s == '9')
 785       if (s == s0) {
 786         k++;
 787         *s++ = '1';
 788         goto ret;
 789         }
 790     ++*s++;
 791     }
 792   else {
 793  chopzeros:
 794     if (b->wds > 1 || b->x[0])
 795       inex = STRTOG_Inexlo;
 796     while(*--s == '0'){}
 797     s++;
 798     }
 799  ret:
 800   Bfree(S);
 801   if (mhi) {
 802     if (mlo && mlo != mhi)
 803       Bfree(mlo);
 804     Bfree(mhi);
 805     }
 806  ret1:
 807   Bfree(b);
 808   *s = 0;
 809   *decpt = k + 1;
 810   if (rve)
 811     *rve = s;
 812   *kindp |= inex;
 813   return s0;
 814   }