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