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 #define Rounding rounding
  70 #undef Check_FLT_ROUNDS
  71 #define Check_FLT_ROUNDS
  72 #else
  73 #define Rounding Flt_Rounds
  74 #endif
  75
  76  char *
  77 dtoa
  78 #ifdef KR_headers
  79         (d, mode, ndigits, decpt, sign, rve)
  80         double d; int mode, ndigits, *decpt, *sign; char **rve;
  81 #else
  82         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  83 #endif
  84 {
  85  /*     Arguments ndigits, decpt, sign are similar to those
  86         of ecvt and fcvt; trailing zeros are suppressed from
  87         the returned string.  If not null, *rve is set to point
  88         to the end of the return value.  If d is +-Infinity or NaN,
  89         then *decpt is set to 9999.
  90
  91         mode:
  92                 0 ==> shortest string that yields d when read in
  93                         and rounded to nearest.
  94                 1 ==> like 0, but with Steele & White stopping rule;
  95                         e.g. with IEEE P754 arithmetic , mode 0 gives
  96                         1e23 whereas mode 1 gives 9.999999999999999e22.
  97                 2 ==> max(1,ndigits) significant digits.  This gives a
  98                         return value similar to that of ecvt, except
  99                         that trailing zeros are suppressed.
 100                 3 ==> through ndigits past the decimal point.  This
 101                         gives a return value similar to that from fcvt,
 102                         except that trailing zeros are suppressed, and
 103                         ndigits can be negative.
 104                 4,5 ==> similar to 2 and 3, respectively, but (in
 105                         round-nearest mode) with the tests of mode 0 to
 106                         possibly return a shorter string that rounds to d.
 107                         With IEEE arithmetic and compilation with
 108                         -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
 109                         as modes 2 and 3 when FLT_ROUNDS != 1.
 110                 6-9 ==> Debugging modes similar to mode - 4:  don't try
 111                         fast floating-point estimate (if applicable).
 112
 113                 Values of mode other than 0-9 are treated as mode 0.
 114
 115                 Sufficient space is allocated to the return value
 116                 to hold the suppressed trailing zeros.
 117         */
 118
 119         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
 120                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
 121                 spec_case, try_quick;
 122         Long L;
 123 #ifndef Sudden_Underflow
 124         int denorm;
 125         ULong x;
 126 #endif
 127         Bigint *b, *b1, *delta, *mlo, *mhi, *S;
 128         double d2, ds, eps;
 129         char *s, *s0;
 130 #ifdef Honor_FLT_ROUNDS
 131         int rounding;
 132 #endif
 133 #ifdef SET_INEXACT
 134         int inexact, oldinexact;
 135 #endif
 136
 137 #ifndef MULTIPLE_THREADS
 138         if (dtoa_result) {
 139                 freedtoa(dtoa_result);
 140                 dtoa_result = 0;
 141                 }
 142 #endif
 143
 144         if (word0(d) & Sign_bit) {
 145                 /* set sign for everything, including 0's and NaNs */
 146                 *sign = 1;
 147                 word0(d) &= ~Sign_bit;  /* clear sign bit */
 148                 }
 149         else
 150                 *sign = 0;
 151
 152 #if defined(IEEE_Arith) + defined(VAX)
 153 #ifdef IEEE_Arith
 154         if ((word0(d) & Exp_mask) == Exp_mask)
 155 #else
 156         if (word0(d)  == 0x8000)
 157 #endif
 158                 {
 159                 /* Infinity or NaN */
 160                 *decpt = 9999;
 161 #ifdef IEEE_Arith
 162                 if (!word1(d) && !(word0(d) & 0xfffff))
 163                         return nrv_alloc("Infinity", rve, 8);
 164 #endif
 165                 return nrv_alloc("NaN", rve, 3);
 166                 }
 167 #endif
 168 #ifdef IBM
 169         dval(d) += 0; /* normalize */
 170 #endif
 171         if (!dval(d)) {
 172                 *decpt = 1;
 173                 return nrv_alloc("0", rve, 1);
 174                 }
 175
 176 #ifdef SET_INEXACT
 177         try_quick = oldinexact = get_inexact();
 178         inexact = 1;
 179 #endif
 180 #ifdef Honor_FLT_ROUNDS
 181         if ((rounding = Flt_Rounds) >= 2) {
 182                 if (*sign)
 183                         rounding = rounding == 2 ? 0 : 2;
 184                 else
 185                         if (rounding != 2)
 186                                 rounding = 0;
 187                 }
 188 #endif
 189
 190         b = d2b(dval(d), &be, &bbits);
 191 #ifdef Sudden_Underflow
 192         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
 193 #else
 194         if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
 195 #endif
 196                 dval(d2) = dval(d);
 197                 word0(d2) &= Frac_mask1;
 198                 word0(d2) |= Exp_11;
 199 #ifdef IBM
 200                 if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
 201                         dval(d2) /= 1 << j;
 202 #endif
 203
 204                 /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
 205                  * log10(x)      =  log(x) / log(10)
 206                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 207                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
 208                  *
 209                  * This suggests computing an approximation k to log10(d) by
 210                  *
 211                  * k = (i - Bias)*0.301029995663981
 212                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 213                  *
 214                  * We want k to be too large rather than too small.
 215                  * The error in the first-order Taylor series approximation
 216                  * is in our favor, so we just round up the constant enough
 217                  * to compensate for any error in the multiplication of
 218                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 219                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 220                  * adding 1e-13 to the constant term more than suffices.
 221                  * Hence we adjust the constant term to 0.1760912590558.
 222                  * (We could get a more accurate k by invoking log10,
 223                  *  but this is probably not worthwhile.)
 224                  */
 225
 226                 i -= Bias;
 227 #ifdef IBM
 228                 i <<= 2;
 229                 i += j;
 230 #endif
 231 #ifndef Sudden_Underflow
 232                 denorm = 0;
 233                 }
 234         else {
 235                 /* d is denormalized */
 236
 237                 i = bbits + be + (Bias + (P-1) - 1);
 238                 x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
 239                             : word1(d) << 32 - i;
 240                 dval(d2) = x;
 241                 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
 242                 i -= (Bias + (P-1) - 1) + 1;
 243                 denorm = 1;
 244                 }
 245 #endif
 246         ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 247         k = (int)ds;
 248         if (ds < 0. && ds != k)
 249                 k--;    /* want k = floor(ds) */
 250         k_check = 1;
 251         if (k >= 0 && k <= Ten_pmax) {
 252                 if (dval(d) < tens[k])
 253                         k--;
 254                 k_check = 0;
 255                 }
 256         j = bbits - i - 1;
 257         if (j >= 0) {
 258                 b2 = 0;
 259                 s2 = j;
 260                 }
 261         else {
 262                 b2 = -j;
 263                 s2 = 0;
 264                 }
 265         if (k >= 0) {
 266                 b5 = 0;
 267                 s5 = k;
 268                 s2 += k;
 269                 }
 270         else {
 271                 b2 -= k;
 272                 b5 = -k;
 273                 s5 = 0;
 274                 }
 275         if (mode < 0 || mode > 9)
 276                 mode = 0;
 277
 278 #ifndef SET_INEXACT
 279 #ifdef Check_FLT_ROUNDS
 280         try_quick = Rounding == 1;
 281 #else
 282         try_quick = 1;
 283 #endif
 284 #endif /*SET_INEXACT*/
 285
 286         if (mode > 5) {
 287                 mode -= 4;
 288                 try_quick = 0;
 289                 }
 290         leftright = 1;
 291         switch(mode) {
 292                 case 0:
 293                 case 1:
 294                         ilim = ilim1 = -1;
 295                         i = 18;
 296                         ndigits = 0;
 297                         break;
 298                 case 2:
 299                         leftright = 0;
 300                         /* no break */
 301                 case 4:
 302                         if (ndigits <= 0)
 303                                 ndigits = 1;
 304                         ilim = ilim1 = i = ndigits;
 305                         break;
 306                 case 3:
 307                         leftright = 0;
 308                         /* no break */
 309                 case 5:
 310                         i = ndigits + k + 1;
 311                         ilim = i;
 312                         ilim1 = i - 1;
 313                         if (i <= 0)
 314                                 i = 1;
 315                 }
 316         s = s0 = rv_alloc(i);
 317
 318 #ifdef Honor_FLT_ROUNDS
 319         if (mode > 1 && rounding != 1)
 320                 leftright = 0;
 321 #endif
 322
 323         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
 324
 325                 /* Try to get by with floating-point arithmetic. */
 326
 327                 i = 0;
 328                 dval(d2) = dval(d);
 329                 k0 = k;
 330                 ilim0 = ilim;
 331                 ieps = 2; /* conservative */
 332                 if (k > 0) {
 333                         ds = tens[k&0xf];
 334                         j = k >> 4;
 335                         if (j & Bletch) {
 336                                 /* prevent overflows */
 337                                 j &= Bletch - 1;
 338                                 dval(d) /= bigtens[n_bigtens-1];
 339                                 ieps++;
 340                                 }
 341                         for(; j; j >>= 1, i++)
 342                                 if (j & 1) {
 343                                         ieps++;
 344                                         ds *= bigtens[i];
 345                                         }
 346                         dval(d) /= ds;
 347                         }
 348                 else if (( j1 = -k )!=0) {
 349                         dval(d) *= tens[j1 & 0xf];
 350                         for(j = j1 >> 4; j; j >>= 1, i++)
 351                                 if (j & 1) {
 352                                         ieps++;
 353                                         dval(d) *= bigtens[i];
 354                                         }
 355                         }
 356                 if (k_check && dval(d) < 1. && ilim > 0) {
 357                         if (ilim1 <= 0)
 358                                 goto fast_failed;
 359                         ilim = ilim1;
 360                         k--;
 361                         dval(d) *= 10.;
 362                         ieps++;
 363                         }
 364                 dval(eps) = ieps*dval(d) + 7.;
 365                 word0(eps) -= (P-1)*Exp_msk1;
 366                 if (ilim == 0) {
 367                         S = mhi = 0;
 368                         dval(d) -= 5.;
 369                         if (dval(d) > dval(eps))
 370                                 goto one_digit;
 371                         if (dval(d) < -dval(eps))
 372                                 goto no_digits;
 373                         goto fast_failed;
 374                         }
 375 #ifndef No_leftright
 376                 if (leftright) {
 377                         /* Use Steele & White method of only
 378                          * generating digits needed.
 379                          */
 380                         dval(eps) = 0.5/tens[ilim-1] - dval(eps);
 381                         for(i = 0;;) {
 382                                 L = dval(d);
 383                                 dval(d) -= L;
 384                                 *s++ = '0' + (int)L;
 385                                 if (dval(d) < dval(eps))
 386                                         goto ret1;
 387                                 if (1. - dval(d) < dval(eps))
 388                                         goto bump_up;
 389                                 if (++i >= ilim)
 390                                         break;
 391                                 dval(eps) *= 10.;
 392                                 dval(d) *= 10.;
 393                                 }
 394                         }
 395                 else {
 396 #endif
 397                         /* Generate ilim digits, then fix them up. */
 398                         dval(eps) *= tens[ilim-1];
 399                         for(i = 1;; i++, dval(d) *= 10.) {
 400                                 L = (Long)(dval(d));
 401                                 if (!(dval(d) -= L))
 402                                         ilim = i;
 403                                 *s++ = '0' + (int)L;
 404                                 if (i == ilim) {
 405                                         if (dval(d) > 0.5 + dval(eps))
 406                                                 goto bump_up;
 407                                         else if (dval(d) < 0.5 - dval(eps)) {
 408                                                 while(*--s == '0');
 409                                                 s++;
 410                                                 goto ret1;
 411                                                 }
 412                                         break;
 413                                         }
 414                                 }
 415 #ifndef No_leftright
 416                         }
 417 #endif
 418  fast_failed:
 419                 s = s0;
 420                 dval(d) = dval(d2);
 421                 k = k0;
 422                 ilim = ilim0;
 423                 }
 424
 425         /* Do we have a "small" integer? */
 426
 427         if (be >= 0 && k <= Int_max) {
 428                 /* Yes. */
 429                 ds = tens[k];
 430                 if (ndigits < 0 && ilim <= 0) {
 431                         S = mhi = 0;
 432                         if (ilim < 0 || dval(d) <= 5*ds)
 433                                 goto no_digits;
 434                         goto one_digit;
 435                         }
 436                 for(i = 1;; i++, dval(d) *= 10.) {
 437                         L = (Long)(dval(d) / ds);
 438                         dval(d) -= L*ds;
 439 #ifdef Check_FLT_ROUNDS
 440                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 441                         if (dval(d) < 0) {
 442                                 L--;
 443                                 dval(d) += ds;
 444                                 }
 445 #endif
 446                         *s++ = '0' + (int)L;
 447                         if (!dval(d)) {
 448 #ifdef SET_INEXACT
 449                                 inexact = 0;
 450 #endif
 451                                 break;
 452                                 }
 453                         if (i == ilim) {
 454 #ifdef Honor_FLT_ROUNDS
 455                                 if (mode > 1)
 456                                 switch(rounding) {
 457                                   case 0: goto ret1;
 458                                   case 2: goto bump_up;
 459                                   }
 460 #endif
 461                                 dval(d) += dval(d);
 462                                 if (dval(d) > ds || dval(d) == ds && L & 1) {
 463  bump_up:
 464                                         while(*--s == '9')
 465                                                 if (s == s0) {
 466                                                         k++;
 467                                                         *s = '0';
 468                                                         break;
 469                                                         }
 470                                         ++*s++;
 471                                         }
 472                                 break;
 473                                 }
 474                         }
 475                 goto ret1;
 476                 }
 477
 478         m2 = b2;
 479         m5 = b5;
 480         mhi = mlo = 0;
 481         if (leftright) {
 482                 i =
 483 #ifndef Sudden_Underflow
 484                         denorm ? be + (Bias + (P-1) - 1 + 1) :
 485 #endif
 486 #ifdef IBM
 487                         1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
 488 #else
 489                         1 + P - bbits;
 490 #endif
 491                 b2 += i;
 492                 s2 += i;
 493                 mhi = i2b(1);
 494                 }
 495         if (m2 > 0 && s2 > 0) {
 496                 i = m2 < s2 ? m2 : s2;
 497                 b2 -= i;
 498                 m2 -= i;
 499                 s2 -= i;
 500                 }
 501         if (b5 > 0) {
 502                 if (leftright) {
 503                         if (m5 > 0) {
 504                                 mhi = pow5mult(mhi, m5);
 505                                 b1 = mult(mhi, b);
 506                                 Bfree(b);
 507                                 b = b1;
 508                                 }
 509                         if (( j = b5 - m5 )!=0)
 510                                 b = pow5mult(b, j);
 511                         }
 512                 else
 513                         b = pow5mult(b, b5);
 514                 }
 515         S = i2b(1);
 516         if (s5 > 0)
 517                 S = pow5mult(S, s5);
 518
 519         /* Check for special case that d is a normalized power of 2. */
 520
 521         spec_case = 0;
 522         if ((mode < 2 || leftright)
 523 #ifdef Honor_FLT_ROUNDS
 524                         && rounding == 1
 525 #endif
 526                                 ) {
 527                 if (!word1(d) && !(word0(d) & Bndry_mask)
 528 #ifndef Sudden_Underflow
 529                  && word0(d) & (Exp_mask & ~Exp_msk1)
 530 #endif
 531                                 ) {
 532                         /* The special case */
 533                         b2 += Log2P;
 534                         s2 += Log2P;
 535                         spec_case = 1;
 536                         }
 537                 }
 538
 539         /* Arrange for convenient computation of quotients:
 540          * shift left if necessary so divisor has 4 leading 0 bits.
 541          *
 542          * Perhaps we should just compute leading 28 bits of S once
 543          * and for all and pass them and a shift to quorem, so it
 544          * can do shifts and ors to compute the numerator for q.
 545          */
 546 #ifdef Pack_32
 547         if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
 548                 i = 32 - i;
 549 #else
 550         if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
 551                 i = 16 - i;
 552 #endif
 553         if (i > 4) {
 554                 i -= 4;
 555                 b2 += i;
 556                 m2 += i;
 557                 s2 += i;
 558                 }
 559         else if (i < 4) {
 560                 i += 28;
 561                 b2 += i;
 562                 m2 += i;
 563                 s2 += i;
 564                 }
 565         if (b2 > 0)
 566                 b = lshift(b, b2);
 567         if (s2 > 0)
 568                 S = lshift(S, s2);
 569         if (k_check) {
 570                 if (cmp(b,S) < 0) {
 571                         k--;
 572                         b = multadd(b, 10, 0);  /* we botched the k estimate */
 573                         if (leftright)
 574                                 mhi = multadd(mhi, 10, 0);
 575                         ilim = ilim1;
 576                         }
 577                 }
 578         if (ilim <= 0 && (mode == 3 || mode == 5)) {
 579                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 580                         /* no digits, fcvt style */
 581  no_digits:
 582                         k = -1 - ndigits;
 583                         goto ret;
 584                         }
 585  one_digit:
 586                 *s++ = '1';
 587                 k++;
 588                 goto ret;
 589                 }
 590         if (leftright) {
 591                 if (m2 > 0)
 592                         mhi = lshift(mhi, m2);
 593
 594                 /* Compute mlo -- check for special case
 595                  * that d is a normalized power of 2.
 596                  */
 597
 598                 mlo = mhi;
 599                 if (spec_case) {
 600                         mhi = Balloc(mhi->k);
 601                         Bcopy(mhi, mlo);
 602                         mhi = lshift(mhi, Log2P);
 603                         }
 604
 605                 for(i = 1;;i++) {
 606                         dig = quorem(b,S) + '0';
 607                         /* Do we yet have the shortest decimal string
 608                          * that will round to d?
 609                          */
 610                         j = cmp(b, mlo);
 611                         delta = diff(S, mhi);
 612                         j1 = delta->sign ? 1 : cmp(b, delta);
 613                         Bfree(delta);
 614 #ifndef ROUND_BIASED
 615                         if (j1 == 0 && mode != 1 && !(word1(d) & 1)
 616 #ifdef Honor_FLT_ROUNDS
 617                                 && rounding >= 1
 618 #endif
 619                                                                    ) {
 620                                 if (dig == '9')
 621                                         goto round_9_up;
 622                                 if (j > 0)
 623                                         dig++;
 624 #ifdef SET_INEXACT
 625                                 else if (!b->x[0] && b->wds <= 1)
 626                                         inexact = 0;
 627 #endif
 628                                 *s++ = dig;
 629                                 goto ret;
 630                                 }
 631 #endif
 632                         if (j < 0 || j == 0 && mode != 1
 633 #ifndef ROUND_BIASED
 634                                                         && !(word1(d) & 1)
 635 #endif
 636                                         ) {
 637                                 if (!b->x[0] && b->wds <= 1) {
 638 #ifdef SET_INEXACT
 639                                         inexact = 0;
 640 #endif
 641                                         goto accept_dig;
 642                                         }
 643 #ifdef Honor_FLT_ROUNDS
 644                                 if (mode > 1)
 645                                  switch(rounding) {
 646                                   case 0: goto accept_dig;
 647                                   case 2: goto keep_dig;
 648                                   }
 649 #endif /*Honor_FLT_ROUNDS*/
 650                                 if (j1 > 0) {
 651                                         b = lshift(b, 1);
 652                                         j1 = cmp(b, S);
 653                                         if ((j1 > 0 || j1 == 0 && dig & 1)
 654                                         && dig++ == '9')
 655                                                 goto round_9_up;
 656                                         }
 657  accept_dig:
 658                                 *s++ = dig;
 659                                 goto ret;
 660                                 }
 661                         if (j1 > 0) {
 662 #ifdef Honor_FLT_ROUNDS
 663                                 if (!rounding)
 664                                         goto accept_dig;
 665 #endif
 666                                 if (dig == '9') { /* possible if i == 1 */
 667  round_9_up:
 668                                         *s++ = '9';
 669                                         goto roundoff;
 670                                         }
 671                                 *s++ = dig + 1;
 672                                 goto ret;
 673                                 }
 674 #ifdef Honor_FLT_ROUNDS
 675  keep_dig:
 676 #endif
 677                         *s++ = dig;
 678                         if (i == ilim)
 679                                 break;
 680                         b = multadd(b, 10, 0);
 681                         if (mlo == mhi)
 682                                 mlo = mhi = multadd(mhi, 10, 0);
 683                         else {
 684                                 mlo = multadd(mlo, 10, 0);
 685                                 mhi = multadd(mhi, 10, 0);
 686                                 }
 687                         }
 688                 }
 689         else
 690                 for(i = 1;; i++) {
 691                         *s++ = dig = quorem(b,S) + '0';
 692                         if (!b->x[0] && b->wds <= 1) {
 693 #ifdef SET_INEXACT
 694                                 inexact = 0;
 695 #endif
 696                                 goto ret;
 697                                 }
 698                         if (i >= ilim)
 699                                 break;
 700                         b = multadd(b, 10, 0);
 701                         }
 702
 703         /* Round off last digit */
 704
 705 #ifdef Honor_FLT_ROUNDS
 706         switch(rounding) {
 707           case 0: goto trimzeros;
 708           case 2: goto roundoff;
 709           }
 710 #endif
 711         b = lshift(b, 1);
 712         j = cmp(b, S);
 713         if (j > 0 || j == 0 && dig & 1) {
 714  roundoff:
 715                 while(*--s == '9')
 716                         if (s == s0) {
 717                                 k++;
 718                                 *s++ = '1';
 719                                 goto ret;
 720                                 }
 721                 ++*s++;
 722                 }
 723         else {
 724  trimzeros:
 725                 while(*--s == '0');
 726                 s++;
 727                 }
 728  ret:
 729         Bfree(S);
 730         if (mhi) {
 731                 if (mlo && mlo != mhi)
 732                         Bfree(mlo);
 733                 Bfree(mhi);
 734                 }
 735  ret1:
 736 #ifdef SET_INEXACT
 737         if (inexact) {
 738                 if (!oldinexact) {
 739                         word0(d) = Exp_1 + (70 << Exp_shift);
 740                         word1(d) = 0;
 741                         dval(d) += 1.;
 742                         }
 743                 }
 744         else if (!oldinexact)
 745                 clear_inexact();
 746 #endif
 747         Bfree(b);
 748         *s = 0;
 749         *decpt = k + 1;
 750         if (rve)
 751                 *rve = s;
 752         return s0;
 753         }