Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

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