[Gmp-commit] /home/hgfiles/gmp: Handle divappr approximation problem more eff...

[mercurial at gmplib.org]

details:   /home/hgfiles/gmp/rev/79782b5e74c8
changeset: 13303:79782b5e74c8
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Sat Jan 02 23:57:51 2010 +0100
description:
Handle divappr approximation problem more efficiently.

diffstat:

 ChangeLog                 |   4 ++++
 mpn/generic/dcpi1_div_q.c |  19 ++++++++++++++++---
 mpn/generic/mu_div_q.c    |  43 +++++++++++++++++++++++++++++++------------
 3 files changed, 51 insertions(+), 15 deletions(-)

diffs (123 lines):

diff -r 31029f8dc100 -r 79782b5e74c8 ChangeLog
--- a/ChangeLog	Sat Jan 02 23:05:33 2010 +0100
+++ b/ChangeLog	Sat Jan 02 23:57:51 2010 +0100
@@ -5,6 +5,10 @@
 
 2010-01-02  Marco Bodrato <bodrato at mail.dm.unipi.it>
 
+	* mpn/generic/dcpi1_div_q.c: Handle divappr approximation problem more
+	efficiently.
+	* mpn/generic/mu_div_q.c: Likewise.
+
 	* mpn/generic/invert.c: Remove duplicated code.
 
 2010-01-01  Torbjorn Granlund  <tege at gmplib.org>
diff -r 31029f8dc100 -r 79782b5e74c8 mpn/generic/dcpi1_div_q.c
--- a/mpn/generic/dcpi1_div_q.c	Sat Jan 02 23:05:33 2010 +0100
+++ b/mpn/generic/dcpi1_div_q.c	Sat Jan 02 23:57:51 2010 +0100
@@ -1,7 +1,7 @@
 /* mpn_dc_div_q -- divide-and-conquer division, returning exact quotient
    only.
 
-   Contributed to the GNU project by Torbjorn Granlund.
+   Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato.
 
    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
@@ -53,8 +53,21 @@
   qh = mpn_dcpi1_divappr_q (wp, tp, nn + 1, dp, dn, dinv);
 
   if (wp[0] == 0)
-    /* FIXME: Should multiply and subtract here, not recompute from scratch.  */
-    qh = mpn_dcpi1_div_qr (qp, np, nn, dp, dn, dinv);
+    {
+      mp_limb_t cy;
+
+      if (qn > dn)
+	mpn_mul (tp, wp + 1, qn, dp, dn);
+      else
+	mpn_mul (tp, dp, dn, wp + 1, qn);
+
+      cy = (qh != 0) ? mpn_add_n (tp + qn, tp + qn, dp, dn) : 0;
+
+      if (cy || mpn_cmp (tp, np, nn) > 0) /* At most is wrong by one, no cycle. */
+	qh -= mpn_sub_1 (qp, wp + 1, qn, 1);
+      else /* Same as below */
+	MPN_COPY (qp, wp + 1, qn);
+    }
   else
     MPN_COPY (qp, wp + 1, qn);
 
diff -r 31029f8dc100 -r 79782b5e74c8 mpn/generic/mu_div_q.c
--- a/mpn/generic/mu_div_q.c	Sat Jan 02 23:05:33 2010 +0100
+++ b/mpn/generic/mu_div_q.c	Sat Jan 02 23:57:51 2010 +0100
@@ -1,6 +1,6 @@
 /* mpn_mu_div_q.
 
-   Contributed to the GNU project by Torbjorn Granlund.
+   Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato.
 
    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
@@ -114,13 +114,10 @@
 			 |________|   divisor  */
       rp = TMP_BALLOC_LIMBS (2 * dn + 1);
 
-	{
-	  this_in = mpn_mu_div_qr_choose_in (qn - dn, dn, 0);
-	  this_ip = ip + in - this_in;
-	  qh = mpn_preinv_mu_div_qr (tp + dn + 1, rp + dn + 1, np + dn, qn, dp, dn,
-				this_ip, this_in, scratch);
-	}
-
+      this_in = mpn_mu_div_qr_choose_in (qn - dn, dn, 0);
+      this_ip = ip + in - this_in;
+      qh = mpn_preinv_mu_div_qr (tp + dn + 1, rp + dn + 1, np + dn, qn, dp, dn,
+				 this_ip, this_in, scratch);
 
       MPN_COPY (rp + 1, np, dn);
       rp[0] = 0;
@@ -136,8 +133,19 @@
 	}
       else
 	{
-	  /* Fall back to plain mpn_mu_div_qr.  */
-	  mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
+	  mp_limb_t cy;
+	  mp_ptr pp;
+
+	  /* FIXME: can we use already allocated space? */
+	  pp = TMP_BALLOC_LIMBS (nn);
+	  mpn_mul (pp, tp + 1, qn, dp, dn);
+
+	  cy = (qh != 0) ? mpn_add_n (pp + qn, pp + qn, dp, dn) : 0;
+
+	  if (cy || mpn_cmp (pp, np, nn) > 0) /* At most is wrong by one, no cycle. */
+	    qh -= mpn_sub_1 (qp, tp + 1, qn, 1);
+	  else /* Same as above */
+	    MPN_COPY (qp, tp + 1, qn);
 	}
     }
   else
@@ -153,8 +161,19 @@
 	}
       else
 	{
-	  rp = TMP_BALLOC_LIMBS (dn);
-	  qh = mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
+	  mp_limb_t cy;
+
+	  /* FIXME: a shorter product should be enough; we may use already
+	     allocated space... */
+	  rp = TMP_BALLOC_LIMBS (nn);
+	  mpn_mul (rp, dp, dn, tp + 1, qn);
+
+	  cy = (qh != 0) ? mpn_add_n (rp + qn, rp + qn, dp, dn) : 0;
+
+	  if (cy || mpn_cmp (rp, np, nn) > 0) /* At most is wrong by one, no cycle. */
+	    qh -= mpn_sub_1 (qp, tp + 1, qn, 1);
+	  else /* Same as above */
+	    MPN_COPY (qp, tp + 1, qn);
 	}
     }