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