libjava/classpath/native/fdlibm/dtoa.c

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