* mpf/set_q.c: Use mpn_tdiv_qr rather than mpn_divrem, so no shifting.

	Don't truncate the divisor, it can make the result inaccurate.

1 files changed, 68 insertions, 79 deletions

diff --git a/mpf/set_q.c b/mpf/set_q.c
index d13b42a63..075f3c3c0 100644
--- a/mpf/set_q.c
+++ b/mpf/set_q.c
@@ -23,17 +23,50 @@ MA 02111-1307, USA. */
 #include "gmp-impl.h"
 #include "longlong.h"
 
+
+/* As usual the aim is to produce PREC(r) limbs, with the high non-zero.
+   The basic mpn_tdiv_qr produces a quotient of nsize-dsize+1 limbs, with
+   either the high or second highest limb non-zero.  We arrange for
+   nsize-dsize+1 to equal prec+1, hence giving either prec or prec+1 result
+   limbs at PTR(r).
+
+   nsize-dsize+1 == prec+1 is achieved by adjusting num(q), either dropping
+   low limbs if it's too big, or padding with low zeros if it's too small.
+   The full given den(q) is always used.
+
+   We cannot truncate den(q), because even when it's much bigger than prec
+   the last limbs can still influence the final quotient.  Often they don't,
+   but we leave optimization of that to a prospective quotient-only mpn
+   division.
+
+   Not done:
+
+   If den(q) is a power of 2 then we may end up with low zero limbs on the
+   result.  But nothing is done about this, since it should be unlikely on
+   random data, and can be left to an application to call mpf_div_2exp if it
+   might occur with any frequency.
+
+   Enhancements:
+
+   The high quotient limb is non-zero when high{np,dsize} >= {dp,dsize}.  We
+   could make that comparison and use qsize==prec instead of qsize==prec+1,
+   to save one limb in the division.
+
+   Future:
+
+   If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of
+   padding n with zeros in temporary space.
+
+   If/when a quotient-only division exists it can be used here immediately.
+   remp is only to satisfy mpn_tdiv_qr, the remainder is not used.  */
+
 void
 mpf_set_q (mpf_t r, mpq_srcptr q)
 {
-  mp_ptr np, dp;
-  mp_ptr rp;
-  mp_size_t nsize, dsize;
-  mp_size_t qsize, rsize;
-  mp_size_t sign_quotient;
-  mp_limb_t qlimb;
-  mp_ptr qp;
-  mp_size_t prec;
+  mp_srcptr np, dp;
+  mp_size_t prec, nsize, dsize, qsize, prospective_qsize, tsize, zeros;
+  mp_size_t sign_quotient, high_zero;
+  mp_ptr qp, tp, remp;
   mp_exp_t exp;
   TMP_DECL (marker);
 
@@ -42,17 +75,16 @@ mpf_set_q (mpf_t r, mpq_srcptr q)
   nsize = SIZ (&q->_mp_num);
   dsize = SIZ (&q->_mp_den);
 
-  if (nsize == 0)
+  if (UNLIKELY (nsize == 0))
     {
       SIZ (r) = 0;
       EXP (r) = 0;
       return;
     }
 
-  prec = PREC (r) + 1;
-
   TMP_MARK (marker);
 
+  prec = PREC (r);
   qp = PTR (r);
 
   sign_quotient = nsize;
@@ -60,83 +92,40 @@ mpf_set_q (mpf_t r, mpq_srcptr q)
   np = PTR (&q->_mp_num);
   dp = PTR (&q->_mp_den);
 
-  exp = nsize - dsize;
+  prospective_qsize = nsize - dsize + 1;  /* q from using given n,d sizes */
+  exp = prospective_qsize;                /* ie. number of integer limbs */
+  qsize = prec + 1;                       /* desired q */
 
-  if (nsize > prec)
-    {
-      np += nsize - prec;
-      nsize = prec;
-    }
-  if (dsize > prec)
-    {
-      dp += dsize - prec;
-      dsize = prec;
-    }
+  /* how many zero limbs to pad n with to get desired qsize */
+  zeros = qsize - prospective_qsize;
+
+  /* space for copy of n, if need zeros */
+  tsize = (zeros > 0 ? nsize + zeros : 1);
 
-  rsize = MAX (nsize, dsize);
-  rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB);
+  TMP_ALLOC_LIMBS_2 (remp, dsize, tp, tsize);
 
-  /* Normalize the denominator, i.e. make its most significant bit set by
-     shifting it NORMALIZATION_STEPS bits to the left.  Also shift the
-     numerator the same number of steps (to keep the quotient the same!).  */
-  if ((dp[dsize - 1] & GMP_NUMB_HIGHBIT) == 0)
+  if (zeros > 0)
     {
-      mp_ptr tp;
-      mp_limb_t nlimb;
-      unsigned normalization_steps;
-
-      count_leading_zeros (normalization_steps, dp[dsize - 1]);
-      normalization_steps -= GMP_NAIL_BITS;
-
-      /* Shift up the denominator setting the most significant bit of
-	 the most significant limb.  Use temporary storage not to clobber
-	 the original contents of the denominator.  */
-      tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
-      mpn_lshift (tp, dp, dsize, normalization_steps);
-      dp = tp;
-
-      if (rsize != nsize)
-	{
-	  MPN_ZERO (rp, rsize - nsize);
-	  nlimb = mpn_lshift (rp + (rsize - nsize),
-			      np, nsize, normalization_steps);
-	}
-      else
-	{
-	  nlimb = mpn_lshift (rp, np, rsize, normalization_steps);
-	}
-      if (nlimb != 0)
-	{
-	  rp[rsize] = nlimb;
-	  exp++;
-	  /* Don't just increase rsize, chop off rp at the low end instead.  */
-	  if (rsize == prec)
-	    rp++;
-	  else
-	    rsize++;
-	}
+      /* pad n with zeros into temporary space */
+      MPN_ZERO (tp, zeros);
+      MPN_COPY (tp+zeros, np, nsize);
+      np = tp;
+      nsize = tsize;
     }
   else
     {
-      if (rsize != nsize)
-	{
-	  MPN_ZERO (rp, rsize - nsize);
-	  MPN_COPY (rp + (rsize - nsize), np, nsize);
-	}
-      else
-	{
-	  MPN_COPY (rp, np, rsize);
-	}
+      /* shorten n to get desired qsize */
+      nsize += zeros;
+      np -= zeros;
     }
 
-  qlimb = mpn_divrem (qp, prec - 1 - (rsize - dsize), rp, rsize, dp, dsize);
-  qsize = prec - 1;
-  if (qlimb)
-    {
-      qp[qsize] = qlimb;
-      qsize++;
-      exp++;
-    }
+  ASSERT (nsize-dsize+1 == qsize);
+  mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize);
+
+  /* strip possible zero high limb */
+  high_zero = (qp[qsize-1] == 0);
+  qsize -= high_zero;
+  exp -= high_zero;
 
   EXP (r) = exp;
   SIZ (r) = sign_quotient >= 0 ? qsize : -qsize;