StdLib/LibC/gdtoa/dtoa.c

   1 /* $NetBSD: dtoa.c,v 1.3.4.1.4.1 2008/04/08 21:10:55 jdc Exp $ */
   2
   3 /****************************************************************
   4
   5 The author of this software is David M. Gay.
   6
   7 Copyright (C) 1998, 1999 by Lucent Technologies
   8 All Rights Reserved
   9
  10 Permission to use, copy, modify, and distribute this software and
  11 its documentation for any purpose and without fee is hereby
  12 granted, provided that the above copyright notice appear in all
  13 copies and that both that the copyright notice and this
  14 permission notice and warranty disclaimer appear in supporting
  15 documentation, and that the name of Lucent or any of its entities
  16 not be used in advertising or publicity pertaining to
  17 distribution of the software without specific, written prior
  18 permission.
  19
  20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  27 THIS SOFTWARE.
  28
  29 ****************************************************************/
  30
  31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  32  * with " at " changed at "@" and " dot " changed to ".").  */
  33 #include  <LibConfig.h>
  34
  35 #include "gdtoaimp.h"
  36
  37 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  38  *
  39  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  40  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
  41  *
  42  * Modifications:
  43  *  1. Rather than iterating, we use a simple numeric overestimate
  44  *     to determine k = floor(log10(d)).  We scale relevant
  45  *     quantities using O(log2(k)) rather than O(k) multiplications.
  46  *  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  47  *     try to generate digits strictly left to right.  Instead, we
  48  *     compute with fewer bits and propagate the carry if necessary
  49  *     when rounding the final digit up.  This is often faster.
  50  *  3. Under the assumption that input will be rounded nearest,
  51  *     mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  52  *     That is, we allow equality in stopping tests when the
  53  *     round-nearest rule will give the same floating-point value
  54  *     as would satisfaction of the stopping test with strict
  55  *     inequality.
  56  *  4. We remove common factors of powers of 2 from relevant
  57  *     quantities.
  58  *  5. When converting floating-point integers less than 1e16,
  59  *     we use floating-point arithmetic rather than resorting
  60  *     to multiple-precision integers.
  61  *  6. When asked to produce fewer than 15 digits, we first try
  62  *     to get by with floating-point arithmetic; we resort to
  63  *     multiple-precision integer arithmetic only if we cannot
  64  *     guarantee that the floating-point calculation has given
  65  *     the correctly rounded result.  For k requested digits and
  66  *     "uniformly" distributed input, the probability is
  67  *     something like 10^(k-15) that we must resort to the Long
  68  *     calculation.
  69  */
  70
  71 #ifdef Honor_FLT_ROUNDS
  72 #define Rounding rounding
  73 #undef Check_FLT_ROUNDS
  74 #define Check_FLT_ROUNDS
  75 #else
  76 #define Rounding Flt_Rounds
  77 #endif
  78
  79 #if defined(_MSC_VER)           /* Handle Microsoft VC++ compiler specifics. */
  80 // Disable: warning C4700: uninitialized local variable 'xx' used
  81 #pragma warning ( disable : 4700 )
  82 #endif  /* defined(_MSC_VER) */
  83
  84  char *
  85 dtoa
  86 #ifdef KR_headers
  87   (d, mode, ndigits, decpt, sign, rve)
  88   double d; int mode, ndigits, *decpt, *sign; char **rve;
  89 #else
  90   (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  91 #endif
  92 {
  93  /* Arguments ndigits, decpt, sign are similar to those
  94   of ecvt and fcvt; trailing zeros are suppressed from
  95   the returned string.  If not null, *rve is set to point
  96   to the end of the return value.  If d is +-Infinity or NaN,
  97   then *decpt is set to 9999.
  98
  99   mode:
 100     0 ==> shortest string that yields d when read in
 101       and rounded to nearest.
 102     1 ==> like 0, but with Steele & White stopping rule;
 103       e.g. with IEEE P754 arithmetic , mode 0 gives
 104       1e23 whereas mode 1 gives 9.999999999999999e22.
 105     2 ==> max(1,ndigits) significant digits.  This gives a
 106       return value similar to that of ecvt, except
 107       that trailing zeros are suppressed.
 108     3 ==> through ndigits past the decimal point.  This
 109       gives a return value similar to that from fcvt,
 110       except that trailing zeros are suppressed, and
 111       ndigits can be negative.
 112     4,5 ==> similar to 2 and 3, respectively, but (in
 113       round-nearest mode) with the tests of mode 0 to
 114       possibly return a shorter string that rounds to d.
 115       With IEEE arithmetic and compilation with
 116       -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
 117       as modes 2 and 3 when FLT_ROUNDS != 1.
 118     6-9 ==> Debugging modes similar to mode - 4:  don't try
 119       fast floating-point estimate (if applicable).
 120
 121     Values of mode other than 0-9 are treated as mode 0.
 122
 123     Sufficient space is allocated to the return value
 124     to hold the suppressed trailing zeros.
 125   */
 126
 127   int bbits, b2, b5, be, dig, i, ieps, ilim0,
 128     j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5,
 129     spec_case, try_quick;
 130   int ilim = 0, ilim1 = 0; /* pacify gcc */
 131   Long L;
 132 #ifndef Sudden_Underflow
 133   int denorm;
 134   ULong x;
 135 #endif
 136   Bigint *b, *b1, *delta, *mhi, *S;
 137   Bigint *mlo = NULL; /* pacify gcc */
 138   double d2, ds, eps;
 139   char *s, *s0;
 140 #ifdef Honor_FLT_ROUNDS
 141   int rounding;
 142 #endif
 143 #ifdef SET_INEXACT
 144   int inexact, oldinexact;
 145 #endif
 146
 147 #ifndef MULTIPLE_THREADS
 148   if (dtoa_result) {
 149     freedtoa(dtoa_result);
 150     dtoa_result = 0;
 151     }
 152 #endif
 153
 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 = Flt_Rounds) >= 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   if (b == NULL)
 202     return NULL;
 203 #ifdef Sudden_Underflow
 204   i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
 205 #else
 206   if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
 207 #endif
 208     dval(d2) = dval(d);
 209     word0(d2) &= Frac_mask1;
 210     word0(d2) |= Exp_11;
 211 #ifdef IBM
 212     if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
 213       dval(d2) /= 1 << j;
 214 #endif
 215
 216     /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
 217      * log10(x)  =  log(x) / log(10)
 218      *    ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 219      * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
 220      *
 221      * This suggests computing an approximation k to log10(d) by
 222      *
 223      * k = (i - Bias)*0.301029995663981
 224      *  + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 225      *
 226      * We want k to be too large rather than too small.
 227      * The error in the first-order Taylor series approximation
 228      * is in our favor, so we just round up the constant enough
 229      * to compensate for any error in the multiplication of
 230      * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 231      * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 232      * adding 1e-13 to the constant term more than suffices.
 233      * Hence we adjust the constant term to 0.1760912590558.
 234      * (We could get a more accurate k by invoking log10,
 235      *  but this is probably not worthwhile.)
 236      */
 237
 238     i -= Bias;
 239 #ifdef IBM
 240     i <<= 2;
 241     i += j;
 242 #endif
 243 #ifndef Sudden_Underflow
 244     denorm = 0;
 245     }
 246   else {
 247     /* d is denormalized */
 248
 249     i = bbits + be + (Bias + (P-1) - 1);
 250     x = i > 32  ? word0(d) << (64 - i) | word1(d) >> (i - 32)
 251           : word1(d) << (32 - i);
 252     dval(d2) = (double)x;
 253     word0(d2) -= 31*Exp_msk1; /* adjust exponent */
 254     i -= (Bias + (P-1) - 1) + 1;
 255     denorm = 1;
 256     }
 257 #endif
 258   ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 259   k = (int)ds;
 260   if (ds < 0. && ds != k)
 261     k--;  /* want k = floor(ds) */
 262   k_check = 1;
 263   if (k >= 0 && k <= Ten_pmax) {
 264     if (dval(d) < tens[k])
 265       k--;
 266     k_check = 0;
 267     }
 268   j = bbits - i - 1;
 269   if (j >= 0) {
 270     b2 = 0;
 271     s2 = j;
 272     }
 273   else {
 274     b2 = -j;
 275     s2 = 0;
 276     }
 277   if (k >= 0) {
 278     b5 = 0;
 279     s5 = k;
 280     s2 += k;
 281     }
 282   else {
 283     b2 -= k;
 284     b5 = -k;
 285     s5 = 0;
 286     }
 287   if (mode < 0 || mode > 9)
 288     mode = 0;
 289
 290 #ifndef SET_INEXACT
 291 #ifdef Check_FLT_ROUNDS
 292   try_quick = Rounding == 1;
 293 #else
 294   try_quick = 1;
 295 #endif
 296 #endif /*SET_INEXACT*/
 297
 298   if (mode > 5) {
 299     mode -= 4;
 300     try_quick = 0;
 301     }
 302   leftright = 1;
 303   switch(mode) {
 304     case 0:
 305     case 1:
 306       ilim = ilim1 = -1;
 307       i = 18;
 308       ndigits = 0;
 309       break;
 310     case 2:
 311       leftright = 0;
 312       /* FALLTHROUGH */
 313     case 4:
 314       if (ndigits <= 0)
 315         ndigits = 1;
 316       ilim = ilim1 = i = ndigits;
 317       break;
 318     case 3:
 319       leftright = 0;
 320       /* FALLTHROUGH */
 321     case 5:
 322       i = ndigits + k + 1;
 323       ilim = i;
 324       ilim1 = i - 1;
 325       if (i <= 0)
 326         i = 1;
 327     }
 328   s = s0 = rv_alloc((size_t)i);
 329   if (s == NULL)
 330     return NULL;
 331
 332 #ifdef Honor_FLT_ROUNDS
 333   if (mode > 1 && rounding != 1)
 334     leftright = 0;
 335 #endif
 336
 337   if (ilim >= 0 && ilim <= Quick_max && try_quick) {
 338
 339     /* Try to get by with floating-point arithmetic. */
 340
 341     i = 0;
 342     dval(d2) = dval(d);
 343     k0 = k;
 344     ilim0 = ilim;
 345     ieps = 2; /* conservative */
 346     if (k > 0) {
 347       ds = tens[k&0xf];
 348       j = (unsigned int)k >> 4;
 349       if (j & Bletch) {
 350         /* prevent overflows */
 351         j &= Bletch - 1;
 352         dval(d) /= bigtens[n_bigtens-1];
 353         ieps++;
 354         }
 355       for(; j; j = (unsigned int)j >> 1, i++)
 356         if (j & 1) {
 357           ieps++;
 358           ds *= bigtens[i];
 359           }
 360       dval(d) /= ds;
 361       }
 362     else if (( jj1 = -k )!=0) {
 363       dval(d) *= tens[jj1 & 0xf];
 364       for(j = jj1 >> 4; j; j >>= 1, i++)
 365         if (j & 1) {
 366           ieps++;
 367           dval(d) *= bigtens[i];
 368           }
 369       }
 370     if (k_check && dval(d) < 1. && ilim > 0) {
 371       if (ilim1 <= 0)
 372         goto fast_failed;
 373       ilim = ilim1;
 374       k--;
 375       dval(d) *= 10.;
 376       ieps++;
 377       }
 378     dval(eps) = ieps*dval(d) + 7.;
 379     word0(eps) -= (P-1)*Exp_msk1;
 380     if (ilim == 0) {
 381       S = mhi = 0;
 382       dval(d) -= 5.;
 383       if (dval(d) > dval(eps))
 384         goto one_digit;
 385       if (dval(d) < -dval(eps))
 386         goto no_digits;
 387       goto fast_failed;
 388       }
 389 #ifndef No_leftright
 390     if (leftright) {
 391       /* Use Steele & White method of only
 392        * generating digits needed.
 393        */
 394       dval(eps) = 0.5/tens[ilim-1] - dval(eps);
 395       for(i = 0;;) {
 396         L = (INT32)dval(d);
 397         dval(d) -= L;
 398         *s++ = (char)('0' + (int)L);
 399         if (dval(d) < dval(eps))
 400           goto ret1;
 401         if (1. - dval(d) < dval(eps))
 402           goto bump_up;
 403         if (++i >= ilim)
 404           break;
 405         dval(eps) *= 10.;
 406         dval(d) *= 10.;
 407         }
 408       }
 409     else {
 410 #endif
 411       /* Generate ilim digits, then fix them up. */
 412       dval(eps) *= tens[ilim-1];
 413       for(i = 1;; i++, dval(d) *= 10.) {
 414         L = (Long)(dval(d));
 415         if (!(dval(d) -= L))
 416           ilim = i;
 417         *s++ = (char)('0' + (int)L);
 418         if (i == ilim) {
 419           if (dval(d) > 0.5 + dval(eps))
 420             goto bump_up;
 421           else if (dval(d) < 0.5 - dval(eps)) {
 422             while(*--s == '0');
 423             s++;
 424             goto ret1;
 425             }
 426           break;
 427           }
 428         }
 429 #ifndef No_leftright
 430       }
 431 #endif
 432  fast_failed:
 433     s = s0;
 434     dval(d) = dval(d2);
 435     k = k0;
 436     ilim = ilim0;
 437     }
 438
 439   /* Do we have a "small" integer? */
 440
 441   if (be >= 0 && k <= Int_max) {
 442     /* Yes. */
 443     ds = tens[k];
 444     if (ndigits < 0 && ilim <= 0) {
 445       S = mhi = 0;
 446       if (ilim < 0 || dval(d) <= 5*ds)
 447         goto no_digits;
 448       goto one_digit;
 449       }
 450     for(i = 1;; i++, dval(d) *= 10.) {
 451       L = (Long)(dval(d) / ds);
 452       dval(d) -= L*ds;
 453 #ifdef Check_FLT_ROUNDS
 454       /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 455       if (dval(d) < 0) {
 456         L--;
 457         dval(d) += ds;
 458         }
 459 #endif
 460       *s++ = (char)('0' + (int)L);
 461       if (!dval(d)) {
 462 #ifdef SET_INEXACT
 463         inexact = 0;
 464 #endif
 465         break;
 466         }
 467       if (i == ilim) {
 468 #ifdef Honor_FLT_ROUNDS
 469         if (mode > 1)
 470         switch(rounding) {
 471           case 0: goto ret1;
 472           case 2: goto bump_up;
 473           }
 474 #endif
 475         dval(d) += dval(d);
 476         if (dval(d) > ds || (dval(d) == ds && L & 1)) {
 477  bump_up:
 478           while(*--s == '9')
 479             if (s == s0) {
 480               k++;
 481               *s = '0';
 482               break;
 483               }
 484           ++*s++;
 485           }
 486         break;
 487         }
 488       }
 489     goto ret1;
 490     }
 491
 492   m2 = b2;
 493   m5 = b5;
 494   mhi = mlo = 0;
 495   if (leftright) {
 496     i =
 497 #ifndef Sudden_Underflow
 498       denorm ? be + (Bias + (P-1) - 1 + 1) :
 499 #endif
 500 #ifdef IBM
 501       1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
 502 #else
 503       1 + P - bbits;
 504 #endif
 505     b2 += i;
 506     s2 += i;
 507     mhi = i2b(1);
 508     if (mhi == NULL)
 509       return NULL;
 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         if (mhi == NULL)
 522           return NULL;
 523         b1 = mult(mhi, b);
 524         if (b1 == NULL)
 525           return NULL;
 526         Bfree(b);
 527         b = b1;
 528         }
 529       if (( j = b5 - m5 )!=0)
 530         b = pow5mult(b, j);
 531         if (b == NULL)
 532           return NULL;
 533       }
 534     else
 535       b = pow5mult(b, b5);
 536       if (b == NULL)
 537         return NULL;
 538     }
 539   S = i2b(1);
 540   if (S == NULL)
 541     return NULL;
 542   if (s5 > 0) {
 543     S = pow5mult(S, s5);
 544     if (S == NULL)
 545       return NULL;
 546   }
 547
 548   /* Check for special case that d is a normalized power of 2. */
 549
 550   spec_case = 0;
 551   if ((mode < 2 || leftright)
 552 #ifdef Honor_FLT_ROUNDS
 553       && rounding == 1
 554 #endif
 555         ) {
 556     if (!word1(d) && !(word0(d) & Bndry_mask)
 557 #ifndef Sudden_Underflow
 558      && word0(d) & (Exp_mask & ~Exp_msk1)
 559 #endif
 560         ) {
 561       /* The special case */
 562       b2 += Log2P;
 563       s2 += Log2P;
 564       spec_case = 1;
 565       }
 566     }
 567
 568   /* Arrange for convenient computation of quotients:
 569    * shift left if necessary so divisor has 4 leading 0 bits.
 570    *
 571    * Perhaps we should just compute leading 28 bits of S once
 572    * and for all and pass them and a shift to quorem, so it
 573    * can do shifts and ors to compute the numerator for q.
 574    */
 575 #ifdef Pack_32
 576   if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
 577     i = 32 - i;
 578 #else
 579   if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
 580     i = 16 - i;
 581 #endif
 582   if (i > 4) {
 583     i -= 4;
 584     b2 += i;
 585     m2 += i;
 586     s2 += i;
 587     }
 588   else if (i < 4) {
 589     i += 28;
 590     b2 += i;
 591     m2 += i;
 592     s2 += i;
 593     }
 594   if (b2 > 0) {
 595     b = lshift(b, b2);
 596     if (b == NULL)
 597       return NULL;
 598   }
 599   if (s2 > 0) {
 600     S = lshift(S, s2);
 601     if (S == NULL)
 602       return NULL;
 603   }
 604   if (k_check) {
 605     if (cmp(b,S) < 0) {
 606       k--;
 607       b = multadd(b, 10, 0);  /* we botched the k estimate */
 608       if (b == NULL)
 609         return NULL;
 610       if (leftright) {
 611         mhi = multadd(mhi, 10, 0);
 612         if (mhi == NULL)
 613           return NULL;
 614       }
 615       ilim = ilim1;
 616       }
 617     }
 618   if (ilim <= 0 && (mode == 3 || mode == 5)) {
 619     if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
 620       /* no digits, fcvt style */
 621  no_digits:
 622       k = -1 - ndigits;
 623       goto ret;
 624       }
 625  one_digit:
 626     *s++ = '1';
 627     k++;
 628     goto ret;
 629     }
 630   if (leftright) {
 631     if (m2 > 0) {
 632       mhi = lshift(mhi, m2);
 633       if (mhi == NULL)
 634         return NULL;
 635     }
 636
 637     /* Compute mlo -- check for special case
 638      * that d is a normalized power of 2.
 639      */
 640
 641     mlo = mhi;
 642     if (spec_case) {
 643       mhi = Balloc(mhi->k);
 644       if (mhi == NULL)
 645         return NULL;
 646       Bcopy(mhi, mlo);
 647       mhi = lshift(mhi, Log2P);
 648       if (mhi == NULL)
 649         return NULL;
 650       }
 651
 652     for(i = 1;;i++) {
 653       dig = quorem(b,S) + '0';
 654       /* Do we yet have the shortest decimal string
 655        * that will round to d?
 656        */
 657       j = cmp(b, mlo);
 658       delta = diff(S, mhi);
 659       if (delta == NULL)
 660         return NULL;
 661       jj1 = delta->sign ? 1 : cmp(b, delta);
 662       Bfree(delta);
 663 #ifndef ROUND_BIASED
 664       if (jj1 == 0 && mode != 1 && !(word1(d) & 1)
 665 #ifdef Honor_FLT_ROUNDS
 666         && rounding >= 1
 667 #endif
 668                    ) {
 669         if (dig == '9')
 670           goto round_9_up;
 671         if (j > 0)
 672           dig++;
 673 #ifdef SET_INEXACT
 674         else if (!b->x[0] && b->wds <= 1)
 675           inexact = 0;
 676 #endif
 677         *s++ = (char)dig;
 678         goto ret;
 679         }
 680 #endif
 681       if (j < 0 || (j == 0 && mode != 1
 682 #ifndef ROUND_BIASED
 683               && !(word1(d) & 1)
 684 #endif
 685           )) {
 686         if (!b->x[0] && b->wds <= 1) {
 687 #ifdef SET_INEXACT
 688           inexact = 0;
 689 #endif
 690           goto accept_dig;
 691           }
 692 #ifdef Honor_FLT_ROUNDS
 693         if (mode > 1)
 694          switch(rounding) {
 695           case 0: goto accept_dig;
 696           case 2: goto keep_dig;
 697           }
 698 #endif /*Honor_FLT_ROUNDS*/
 699         if (jj1 > 0) {
 700           b = lshift(b, 1);
 701           if (b == NULL)
 702             return NULL;
 703           jj1 = cmp(b, S);
 704           if ((jj1 > 0 || (jj1 == 0 && dig & 1))
 705           && dig++ == '9')
 706             goto round_9_up;
 707           }
 708  accept_dig:
 709         *s++ = (char)dig;
 710         goto ret;
 711         }
 712       if (jj1 > 0) {
 713 #ifdef Honor_FLT_ROUNDS
 714         if (!rounding)
 715           goto accept_dig;
 716 #endif
 717         if (dig == '9') { /* possible if i == 1 */
 718  round_9_up:
 719           *s++ = '9';
 720           goto roundoff;
 721           }
 722         *s++ = (char)(dig + 1);
 723         goto ret;
 724         }
 725 #ifdef Honor_FLT_ROUNDS
 726  keep_dig:
 727 #endif
 728       *s++ = (char)dig;
 729       if (i == ilim)
 730         break;
 731       b = multadd(b, 10, 0);
 732       if (b == NULL)
 733         return NULL;
 734       if (mlo == mhi) {
 735         mlo = mhi = multadd(mhi, 10, 0);
 736         if (mlo == NULL)
 737           return NULL;
 738         }
 739       else {
 740         mlo = multadd(mlo, 10, 0);
 741         if (mlo == NULL)
 742           return NULL;
 743         mhi = multadd(mhi, 10, 0);
 744         if (mhi == NULL)
 745           return NULL;
 746         }
 747       }
 748     }
 749   else
 750     for(i = 1;; i++) {
 751       *s++ = (char)(dig = (int)(quorem(b,S) + '0'));
 752       if (!b->x[0] && b->wds <= 1) {
 753 #ifdef SET_INEXACT
 754         inexact = 0;
 755 #endif
 756         goto ret;
 757         }
 758       if (i >= ilim)
 759         break;
 760       b = multadd(b, 10, 0);
 761       if (b == NULL)
 762         return NULL;
 763       }
 764
 765   /* Round off last digit */
 766
 767 #ifdef Honor_FLT_ROUNDS
 768   switch(rounding) {
 769     case 0: goto trimzeros;
 770     case 2: goto roundoff;
 771     }
 772 #endif
 773   b = lshift(b, 1);
 774   j = cmp(b, S);
 775   if (j > 0 || (j == 0 && dig & 1)) {
 776  roundoff:
 777     while(*--s == '9')
 778       if (s == s0) {
 779         k++;
 780         *s++ = '1';
 781         goto ret;
 782         }
 783     ++*s++;
 784     }
 785   else {
 786 #ifdef Honor_FLT_ROUNDS
 787  trimzeros:
 788 #endif
 789     while(*--s == '0');
 790     s++;
 791     }
 792  ret:
 793   Bfree(S);
 794   if (mhi) {
 795     if (mlo && mlo != mhi)
 796       Bfree(mlo);
 797     Bfree(mhi);
 798     }
 799  ret1:
 800 #ifdef SET_INEXACT
 801   if (inexact) {
 802     if (!oldinexact) {
 803       word0(d) = Exp_1 + (70 << Exp_shift);
 804       word1(d) = 0;
 805       dval(d) += 1.;
 806       }
 807     }
 808   else if (!oldinexact)
 809     clear_inexact();
 810 #endif
 811   Bfree(b);
 812   if (s == s0) {      /* don't return empty string */
 813     *s++ = '0';
 814     k = 0;
 815   }
 816   *s = 0;
 817   *decpt = k + 1;
 818   if (rve)
 819     *rve = s;
 820   return s0;
 821   }