contrib/gdtoa/dtoa.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 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  35  *
  36  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  37  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
  38  *
  39  * Modifications:
  40  *      1. Rather than iterating, we use a simple numeric overestimate
  41  *         to determine k = floor(log10(d)).  We scale relevant
  42  *         quantities using O(log2(k)) rather than O(k) multiplications.
  43  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  44  *         try to generate digits strictly left to right.  Instead, we
  45  *         compute with fewer bits and propagate the carry if necessary
  46  *         when rounding the final digit up.  This is often faster.
  47  *      3. Under the assumption that input will be rounded nearest,
  48  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  49  *         That is, we allow equality in stopping tests when the
  50  *         round-nearest rule will give the same floating-point value
  51  *         as would satisfaction of the stopping test with strict
  52  *         inequality.
  53  *      4. We remove common factors of powers of 2 from relevant
  54  *         quantities.
  55  *      5. When converting floating-point integers less than 1e16,
  56  *         we use floating-point arithmetic rather than resorting
  57  *         to multiple-precision integers.
  58  *      6. When asked to produce fewer than 15 digits, we first try
  59  *         to get by with floating-point arithmetic; we resort to
  60  *         multiple-precision integer arithmetic only if we cannot
  61  *         guarantee that the floating-point calculation has given
  62  *         the correctly rounded result.  For k requested digits and
  63  *         "uniformly" distributed input, the probability is
  64  *         something like 10^(k-15) that we must resort to the Long
  65  *         calculation.
  66  */
  67
  68 #ifdef Honor_FLT_ROUNDS
  69 #undef Check_FLT_ROUNDS
  70 #define Check_FLT_ROUNDS
  71 #else
  72 #define Rounding Flt_Rounds
  73 #endif
  74
  75  char *
  76 dtoa
  77 #ifdef KR_headers
  78         (d, mode, ndigits, decpt, sign, rve)
  79         double d; int mode, ndigits, *decpt, *sign; char **rve;
  80 #else
  81         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  82 #endif
  83 {
  84  /*     Arguments ndigits, decpt, sign are similar to those
  85         of ecvt and fcvt; trailing zeros are suppressed from
  86         the returned string.  If not null, *rve is set to point
  87         to the end of the return value.  If d is +-Infinity or NaN,
  88         then *decpt is set to 9999.
  89
  90         mode:
  91                 0 ==> shortest string that yields d when read in
  92                         and rounded to nearest.
  93                 1 ==> like 0, but with Steele & White stopping rule;
  94                         e.g. with IEEE P754 arithmetic , mode 0 gives
  95                         1e23 whereas mode 1 gives 9.999999999999999e22.
  96                 2 ==> max(1,ndigits) significant digits.  This gives a
  97                         return value similar to that of ecvt, except
  98                         that trailing zeros are suppressed.
  99                 3 ==> through ndigits past the decimal point.  This
 100                         gives a return value similar to that from fcvt,
 101                         except that trailing zeros are suppressed, and
 102                         ndigits can be negative.
 103                 4,5 ==> similar to 2 and 3, respectively, but (in
 104                         round-nearest mode) with the tests of mode 0 to
 105                         possibly return a shorter string that rounds to d.
 106                         With IEEE arithmetic and compilation with
 107                         -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
 108                         as modes 2 and 3 when FLT_ROUNDS != 1.
 109                 6-9 ==> Debugging modes similar to mode - 4:  don't try
 110                         fast floating-point estimate (if applicable).
 111
 112                 Values of mode other than 0-9 are treated as mode 0.
 113
 114                 Sufficient space is allocated to the return value
 115                 to hold the suppressed trailing zeros.
 116         */
 117
 118         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
 119                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
 120                 spec_case, try_quick;
 121         Long L;
 122 #ifndef Sudden_Underflow
 123         int denorm;
 124         ULong x;
 125 #endif
 126         Bigint *b, *b1, *delta, *mlo, *mhi, *S;
 127         double d2, ds, eps;
 128         char *s, *s0;
 129 #ifdef SET_INEXACT
 130         int inexact, oldinexact;
 131 #endif
 132 #ifdef Honor_FLT_ROUNDS /*{*/
 133         int Rounding;
 134 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
 135         Rounding = Flt_Rounds;
 136 #else /*}{*/
 137         Rounding = 1;
 138         switch(fegetround()) {
 139           case FE_TOWARDZERO:   Rounding = 0; break;
 140           case FE_UPWARD:       Rounding = 2; break;
 141           case FE_DOWNWARD:     Rounding = 3;
 142           }
 143 #endif /*}}*/
 144 #endif /*}*/
 145
 146 #ifndef MULTIPLE_THREADS
 147         if (dtoa_result) {
 148                 freedtoa(dtoa_result);
 149                 dtoa_result = 0;
 150                 }
 151 #endif
 152
 153         if (word0(d) & Sign_bit) {
 154                 /* set sign for everything, including 0's and NaNs */
 155                 *sign = 1;
 156                 word0(d) &= ~Sign_bit;  /* clear sign bit */
 157                 }
 158         else
 159                 *sign = 0;
 160
 161 #if defined(IEEE_Arith) + defined(VAX)
 162 #ifdef IEEE_Arith
 163         if ((word0(d) & Exp_mask) == Exp_mask)
 164 #else
 165         if (word0(d)  == 0x8000)
 166 #endif
 167                 {
 168                 /* Infinity or NaN */
 169                 *decpt = 9999;
 170 #ifdef IEEE_Arith
 171                 if (!word1(d) && !(word0(d) & 0xfffff))
 172                         return nrv_alloc("Infinity", rve, 8);
 173 #endif
 174                 return nrv_alloc("NaN", rve, 3);
 175                 }
 176 #endif
 177 #ifdef IBM
 178         dval(d) += 0; /* normalize */
 179 #endif
 180         if (!dval(d)) {
 181                 *decpt = 1;
 182                 return nrv_alloc("0", rve, 1);
 183                 }
 184
 185 #ifdef SET_INEXACT
 186         try_quick = oldinexact = get_inexact();
 187         inexact = 1;
 188 #endif
 189 #ifdef Honor_FLT_ROUNDS
 190         if (Rounding >= 2) {
 191                 if (*sign)
 192                         Rounding = Rounding == 2 ? 0 : 2;
 193                 else
 194                         if (Rounding != 2)
 195                                 Rounding = 0;
 196                 }
 197 #endif
 198
 199         b = d2b(dval(d), &be, &bbits);
 200 #ifdef Sudden_Underflow
 201         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
 202 #else
 203         if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
 204 #endif
 205                 dval(d2) = dval(d);
 206                 word0(d2) &= Frac_mask1;
 207                 word0(d2) |= Exp_11;
 208 #ifdef IBM
 209                 if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
 210                         dval(d2) /= 1 << j;
 211 #endif
 212
 213                 /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
 214                  * log10(x)      =  log(x) / log(10)
 215                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 216                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
 217                  *
 218                  * This suggests computing an approximation k to log10(d) by
 219                  *
 220                  * k = (i - Bias)*0.301029995663981
 221                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 222                  *
 223                  * We want k to be too large rather than too small.
 224                  * The error in the first-order Taylor series approximation
 225                  * is in our favor, so we just round up the constant enough
 226                  * to compensate for any error in the multiplication of
 227                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 228                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 229                  * adding 1e-13 to the constant term more than suffices.
 230                  * Hence we adjust the constant term to 0.1760912590558.
 231                  * (We could get a more accurate k by invoking log10,
 232                  *  but this is probably not worthwhile.)
 233                  */
 234
 235                 i -= Bias;
 236 #ifdef IBM
 237                 i <<= 2;
 238                 i += j;
 239 #endif
 240 #ifndef Sudden_Underflow
 241                 denorm = 0;
 242                 }
 243         else {
 244                 /* d is denormalized */
 245
 246                 i = bbits + be + (Bias + (P-1) - 1);
 247                 x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
 248                             : word1(d) << 32 - i;
 249                 dval(d2) = x;
 250                 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
 251                 i -= (Bias + (P-1) - 1) + 1;
 252                 denorm = 1;
 253                 }
 254 #endif
 255         ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 256         k = (int)ds;
 257         if (ds < 0. && ds != k)
 258                 k--;    /* want k = floor(ds) */
 259         k_check = 1;
 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
 287 #ifndef SET_INEXACT
 288 #ifdef Check_FLT_ROUNDS
 289         try_quick = Rounding == 1;
 290 #else
 291         try_quick = 1;
 292 #endif
 293 #endif /*SET_INEXACT*/
 294
 295         if (mode > 5) {
 296                 mode -= 4;
 297                 try_quick = 0;
 298                 }
 299         leftright = 1;
 300         switch(mode) {
 301                 case 0:
 302                 case 1:
 303                         ilim = ilim1 = -1;
 304                         i = 18;
 305                         ndigits = 0;
 306                         break;
 307                 case 2:
 308                         leftright = 0;
 309                         /* no break */
 310                 case 4:
 311                         if (ndigits <= 0)
 312                                 ndigits = 1;
 313                         ilim = ilim1 = i = ndigits;
 314                         break;
 315                 case 3:
 316                         leftright = 0;
 317                         /* no break */
 318                 case 5:
 319                         i = ndigits + k + 1;
 320                         ilim = i;
 321                         ilim1 = i - 1;
 322                         if (i <= 0)
 323                                 i = 1;
 324                 }
 325         s = s0 = rv_alloc(i);
 326
 327 #ifdef Honor_FLT_ROUNDS
 328         if (mode > 1 && Rounding != 1)
 329                 leftright = 0;
 330 #endif
 331
 332         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
 333
 334                 /* Try to get by with floating-point arithmetic. */
 335
 336                 i = 0;
 337                 dval(d2) = dval(d);
 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                         dval(d) /= ds;
 356                         }
 357                 else if (( j1 = -k )!=0) {
 358                         dval(d) *= tens[j1 & 0xf];
 359                         for(j = j1 >> 4; j; j >>= 1, i++)
 360                                 if (j & 1) {
 361                                         ieps++;
 362                                         dval(d) *= bigtens[i];
 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) = 0.5/tens[ilim-1] - dval(eps);
 390                         for(i = 0;;) {
 391                                 L = dval(d);
 392                                 dval(d) -= L;
 393                                 *s++ = '0' + (int)L;
 394                                 if (dval(d) < dval(eps))
 395                                         goto ret1;
 396                                 if (1. - dval(d) < dval(eps))
 397                                         goto bump_up;
 398                                 if (++i >= ilim)
 399                                         break;
 400                                 dval(eps) *= 10.;
 401                                 dval(d) *= 10.;
 402                                 }
 403                         }
 404                 else {
 405 #endif
 406                         /* Generate ilim digits, then fix them up. */
 407                         dval(eps) *= tens[ilim-1];
 408                         for(i = 1;; i++, dval(d) *= 10.) {
 409                                 L = (Long)(dval(d));
 410                                 if (!(dval(d) -= L))
 411                                         ilim = i;
 412                                 *s++ = '0' + (int)L;
 413                                 if (i == ilim) {
 414                                         if (dval(d) > 0.5 + dval(eps))
 415                                                 goto bump_up;
 416                                         else if (dval(d) < 0.5 - dval(eps)) {
 417                                                 while(*--s == '0');
 418                                                 s++;
 419                                                 goto ret1;
 420                                                 }
 421                                         break;
 422                                         }
 423                                 }
 424 #ifndef No_leftright
 425                         }
 426 #endif
 427  fast_failed:
 428                 s = s0;
 429                 dval(d) = dval(d2);
 430                 k = k0;
 431                 ilim = ilim0;
 432                 }
 433
 434         /* Do we have a "small" integer? */
 435
 436         if (be >= 0 && k <= Int_max) {
 437                 /* Yes. */
 438                 ds = tens[k];
 439                 if (ndigits < 0 && ilim <= 0) {
 440                         S = mhi = 0;
 441                         if (ilim < 0 || dval(d) <= 5*ds)
 442                                 goto no_digits;
 443                         goto one_digit;
 444                         }
 445                 for(i = 1;; i++, dval(d) *= 10.) {
 446                         L = (Long)(dval(d) / ds);
 447                         dval(d) -= L*ds;
 448 #ifdef Check_FLT_ROUNDS
 449                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 450                         if (dval(d) < 0) {
 451                                 L--;
 452                                 dval(d) += ds;
 453                                 }
 454 #endif
 455                         *s++ = '0' + (int)L;
 456                         if (!dval(d)) {
 457 #ifdef SET_INEXACT
 458                                 inexact = 0;
 459 #endif
 460                                 break;
 461                                 }
 462                         if (i == ilim) {
 463 #ifdef Honor_FLT_ROUNDS
 464                                 if (mode > 1)
 465                                 switch(Rounding) {
 466                                   case 0: goto ret1;
 467                                   case 2: goto bump_up;
 468                                   }
 469 #endif
 470                                 dval(d) += dval(d);
 471                                 if (dval(d) > ds || dval(d) == ds && L & 1) {
 472  bump_up:
 473                                         while(*--s == '9')
 474                                                 if (s == s0) {
 475                                                         k++;
 476                                                         *s = '0';
 477                                                         break;
 478                                                         }
 479                                         ++*s++;
 480                                         }
 481                                 break;
 482                                 }
 483                         }
 484                 goto ret1;
 485                 }
 486
 487         m2 = b2;
 488         m5 = b5;
 489         mhi = mlo = 0;
 490         if (leftright) {
 491                 i =
 492 #ifndef Sudden_Underflow
 493                         denorm ? be + (Bias + (P-1) - 1 + 1) :
 494 #endif
 495 #ifdef IBM
 496                         1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
 497 #else
 498                         1 + P - bbits;
 499 #endif
 500                 b2 += i;
 501                 s2 += i;
 502                 mhi = i2b(1);
 503                 }
 504         if (m2 > 0 && s2 > 0) {
 505                 i = m2 < s2 ? m2 : s2;
 506                 b2 -= i;
 507                 m2 -= i;
 508                 s2 -= i;
 509                 }
 510         if (b5 > 0) {
 511                 if (leftright) {
 512                         if (m5 > 0) {
 513                                 mhi = pow5mult(mhi, m5);
 514                                 b1 = mult(mhi, b);
 515                                 Bfree(b);
 516                                 b = b1;
 517                                 }
 518                         if (( j = b5 - m5 )!=0)
 519                                 b = pow5mult(b, j);
 520                         }
 521                 else
 522                         b = pow5mult(b, b5);
 523                 }
 524         S = i2b(1);
 525         if (s5 > 0)
 526                 S = pow5mult(S, s5);
 527
 528         /* Check for special case that d is a normalized power of 2. */
 529
 530         spec_case = 0;
 531         if ((mode < 2 || leftright)
 532 #ifdef Honor_FLT_ROUNDS
 533                         && Rounding == 1
 534 #endif
 535                                 ) {
 536                 if (!word1(d) && !(word0(d) & Bndry_mask)
 537 #ifndef Sudden_Underflow
 538                  && word0(d) & (Exp_mask & ~Exp_msk1)
 539 #endif
 540                                 ) {
 541                         /* The special case */
 542                         b2 += Log2P;
 543                         s2 += Log2P;
 544                         spec_case = 1;
 545                         }
 546                 }
 547
 548         /* Arrange for convenient computation of quotients:
 549          * shift left if necessary so divisor has 4 leading 0 bits.
 550          *
 551          * Perhaps we should just compute leading 28 bits of S once
 552          * and for all and pass them and a shift to quorem, so it
 553          * can do shifts and ors to compute the numerator for q.
 554          */
 555 #ifdef Pack_32
 556         if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
 557                 i = 32 - i;
 558 #else
 559         if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
 560                 i = 16 - i;
 561 #endif
 562         if (i > 4) {
 563                 i -= 4;
 564                 b2 += i;
 565                 m2 += i;
 566                 s2 += i;
 567                 }
 568         else if (i < 4) {
 569                 i += 28;
 570                 b2 += i;
 571                 m2 += i;
 572                 s2 += i;
 573                 }
 574         if (b2 > 0)
 575                 b = lshift(b, b2);
 576         if (s2 > 0)
 577                 S = lshift(S, s2);
 578         if (k_check) {
 579                 if (cmp(b,S) < 0) {
 580                         k--;
 581                         b = multadd(b, 10, 0);  /* we botched the k estimate */
 582                         if (leftright)
 583                                 mhi = multadd(mhi, 10, 0);
 584                         ilim = ilim1;
 585                         }
 586                 }
 587         if (ilim <= 0 && (mode == 3 || mode == 5)) {
 588                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 589                         /* no digits, fcvt style */
 590  no_digits:
 591                         k = -1 - ndigits;
 592                         goto ret;
 593                         }
 594  one_digit:
 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, Log2P);
 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 != 1 && !(word1(d) & 1)
 625 #ifdef Honor_FLT_ROUNDS
 626                                 && Rounding >= 1
 627 #endif
 628                                                                    ) {
 629                                 if (dig == '9')
 630                                         goto round_9_up;
 631                                 if (j > 0)
 632                                         dig++;
 633 #ifdef SET_INEXACT
 634                                 else if (!b->x[0] && b->wds <= 1)
 635                                         inexact = 0;
 636 #endif
 637                                 *s++ = dig;
 638                                 goto ret;
 639                                 }
 640 #endif
 641                         if (j < 0 || j == 0 && mode != 1
 642 #ifndef ROUND_BIASED
 643                                                         && !(word1(d) & 1)
 644 #endif
 645                                         ) {
 646                                 if (!b->x[0] && b->wds <= 1) {
 647 #ifdef SET_INEXACT
 648                                         inexact = 0;
 649 #endif
 650                                         goto accept_dig;
 651                                         }
 652 #ifdef Honor_FLT_ROUNDS
 653                                 if (mode > 1)
 654                                  switch(Rounding) {
 655                                   case 0: goto accept_dig;
 656                                   case 2: goto keep_dig;
 657                                   }
 658 #endif /*Honor_FLT_ROUNDS*/
 659                                 if (j1 > 0) {
 660                                         b = lshift(b, 1);
 661                                         j1 = cmp(b, S);
 662                                         if ((j1 > 0 || j1 == 0 && dig & 1)
 663                                         && dig++ == '9')
 664                                                 goto round_9_up;
 665                                         }
 666  accept_dig:
 667                                 *s++ = dig;
 668                                 goto ret;
 669                                 }
 670                         if (j1 > 0) {
 671 #ifdef Honor_FLT_ROUNDS
 672                                 if (!Rounding)
 673                                         goto accept_dig;
 674 #endif
 675                                 if (dig == '9') { /* possible if i == 1 */
 676  round_9_up:
 677                                         *s++ = '9';
 678                                         goto roundoff;
 679                                         }
 680                                 *s++ = dig + 1;
 681                                 goto ret;
 682                                 }
 683 #ifdef Honor_FLT_ROUNDS
 684  keep_dig:
 685 #endif
 686                         *s++ = dig;
 687                         if (i == ilim)
 688                                 break;
 689                         b = multadd(b, 10, 0);
 690                         if (mlo == mhi)
 691                                 mlo = mhi = multadd(mhi, 10, 0);
 692                         else {
 693                                 mlo = multadd(mlo, 10, 0);
 694                                 mhi = multadd(mhi, 10, 0);
 695                                 }
 696                         }
 697                 }
 698         else
 699                 for(i = 1;; i++) {
 700                         *s++ = dig = quorem(b,S) + '0';
 701                         if (!b->x[0] && b->wds <= 1) {
 702 #ifdef SET_INEXACT
 703                                 inexact = 0;
 704 #endif
 705                                 goto ret;
 706                                 }
 707                         if (i >= ilim)
 708                                 break;
 709                         b = multadd(b, 10, 0);
 710                         }
 711
 712         /* Round off last digit */
 713
 714 #ifdef Honor_FLT_ROUNDS
 715         switch(Rounding) {
 716           case 0: goto trimzeros;
 717           case 2: goto roundoff;
 718           }
 719 #endif
 720         b = lshift(b, 1);
 721         j = cmp(b, S);
 722         if (j > 0 || j == 0 && dig & 1) {
 723  roundoff:
 724                 while(*--s == '9')
 725                         if (s == s0) {
 726                                 k++;
 727                                 *s++ = '1';
 728                                 goto ret;
 729                                 }
 730                 ++*s++;
 731                 }
 732         else {
 733  trimzeros:
 734                 while(*--s == '0');
 735                 s++;
 736                 }
 737  ret:
 738         Bfree(S);
 739         if (mhi) {
 740                 if (mlo && mlo != mhi)
 741                         Bfree(mlo);
 742                 Bfree(mhi);
 743                 }
 744  ret1:
 745 #ifdef SET_INEXACT
 746         if (inexact) {
 747                 if (!oldinexact) {
 748                         word0(d) = Exp_1 + (70 << Exp_shift);
 749                         word1(d) = 0;
 750                         dval(d) += 1.;
 751                         }
 752                 }
 753         else if (!oldinexact)
 754                 clear_inexact();
 755 #endif
 756         Bfree(b);
 757         *s = 0;
 758         *decpt = k + 1;
 759         if (rve)
 760                 *rve = s;
 761         return s0;
 762         }