libquadmath/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                                 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  clear_trailing0:
 485                                         while(*--s == '0'){}
 486                                         ++s;
 487                                         }
 488                                 break;
 489                                 }
 490                         }
 491                 goto ret1;
 492                 }
 493
 494         m2 = b2;
 495         m5 = b5;
 496         mhi = mlo = 0;
 497         if (leftright) {
 498                 if (mode < 2) {
 499                         i = nbits - bbits;
 500                         if (be - i++ < fpi->emin)
 501                                 /* denormal */
 502                                 i = be - fpi->emin + 1;
 503                         }
 504                 else {
 505                         j = ilim - 1;
 506                         if (m5 >= j)
 507                                 m5 -= j;
 508                         else {
 509                                 s5 += j -= m5;
 510                                 b5 += j;
 511                                 m5 = 0;
 512                                 }
 513                         if ((i = ilim) < 0) {
 514                                 m2 -= i;
 515                                 i = 0;
 516                                 }
 517                         }
 518                 b2 += i;
 519                 s2 += i;
 520                 mhi = i2b(1);
 521                 }
 522         if (m2 > 0 && s2 > 0) {
 523                 i = m2 < s2 ? m2 : s2;
 524                 b2 -= i;
 525                 m2 -= i;
 526                 s2 -= i;
 527                 }
 528         if (b5 > 0) {
 529                 if (leftright) {
 530                         if (m5 > 0) {
 531                                 mhi = pow5mult(mhi, m5);
 532                                 b1 = mult(mhi, b);
 533                                 Bfree(b);
 534                                 b = b1;
 535                                 }
 536                         if ( (j = b5 - m5) !=0)
 537                                 b = pow5mult(b, j);
 538                         }
 539                 else
 540                         b = pow5mult(b, b5);
 541                 }
 542         S = i2b(1);
 543         if (s5 > 0)
 544                 S = pow5mult(S, s5);
 545
 546         /* Check for special case that d is a normalized power of 2. */
 547
 548         spec_case = 0;
 549         if (mode < 2) {
 550                 if (bbits == 1 && be0 > fpi->emin + 1) {
 551                         /* The special case */
 552                         b2++;
 553                         s2++;
 554                         spec_case = 1;
 555                         }
 556                 }
 557
 558         /* Arrange for convenient computation of quotients:
 559          * shift left if necessary so divisor has 4 leading 0 bits.
 560          *
 561          * Perhaps we should just compute leading 28 bits of S once
 562          * and for all and pass them and a shift to quorem, so it
 563          * can do shifts and ors to compute the numerator for q.
 564          */
 565         i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
 566         m2 += i;
 567         if ((b2 += i) > 0)
 568                 b = lshift(b, b2);
 569         if ((s2 += i) > 0)
 570                 S = lshift(S, s2);
 571         if (k_check) {
 572                 if (cmp(b,S) < 0) {
 573                         k--;
 574                         b = multadd(b, 10, 0);  /* we botched the k estimate */
 575                         if (leftright)
 576                                 mhi = multadd(mhi, 10, 0);
 577                         ilim = ilim1;
 578                         }
 579                 }
 580         if (ilim <= 0 && mode > 2) {
 581                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 582                         /* no digits, fcvt style */
 583  no_digits:
 584                         k = -1 - ndigits;
 585                         inex = STRTOG_Inexlo;
 586                         goto ret;
 587                         }
 588  one_digit:
 589                 inex = STRTOG_Inexhi;
 590                 *s++ = '1';
 591                 k++;
 592                 goto ret;
 593                 }
 594         if (leftright) {
 595                 if (m2 > 0)
 596                         mhi = lshift(mhi, m2);
 597
 598                 /* Compute mlo -- check for special case
 599                  * that d is a normalized power of 2.
 600                  */
 601
 602                 mlo = mhi;
 603                 if (spec_case) {
 604                         mhi = Balloc(mhi->k);
 605                         Bcopy(mhi, mlo);
 606                         mhi = lshift(mhi, 1);
 607                         }
 608
 609                 for(i = 1;;i++) {
 610                         dig = quorem(b,S) + '0';
 611                         /* Do we yet have the shortest decimal string
 612                          * that will round to d?
 613                          */
 614                         j = cmp(b, mlo);
 615                         delta = diff(S, mhi);
 616                         j1 = delta->sign ? 1 : cmp(b, delta);
 617                         Bfree(delta);
 618 #ifndef ROUND_BIASED
 619                         if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
 620                                 if (dig == '9')
 621                                         goto round_9_up;
 622                                 if (j <= 0) {
 623                                         if (b->wds > 1 || b->x[0])
 624                                                 inex = STRTOG_Inexlo;
 625                                         }
 626                                 else {
 627                                         dig++;
 628                                         inex = STRTOG_Inexhi;
 629                                         }
 630                                 *s++ = dig;
 631                                 goto ret;
 632                                 }
 633 #endif
 634                         if (j < 0 || (j == 0 && !mode
 635 #ifndef ROUND_BIASED
 636                                                         && !(bits[0] & 1)
 637 #endif
 638                                         )) {
 639                                 if (rdir && (b->wds > 1 || b->x[0])) {
 640                                         if (rdir == 2) {
 641                                                 inex = STRTOG_Inexlo;
 642                                                 goto accept;
 643                                                 }
 644                                         while (cmp(S,mhi) > 0) {
 645                                                 *s++ = dig;
 646                                                 mhi1 = multadd(mhi, 10, 0);
 647                                                 if (mlo == mhi)
 648                                                         mlo = mhi1;
 649                                                 mhi = mhi1;
 650                                                 b = multadd(b, 10, 0);
 651                                                 dig = quorem(b,S) + '0';
 652                                                 }
 653                                         if (dig++ == '9')
 654                                                 goto round_9_up;
 655                                         inex = STRTOG_Inexhi;
 656                                         goto accept;
 657                                         }
 658                                 if (j1 > 0) {
 659                                         b = lshift(b, 1);
 660                                         j1 = cmp(b, S);
 661                                         if ((j1 > 0 || (j1 == 0 && dig & 1))
 662                                         && dig++ == '9')
 663                                                 goto round_9_up;
 664                                         inex = STRTOG_Inexhi;
 665                                         }
 666                                 if (b->wds > 1 || b->x[0])
 667                                         inex = STRTOG_Inexlo;
 668  accept:
 669                                 *s++ = dig;
 670                                 goto ret;
 671                                 }
 672                         if (j1 > 0 && rdir != 2) {
 673                                 if (dig == '9') { /* possible if i == 1 */
 674  round_9_up:
 675                                         *s++ = '9';
 676                                         inex = STRTOG_Inexhi;
 677                                         goto roundoff;
 678                                         }
 679                                 inex = STRTOG_Inexhi;
 680                                 *s++ = dig + 1;
 681                                 goto ret;
 682                                 }
 683                         *s++ = dig;
 684                         if (i == ilim)
 685                                 break;
 686                         b = multadd(b, 10, 0);
 687                         if (mlo == mhi)
 688                                 mlo = mhi = multadd(mhi, 10, 0);
 689                         else {
 690                                 mlo = multadd(mlo, 10, 0);
 691                                 mhi = multadd(mhi, 10, 0);
 692                                 }
 693                         }
 694                 }
 695         else
 696                 for(i = 1;; i++) {
 697                         *s++ = dig = quorem(b,S) + '0';
 698                         if (i >= ilim)
 699                                 break;
 700                         b = multadd(b, 10, 0);
 701                         }
 702
 703         /* Round off last digit */
 704
 705         if (rdir) {
 706                 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
 707                         goto chopzeros;
 708                 goto roundoff;
 709                 }
 710         b = lshift(b, 1);
 711         j = cmp(b, S);
 712         if (j > 0 || (j == 0 && dig & 1)) {
 713  roundoff:
 714                 inex = STRTOG_Inexhi;
 715                 while(*--s == '9')
 716                         if (s == s0) {
 717                                 k++;
 718                                 *s++ = '1';
 719                                 goto ret;
 720                                 }
 721                 ++*s++;
 722                 }
 723         else {
 724  chopzeros:
 725                 if (b->wds > 1 || b->x[0])
 726                         inex = STRTOG_Inexlo;
 727                 while(*--s == '0'){}
 728                 ++s;
 729                 }
 730  ret:
 731         Bfree(S);
 732         if (mhi) {
 733                 if (mlo && mlo != mhi)
 734                         Bfree(mlo);
 735                 Bfree(mhi);
 736                 }
 737  ret1:
 738         Bfree(b);
 739         *s = 0;
 740         *decpt = k + 1;
 741         if (rve)
 742                 *rve = s;
 743         *kindp |= inex;
 744         return s0;
 745         }