newlib/libc/stdlib/dtoa.c

   1 /****************************************************************
   2  *
   3  * The author of this software is David M. Gay.
   4  *
   5  * Copyright (c) 1991 by AT&T.
   6  *
   7  * Permission to use, copy, modify, and distribute this software for any
   8  * purpose without fee is hereby granted, provided that this entire notice
   9  * is included in all copies of any software which is or includes a copy
  10  * or modification of this software and in all copies of the supporting
  11  * documentation for such software.
  12  *
  13  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  15  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17  *
  18  ***************************************************************/
  19
  20 /* Please send bug reports to
  21         David M. Gay
  22         AT&T Bell Laboratories, Room 2C-463
  23         600 Mountain Avenue
  24         Murray Hill, NJ 07974-2070
  25         U.S.A.
  26         dmg@research.att.com or research!dmg
  27  */
  28
  29 #include <_ansi.h>
  30 #include <stdlib.h>
  31 #include <reent.h>
  32 #include <string.h>
  33 #include "mprec.h"
  34
  35 #ifdef _REENT_THREAD_LOCAL
  36 _Thread_local struct _Bigint *_tls_mp_result;
  37 _Thread_local int _tls_mp_result_k;
  38 #endif
  39
  40 static int
  41 quorem (_Bigint * b, _Bigint * S)
  42 {
  43   int n;
  44   __Long borrow, y;
  45   __ULong carry, q, ys;
  46   __ULong *bx, *bxe, *sx, *sxe;
  47 #ifdef Pack_32
  48   __Long z;
  49   __ULong si, zs;
  50 #endif
  51
  52   n = S->_wds;
  53 #ifdef DEBUG
  54   /*debug*/ if (b->_wds > n)
  55     /*debug*/ Bug ("oversize b in quorem");
  56 #endif
  57   if (b->_wds < n)
  58     return 0;
  59   sx = S->_x;
  60   sxe = sx + --n;
  61   bx = b->_x;
  62   bxe = bx + n;
  63   q = *bxe / (*sxe + 1);        /* ensure q <= true quotient */
  64 #ifdef DEBUG
  65   /*debug*/ if (q > 9)
  66     /*debug*/ Bug ("oversized quotient in quorem");
  67 #endif
  68   if (q)
  69     {
  70       borrow = 0;
  71       carry = 0;
  72       do
  73         {
  74 #ifdef Pack_32
  75           si = *sx++;
  76           ys = (si & 0xffff) * q + carry;
  77           zs = (si >> 16) * q + (ys >> 16);
  78           carry = zs >> 16;
  79           y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  80           borrow = y >> 16;
  81           Sign_Extend (borrow, y);
  82           z = (*bx >> 16) - (zs & 0xffff) + borrow;
  83           borrow = z >> 16;
  84           Sign_Extend (borrow, z);
  85           Storeinc (bx, z, y);
  86 #else
  87           ys = *sx++ * q + carry;
  88           carry = ys >> 16;
  89           y = *bx - (ys & 0xffff) + borrow;
  90           borrow = y >> 16;
  91           Sign_Extend (borrow, y);
  92           *bx++ = y & 0xffff;
  93 #endif
  94         }
  95       while (sx <= sxe);
  96       if (!*bxe)
  97         {
  98           bx = b->_x;
  99           while (--bxe > bx && !*bxe)
 100             --n;
 101           b->_wds = n;
 102         }
 103     }
 104   if (cmp (b, S) >= 0)
 105     {
 106       q++;
 107       borrow = 0;
 108       carry = 0;
 109       bx = b->_x;
 110       sx = S->_x;
 111       do
 112         {
 113 #ifdef Pack_32
 114           si = *sx++;
 115           ys = (si & 0xffff) + carry;
 116           zs = (si >> 16) + (ys >> 16);
 117           carry = zs >> 16;
 118           y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
 119           borrow = y >> 16;
 120           Sign_Extend (borrow, y);
 121           z = (*bx >> 16) - (zs & 0xffff) + borrow;
 122           borrow = z >> 16;
 123           Sign_Extend (borrow, z);
 124           Storeinc (bx, z, y);
 125 #else
 126           ys = *sx++ + carry;
 127           carry = ys >> 16;
 128           y = *bx - (ys & 0xffff) + borrow;
 129           borrow = y >> 16;
 130           Sign_Extend (borrow, y);
 131           *bx++ = y & 0xffff;
 132 #endif
 133         }
 134       while (sx <= sxe);
 135       bx = b->_x;
 136       bxe = bx + n;
 137       if (!*bxe)
 138         {
 139           while (--bxe > bx && !*bxe)
 140             --n;
 141           b->_wds = n;
 142         }
 143     }
 144   return q;
 145 }
 146
 147 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
 148  *
 149  * Inspired by "How to Print Floating-Point Numbers Accurately" by
 150  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
 151  *
 152  * Modifications:
 153  *      1. Rather than iterating, we use a simple numeric overestimate
 154  *         to determine k = floor(log10(d)).  We scale relevant
 155  *         quantities using O(log2(k)) rather than O(k) multiplications.
 156  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
 157  *         try to generate digits strictly left to right.  Instead, we
 158  *         compute with fewer bits and propagate the carry if necessary
 159  *         when rounding the final digit up.  This is often faster.
 160  *      3. Under the assumption that input will be rounded nearest,
 161  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
 162  *         That is, we allow equality in stopping tests when the
 163  *         round-nearest rule will give the same floating-point value
 164  *         as would satisfaction of the stopping test with strict
 165  *         inequality.
 166  *      4. We remove common factors of powers of 2 from relevant
 167  *         quantities.
 168  *      5. When converting floating-point integers less than 1e16,
 169  *         we use floating-point arithmetic rather than resorting
 170  *         to multiple-precision integers.
 171  *      6. When asked to produce fewer than 15 digits, we first try
 172  *         to get by with floating-point arithmetic; we resort to
 173  *         multiple-precision integer arithmetic only if we cannot
 174  *         guarantee that the floating-point calculation has given
 175  *         the correctly rounded result.  For k requested digits and
 176  *         "uniformly" distributed input, the probability is
 177  *         something like 10^(k-15) that we must resort to the long
 178  *         calculation.
 179  */
 180
 181
 182 char *
 183 _dtoa_r (struct _reent *ptr,
 184         double _d,
 185         int mode,
 186         int ndigits,
 187         int *decpt,
 188         int *sign,
 189         char **rve)
 190 {
 191   /*    Arguments ndigits, decpt, sign are similar to those
 192         of ecvt and fcvt; trailing zeros are suppressed from
 193         the returned string.  If not null, *rve is set to point
 194         to the end of the return value.  If d is +-Infinity or NaN,
 195         then *decpt is set to 9999.
 196
 197         mode:
 198                 0 ==> shortest string that yields d when read in
 199                         and rounded to nearest.
 200                 1 ==> like 0, but with Steele & White stopping rule;
 201                         e.g. with IEEE P754 arithmetic , mode 0 gives
 202                         1e23 whereas mode 1 gives 9.999999999999999e22.
 203                 2 ==> max(1,ndigits) significant digits.  This gives a
 204                         return value similar to that of ecvt, except
 205                         that trailing zeros are suppressed.
 206                 3 ==> through ndigits past the decimal point.  This
 207                         gives a return value similar to that from fcvt,
 208                         except that trailing zeros are suppressed, and
 209                         ndigits can be negative.
 210                 4-9 should give the same return values as 2-3, i.e.,
 211                         4 <= mode <= 9 ==> same return as mode
 212                         2 + (mode & 1).  These modes are mainly for
 213                         debugging; often they run slower but sometimes
 214                         faster than modes 2-3.
 215                 4,5,8,9 ==> left-to-right digit generation.
 216                 6-9 ==> don't try fast floating-point estimate
 217                         (if applicable).
 218
 219                 Values of mode other than 0-9 are treated as mode 0.
 220
 221                 Sufficient space is allocated to the return value
 222                 to hold the suppressed trailing zeros.
 223         */
 224
 225   int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, j, j1, k, k0,
 226     k_check, leftright, m2, m5, s2, s5, spec_case, try_quick;
 227   union double_union d, d2, eps;
 228   __Long L;
 229 #ifndef Sudden_Underflow
 230   int denorm;
 231   __ULong x;
 232 #endif
 233   _Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
 234   double ds;
 235   char *s, *s0;
 236
 237   d.d = _d;
 238
 239   _REENT_CHECK_MP(ptr);
 240   if (_REENT_MP_RESULT(ptr))
 241     {
 242       _REENT_MP_RESULT(ptr)->_k = _REENT_MP_RESULT_K(ptr);
 243       _REENT_MP_RESULT(ptr)->_maxwds = 1 << _REENT_MP_RESULT_K(ptr);
 244       Bfree (ptr, _REENT_MP_RESULT(ptr));
 245       _REENT_MP_RESULT(ptr) = 0;
 246     }
 247
 248   if (word0 (d) & Sign_bit)
 249     {
 250       /* set sign for everything, including 0's and NaNs */
 251       *sign = 1;
 252       word0 (d) &= ~Sign_bit;   /* clear sign bit */
 253     }
 254   else
 255     *sign = 0;
 256
 257 #if defined(IEEE_Arith) + defined(VAX)
 258 #ifdef IEEE_Arith
 259   if ((word0 (d) & Exp_mask) == Exp_mask)
 260 #else
 261   if (word0 (d) == 0x8000)
 262 #endif
 263     {
 264       /* Infinity or NaN */
 265       *decpt = 9999;
 266       s =
 267 #ifdef IEEE_Arith
 268         !word1 (d) && !(word0 (d) & 0xfffff) ? "Infinity" :
 269 #endif
 270         "NaN";
 271       if (rve)
 272         *rve =
 273 #ifdef IEEE_Arith
 274           s[3] ? s + 8 :
 275 #endif
 276           s + 3;
 277       return s;
 278     }
 279 #endif
 280 #ifdef IBM
 281   d.d += 0;                     /* normalize */
 282 #endif
 283   if (!d.d)
 284     {
 285       *decpt = 1;
 286       s = "0";
 287       if (rve)
 288         *rve = s + 1;
 289       return s;
 290     }
 291
 292   b = d2b (ptr, d.d, &be, &bbits);
 293 #ifdef Sudden_Underflow
 294   i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
 295 #else
 296   if ((i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))) != 0)
 297     {
 298 #endif
 299       d2.d = d.d;
 300       word0 (d2) &= Frac_mask1;
 301       word0 (d2) |= Exp_11;
 302 #ifdef IBM
 303       if (j = 11 - hi0bits (word0 (d2) & Frac_mask))
 304         d2.d /= 1 << j;
 305 #endif
 306
 307       /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
 308                  * log10(x)      =  log(x) / log(10)
 309                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 310                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
 311                  *
 312                  * This suggests computing an approximation k to log10(d) by
 313                  *
 314                  * k = (i - Bias)*0.301029995663981
 315                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 316                  *
 317                  * We want k to be too large rather than too small.
 318                  * The error in the first-order Taylor series approximation
 319                  * is in our favor, so we just round up the constant enough
 320                  * to compensate for any error in the multiplication of
 321                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 322                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 323                  * adding 1e-13 to the constant term more than suffices.
 324                  * Hence we adjust the constant term to 0.1760912590558.
 325                  * (We could get a more accurate k by invoking log10,
 326                  *  but this is probably not worthwhile.)
 327                  */
 328
 329       i -= Bias;
 330 #ifdef IBM
 331       i <<= 2;
 332       i += j;
 333 #endif
 334 #ifndef Sudden_Underflow
 335       denorm = 0;
 336     }
 337   else
 338     {
 339       /* d is denormalized */
 340
 341       i = bbits + be + (Bias + (P - 1) - 1);
 342 #if defined (_DOUBLE_IS_32BITS)
 343       x = word0 (d) << (32 - i);
 344 #else
 345       x = (i > 32) ? (word0 (d) << (64 - i)) | (word1 (d) >> (i - 32))
 346        : (word1 (d) << (32 - i));
 347 #endif
 348       d2.d = x;
 349       word0 (d2) -= 31 * Exp_msk1;      /* adjust exponent */
 350       i -= (Bias + (P - 1) - 1) + 1;
 351       denorm = 1;
 352     }
 353 #endif
 354 #if defined (_DOUBLE_IS_32BITS)
 355   ds = (d2.d - 1.5) * 0.289529651 + 0.176091269 + i * 0.30103001;
 356 #else
 357   ds = (d2.d - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
 358 #endif
 359   k = (int) ds;
 360   if (ds < 0. && ds != k)
 361     k--;                        /* want k = floor(ds) */
 362   k_check = 1;
 363   if (k >= 0 && k <= Ten_pmax)
 364     {
 365       if (d.d < tens[k])
 366         k--;
 367       k_check = 0;
 368     }
 369   j = bbits - i - 1;
 370   if (j >= 0)
 371     {
 372       b2 = 0;
 373       s2 = j;
 374     }
 375   else
 376     {
 377       b2 = -j;
 378       s2 = 0;
 379     }
 380   if (k >= 0)
 381     {
 382       b5 = 0;
 383       s5 = k;
 384       s2 += k;
 385     }
 386   else
 387     {
 388       b2 -= k;
 389       b5 = -k;
 390       s5 = 0;
 391     }
 392   if (mode < 0 || mode > 9)
 393     mode = 0;
 394   try_quick = 1;
 395   if (mode > 5)
 396     {
 397       mode -= 4;
 398       try_quick = 0;
 399     }
 400   leftright = 1;
 401   ilim = ilim1 = -1;
 402   switch (mode)
 403     {
 404     case 0:
 405     case 1:
 406       i = 18;
 407       ndigits = 0;
 408       break;
 409     case 2:
 410       leftright = 0;
 411       /* no break */
 412     case 4:
 413       if (ndigits <= 0)
 414         ndigits = 1;
 415       ilim = ilim1 = i = ndigits;
 416       break;
 417     case 3:
 418       leftright = 0;
 419       /* no break */
 420     case 5:
 421       i = ndigits + k + 1;
 422       ilim = i;
 423       ilim1 = i - 1;
 424       if (i <= 0)
 425         i = 1;
 426     }
 427   j = sizeof (__ULong);
 428   for (_REENT_MP_RESULT_K(ptr) = 0; sizeof (_Bigint) - sizeof (__ULong) + j <= i;
 429        j <<= 1)
 430     _REENT_MP_RESULT_K(ptr)++;
 431   _REENT_MP_RESULT(ptr) = eBalloc (ptr, _REENT_MP_RESULT_K(ptr));
 432   s = s0 = (char *) _REENT_MP_RESULT(ptr);
 433
 434   if (ilim >= 0 && ilim <= Quick_max && try_quick)
 435     {
 436       /* Try to get by with floating-point arithmetic. */
 437
 438       i = 0;
 439       d2.d = d.d;
 440       k0 = k;
 441       ilim0 = ilim;
 442       ieps = 2;                 /* conservative */
 443       if (k > 0)
 444         {
 445           ds = tens[k & 0xf];
 446           j = k >> 4;
 447           if (j & Bletch)
 448             {
 449               /* prevent overflows */
 450               j &= Bletch - 1;
 451               d.d /= bigtens[n_bigtens - 1];
 452               ieps++;
 453             }
 454           for (; j; j >>= 1, i++)
 455             if (j & 1)
 456               {
 457                 ieps++;
 458                 ds *= bigtens[i];
 459               }
 460           d.d /= ds;
 461         }
 462       else if ((j1 = -k) != 0)
 463         {
 464           d.d *= tens[j1 & 0xf];
 465           for (j = j1 >> 4; j; j >>= 1, i++)
 466             if (j & 1)
 467               {
 468                 ieps++;
 469                 d.d *= bigtens[i];
 470               }
 471         }
 472       if (k_check && d.d < 1. && ilim > 0)
 473         {
 474           if (ilim1 <= 0)
 475             goto fast_failed;
 476           ilim = ilim1;
 477           k--;
 478           d.d *= 10.;
 479           ieps++;
 480         }
 481       eps.d = ieps * d.d + 7.;
 482       word0 (eps) -= (P - 1) * Exp_msk1;
 483       if (ilim == 0)
 484         {
 485           S = mhi = 0;
 486           d.d -= 5.;
 487           if (d.d > eps.d)
 488             goto one_digit;
 489           if (d.d < -eps.d)
 490             goto no_digits;
 491           goto fast_failed;
 492         }
 493 #ifndef No_leftright
 494       if (leftright)
 495         {
 496           /* Use Steele & White method of only
 497            * generating digits needed.
 498            */
 499           eps.d = 0.5 / tens[ilim - 1] - eps.d;
 500           for (i = 0;;)
 501             {
 502               L = d.d;
 503               d.d -= L;
 504               *s++ = '0' + (int) L;
 505               if (d.d < eps.d)
 506                 goto ret1;
 507               if (1. - d.d < eps.d)
 508                 goto bump_up;
 509               if (++i >= ilim)
 510                 break;
 511               eps.d *= 10.;
 512               d.d *= 10.;
 513             }
 514         }
 515       else
 516         {
 517 #endif
 518           /* Generate ilim digits, then fix them up. */
 519           eps.d *= tens[ilim - 1];
 520           for (i = 1;; i++, d.d *= 10.)
 521             {
 522               L = d.d;
 523               d.d -= L;
 524               *s++ = '0' + (int) L;
 525               if (i == ilim)
 526                 {
 527                   if (d.d > 0.5 + eps.d)
 528                     goto bump_up;
 529                   else if (d.d < 0.5 - eps.d)
 530                     {
 531                       while (*--s == '0');
 532                       s++;
 533                       goto ret1;
 534                     }
 535                   break;
 536                 }
 537             }
 538 #ifndef No_leftright
 539         }
 540 #endif
 541     fast_failed:
 542       s = s0;
 543       d.d = d2.d;
 544       k = k0;
 545       ilim = ilim0;
 546     }
 547
 548   /* Do we have a "small" integer? */
 549
 550   if (be >= 0 && k <= Int_max)
 551     {
 552       /* Yes. */
 553       ds = tens[k];
 554       if (ndigits < 0 && ilim <= 0)
 555         {
 556           S = mhi = 0;
 557           if (ilim < 0 || d.d <= 5 * ds)
 558             goto no_digits;
 559           goto one_digit;
 560         }
 561       for (i = 1;; i++)
 562         {
 563           L = d.d / ds;
 564           d.d -= L * ds;
 565 #ifdef Check_FLT_ROUNDS
 566           /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 567           if (d.d < 0)
 568             {
 569               L--;
 570               d.d += ds;
 571             }
 572 #endif
 573           *s++ = '0' + (int) L;
 574           if (i == ilim)
 575             {
 576               d.d += d.d;
 577              if ((d.d > ds) || ((d.d == ds) && (L & 1)))
 578                 {
 579                 bump_up:
 580                   while (*--s == '9')
 581                     if (s == s0)
 582                       {
 583                         k++;
 584                         *s = '0';
 585                         break;
 586                       }
 587                   ++*s++;
 588                 }
 589               break;
 590             }
 591           if (!(d.d *= 10.))
 592             break;
 593         }
 594       goto ret1;
 595     }
 596
 597   m2 = b2;
 598   m5 = b5;
 599   mhi = mlo = 0;
 600   if (leftright)
 601     {
 602       if (mode < 2)
 603         {
 604           i =
 605 #ifndef Sudden_Underflow
 606             denorm ? be + (Bias + (P - 1) - 1 + 1) :
 607 #endif
 608 #ifdef IBM
 609             1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3);
 610 #else
 611             1 + P - bbits;
 612 #endif
 613         }
 614       else
 615         {
 616           j = ilim - 1;
 617           if (m5 >= j)
 618             m5 -= j;
 619           else
 620             {
 621               s5 += j -= m5;
 622               b5 += j;
 623               m5 = 0;
 624             }
 625           if ((i = ilim) < 0)
 626             {
 627               m2 -= i;
 628               i = 0;
 629             }
 630         }
 631       b2 += i;
 632       s2 += i;
 633       mhi = i2b (ptr, 1);
 634     }
 635   if (m2 > 0 && s2 > 0)
 636     {
 637       i = m2 < s2 ? m2 : s2;
 638       b2 -= i;
 639       m2 -= i;
 640       s2 -= i;
 641     }
 642   if (b5 > 0)
 643     {
 644       if (leftright)
 645         {
 646           if (m5 > 0)
 647             {
 648               mhi = pow5mult (ptr, mhi, m5);
 649               b1 = mult (ptr, mhi, b);
 650               Bfree (ptr, b);
 651               b = b1;
 652             }
 653          if ((j = b5 - m5) != 0)
 654             b = pow5mult (ptr, b, j);
 655         }
 656       else
 657         b = pow5mult (ptr, b, b5);
 658     }
 659   S = i2b (ptr, 1);
 660   if (s5 > 0)
 661     S = pow5mult (ptr, S, s5);
 662
 663   /* Check for special case that d is a normalized power of 2. */
 664
 665   spec_case = 0;
 666   if (mode < 2)
 667     {
 668       if (!word1 (d) && !(word0 (d) & Bndry_mask)
 669 #ifndef Sudden_Underflow
 670           && word0 (d) & Exp_mask
 671 #endif
 672         )
 673         {
 674           /* The special case */
 675           b2 += Log2P;
 676           s2 += Log2P;
 677           spec_case = 1;
 678         }
 679     }
 680
 681   /* Arrange for convenient computation of quotients:
 682    * shift left if necessary so divisor has 4 leading 0 bits.
 683    *
 684    * Perhaps we should just compute leading 28 bits of S once
 685    * and for all and pass them and a shift to quorem, so it
 686    * can do shifts and ors to compute the numerator for q.
 687    */
 688
 689 #ifdef Pack_32
 690   if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0x1f) != 0)
 691     i = 32 - i;
 692 #else
 693   if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0xf) != 0)
 694     i = 16 - i;
 695 #endif
 696   if (i > 4)
 697     {
 698       i -= 4;
 699       b2 += i;
 700       m2 += i;
 701       s2 += i;
 702     }
 703   else if (i < 4)
 704     {
 705       i += 28;
 706       b2 += i;
 707       m2 += i;
 708       s2 += i;
 709     }
 710   if (b2 > 0)
 711     b = lshift (ptr, b, b2);
 712   if (s2 > 0)
 713     S = lshift (ptr, S, s2);
 714   if (k_check)
 715     {
 716       if (cmp (b, S) < 0)
 717         {
 718           k--;
 719           b = multadd (ptr, b, 10, 0);  /* we botched the k estimate */
 720           if (leftright)
 721             mhi = multadd (ptr, mhi, 10, 0);
 722           ilim = ilim1;
 723         }
 724     }
 725   if (ilim <= 0 && mode > 2)
 726     {
 727       if (ilim < 0 || cmp (b, S = multadd (ptr, S, 5, 0)) <= 0)
 728         {
 729           /* no digits, fcvt style */
 730         no_digits:
 731           k = -1 - ndigits;
 732           goto ret;
 733         }
 734     one_digit:
 735       *s++ = '1';
 736       k++;
 737       goto ret;
 738     }
 739   if (leftright)
 740     {
 741       if (m2 > 0)
 742         mhi = lshift (ptr, mhi, m2);
 743
 744       /* Compute mlo -- check for special case
 745        * that d is a normalized power of 2.
 746        */
 747
 748       mlo = mhi;
 749       if (spec_case)
 750         {
 751           mhi = eBalloc (ptr, mhi->_k);
 752           Bcopy (mhi, mlo);
 753           mhi = lshift (ptr, mhi, Log2P);
 754         }
 755
 756       for (i = 1;; i++)
 757         {
 758           dig = quorem (b, S) + '0';
 759           /* Do we yet have the shortest decimal string
 760            * that will round to d?
 761            */
 762           j = cmp (b, mlo);
 763           delta = diff (ptr, S, mhi);
 764           j1 = delta->_sign ? 1 : cmp (b, delta);
 765           Bfree (ptr, delta);
 766 #ifndef ROUND_BIASED
 767           if (j1 == 0 && !mode && !(word1 (d) & 1))
 768             {
 769               if (dig == '9')
 770                 goto round_9_up;
 771               if (j > 0)
 772                 dig++;
 773               *s++ = dig;
 774               goto ret;
 775             }
 776 #endif
 777          if ((j < 0) || ((j == 0) && !mode
 778 #ifndef ROUND_BIASED
 779               && !(word1 (d) & 1)
 780 #endif
 781            ))
 782             {
 783               if (j1 > 0)
 784                 {
 785                   b = lshift (ptr, b, 1);
 786                   j1 = cmp (b, S);
 787                  if (((j1 > 0) || ((j1 == 0) && (dig & 1)))
 788                       && dig++ == '9')
 789                     goto round_9_up;
 790                 }
 791               *s++ = dig;
 792               goto ret;
 793             }
 794           if (j1 > 0)
 795             {
 796               if (dig == '9')
 797                 {               /* possible if i == 1 */
 798                 round_9_up:
 799                   *s++ = '9';
 800                   goto roundoff;
 801                 }
 802               *s++ = dig + 1;
 803               goto ret;
 804             }
 805           *s++ = dig;
 806           if (i == ilim)
 807             break;
 808           b = multadd (ptr, b, 10, 0);
 809           if (mlo == mhi)
 810             mlo = mhi = multadd (ptr, mhi, 10, 0);
 811           else
 812             {
 813               mlo = multadd (ptr, mlo, 10, 0);
 814               mhi = multadd (ptr, mhi, 10, 0);
 815             }
 816         }
 817     }
 818   else
 819     for (i = 1;; i++)
 820       {
 821         *s++ = dig = quorem (b, S) + '0';
 822         if (i >= ilim)
 823           break;
 824         b = multadd (ptr, b, 10, 0);
 825       }
 826
 827   /* Round off last digit */
 828
 829   b = lshift (ptr, b, 1);
 830   j = cmp (b, S);
 831   if ((j > 0) || ((j == 0) && (dig & 1)))
 832     {
 833     roundoff:
 834       while (*--s == '9')
 835         if (s == s0)
 836           {
 837             k++;
 838             *s++ = '1';
 839             goto ret;
 840           }
 841       ++*s++;
 842     }
 843   else
 844     {
 845       while (*--s == '0');
 846       s++;
 847     }
 848 ret:
 849   Bfree (ptr, S);
 850   if (mhi)
 851     {
 852       if (mlo && mlo != mhi)
 853         Bfree (ptr, mlo);
 854       Bfree (ptr, mhi);
 855     }
 856 ret1:
 857   Bfree (ptr, b);
 858   *s = 0;
 859   *decpt = k + 1;
 860   if (rve)
 861     *rve = s;
 862   return s0;
 863 }