(mpn_gcd_subdiv_step): Reorganized to use a single hook function.

6 files changed, 139 insertions, 134 deletions

diff --git a/ChangeLog b/ChangeLog
index 22ed83a18..0f302fe9b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,17 @@
+2010-05-03  Niels Möller  <nisse@lysator.liu.se>
+
+	* mpn/generic/gcd_subdiv_step.c (mpn_gcd_subdiv_step): Reorganized
+	to use a single hook function.
+	* mpn/generic/gcdext.c (mpn_gcdext): Adapted to new hook
+	interface.
+	* mpn/generic/gcdext_lehmer.c (mpn_gcdext_hook): New unified hook
+	function.
+	* mpn/generic/gcd.c (gcd_hook): Renamed from gcd_done, and adapted
+	to new hook interface.
+	* gmp-impl.h (gcd_subdiv_step_hook): New typedef, now a function
+	type, not a struct.
+	(mpn_gcdext_hook): Declare.
+
 2010-05-03  Torbjorn Granlund  <tege@gmplib.org>
 
 	* mpn/powerpc64/sqr_diagonal.asm: Complete rewrite.
diff --git a/gmp-impl.h b/gmp-impl.h
index ee311f172..8561c7192 100644
--- a/gmp-impl.h
+++ b/gmp-impl.h
@@ -3859,18 +3859,13 @@ __GMP_DECLSPEC mp_size_t mpn_hgcd_itch __GMP_PROTO ((mp_size_t));
 #define mpn_hgcd __MPN (hgcd)
 __GMP_DECLSPEC mp_size_t mpn_hgcd __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, struct hgcd_matrix *, mp_ptr));
 
-struct gcd_subdiv_step_hook
-{
-  /* Passes quotient. */
-  void (*update)(void *, mp_srcptr, mp_size_t, unsigned);
-  /* Passes final gcd (must be copied if it is to be retained). */
-  void (*done)(void *, mp_srcptr, mp_size_t, unsigned);
-};
+typedef void gcd_subdiv_step_hook(void *, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, int);
+
 /* Needs storage for the quotient */
 #define MPN_GCD_SUBDIV_STEP_ITCH(n) (n)
 
 #define mpn_gcd_subdiv_step __MPN(gcd_subdiv_step)
-__GMP_DECLSPEC mp_size_t mpn_gcd_subdiv_step __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, const struct gcd_subdiv_step_hook *, void *, mp_ptr));
+__GMP_DECLSPEC mp_size_t mpn_gcd_subdiv_step __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, gcd_subdiv_step_hook *, void *, mp_ptr));
 
 struct gcdext_ctx
 {
@@ -3885,9 +3880,8 @@ struct gcdext_ctx
   mp_ptr u0, u1, tp;
 };
 
-#define gcdext_hook __gmp_gcdext_hook
-extern const struct gcd_subdiv_step_hook
-gcdext_hook;
+#define mpn_gcdext_hook __MPN (gcdext_hook)
+gcd_subdiv_step_hook mpn_gcdext_hook;
 
 #define MPN_GCDEXT_LEHMER_N_ITCH(n) (4*(n) + 3)
 
diff --git a/mpn/generic/gcd.c b/mpn/generic/gcd.c
index 2d1e41266..94ab21742 100644
--- a/mpn/generic/gcd.c
+++ b/mpn/generic/gcd.c
@@ -18,8 +18,6 @@ License for more details.
 You should have received a copy of the GNU Lesser General Public License
 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
 
-#include <stdlib.h>		/* for NULL */
-
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "longlong.h"
@@ -60,18 +58,13 @@ struct gcd_ctx
 };
 
 static void
-gcd_done (void *p, mp_srcptr gp, mp_size_t gn, unsigned swapped)
+gcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
+	  mp_srcptr qp, mp_size_t qn, int d)
 {
   struct gcd_ctx *ctx = (struct gcd_ctx *) p;
   MPN_COPY (ctx->gp, gp, gn);
   ctx->gn = gn;
 }
-
-static const struct gcd_subdiv_step_hook
-gcd_hook = {
-  NULL,
-  gcd_done,
-};
   
 #if GMP_NAIL_BITS > 0
 /* Nail supports should be easy, replacing the sub_ddmmss with nails
@@ -213,7 +206,7 @@ mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
       else
 	{
 	  /* Temporary storage n */
-	  n = mpn_gcd_subdiv_step (up, vp, n, &gcd_hook, &ctx, tp);
+	  n = mpn_gcd_subdiv_step (up, vp, n, gcd_hook, &ctx, tp);
 	  if (n == 0)
 	    goto done;
 	}
diff --git a/mpn/generic/gcd_subdiv_step.c b/mpn/generic/gcd_subdiv_step.c
index 11ccae356..0a5c65bf6 100644
--- a/mpn/generic/gcd_subdiv_step.c
+++ b/mpn/generic/gcd_subdiv_step.c
@@ -21,6 +21,8 @@ License for more details.
 You should have received a copy of the GNU Lesser General Public License
 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
 
+#include <stdlib.h>		/* for NULL */
+
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "longlong.h"
@@ -30,29 +32,39 @@ along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
    followed by one division. Returns zero if the gcd is found.
    Otherwise, compute the reduced a and b, and return the new size. */
 
-/* About the hook functions.
+/* The hook function is called as
+
+     hook(ctx, gp, gn, qp, qn, d)
+
+   in the following cases:
 
-     hook->update (ctx, qp, qn, d)
+   + If A == B at the start, G is the gcd, Q is NULL, d = -1.
 
-   is called when one of the numbers, or a multiple of it, is
-   subtracted from the other. d == 0 means q a is subtracted from b, d
-   == 1 means that q b is subtracted from a.
+   + If one input is zero at the start, G is the gcd, Q is NULL,
+     d = 0 if A == G and d = 1 if B == G.
 
-     hook->done (ctx, gp, gn, d)
+   Otherwise, if d = 0 we have just subtracted a multiple of A from B,
+   and if d = 1 we have subtracted a multiple of B from A.
+   
+   + If A == B after subtraction, G is the gcd, Q is NULL.
+   
+   + If we get a zero remainder after division, G is the gcd, Q is the
+     quotient.
 
-   is called when the gcd is found. d == 0 if the last reduction
-   subtracted a from b, d == 1 if it subtracted b from a, and d == 2
-   if it is unknown which was the most recent reduction. */
+   + Otherwise, G is NULL, Q is the quotient (often 1).
+
+ */
 
 mp_size_t
 mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n,
-		     const struct gcd_subdiv_step_hook *hook, void *ctx,
+		     gcd_subdiv_step_hook *hook, void *ctx,
 		     mp_ptr tp)
 {
-  mp_size_t an, bn;
+  static const mp_limb_t one = CNST_LIMB(1);
+  mp_size_t an, bn, qn;
 
   int swapped;
-  
+
   ASSERT (n > 0);
   ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
 
@@ -71,7 +83,7 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n,
       if (UNLIKELY (c == 0))
 	{
 	  /* For gcdext, return the smallest of the two cofactors. */
-	  hook->done (ctx, ap, an, 2);
+	  hook (ctx, ap, an, NULL, 0, -1);
 	  return 0;
 	}
       else if (c > 0)
@@ -89,7 +101,7 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n,
 	}
       if (an == 0)
 	{
-	  hook->done (ctx, bp, bn, swapped ^ 1);
+	  hook (ctx, bp, bn, NULL, 0, swapped ^ 1);
 	  return 0;
 	}
     }
@@ -97,47 +109,46 @@ mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n,
   MPN_NORMALIZE (bp, bn);
   ASSERT (bn > 0);
 
-  if (hook->update)
-    {
-      static const mp_limb_t one = CNST_LIMB(1);
-      hook->update(ctx, &one, 1, swapped);
-    }
-
-  /* Arrange so that a < b, and divide b = q a + r */
+  /* Arrange so that a < b */
   if (an == bn)
     {
       int c;
       MPN_CMP (c, ap, bp, an);
       if (UNLIKELY (c == 0))
 	{
-	  hook->done (ctx, bp, bn, swapped);
+	  hook (ctx, bp, bn, NULL, 0, swapped);
 	  return 0;	  
 	}
-      else if (c > 0)
+      
+      hook (ctx, NULL, 0, &one, 1, swapped);
+
+      if (c > 0)
 	{
 	  MP_PTR_SWAP (ap, bp);
 	  swapped ^= 1;
 	}
     }
-  else if (an > bn)
+  else
     {
-      MPN_PTR_SWAP (ap, an, bp, bn);
-      swapped ^= 1;
+      hook (ctx, NULL, 0, &one, 1, swapped);
+
+      if (an > bn)
+	{
+	  MPN_PTR_SWAP (ap, an, bp, bn);
+	  swapped ^= 1;
+	}
     }
 
   mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an);
-
-  /* FIXME: For jacobi, we need to call update before calling done.
-     While for gcdext, calling update in this case would do useless
-     work. */
+  qn = bn - an + 1;
   if (mpn_zero_p (bp, an))
     {
-      hook->done (ctx, ap, an, swapped);
+      hook (ctx, ap, an, tp, qn, swapped);
       return 0;
     }
-
-  if (hook->update)
-    hook->update(ctx, tp, bn - an + 1, swapped);
-  
-  return an;
+  else
+    {
+      hook (ctx, NULL, 0, tp, qn, swapped);
+      return an;
+    }
 }
diff --git a/mpn/generic/gcdext.c b/mpn/generic/gcdext.c
index 2cc268113..d314ab94c 100644
--- a/mpn/generic/gcdext.c
+++ b/mpn/generic/gcdext.c
@@ -318,7 +318,7 @@ mpn_gcdext (mp_ptr gp, mp_ptr up, mp_size_t *usizep,
 	ctx.un = 1;
 
 	/* Temporary storage n */
-	n = mpn_gcd_subdiv_step (ap, bp, n, &gcdext_hook, &ctx, tp);
+	n = mpn_gcd_subdiv_step (ap, bp, n, mpn_gcdext_hook, &ctx, tp);
 	if (n == 0)
 	  {
 	    TMP_FREE;
@@ -373,7 +373,7 @@ mpn_gcdext (mp_ptr gp, mp_ptr up, mp_size_t *usizep,
  	  ctx.un = un;
 
 	  /* Temporary storage n */
-	  n = mpn_gcd_subdiv_step (ap, bp, n, &gcdext_hook, &ctx, tp);
+	  n = mpn_gcd_subdiv_step (ap, bp, n, mpn_gcdext_hook, &ctx, tp);
 	  if (n == 0)
 	    {
 	      TMP_FREE;
diff --git a/mpn/generic/gcdext_lehmer.c b/mpn/generic/gcdext_lehmer.c
index 837cd7674..b7288c9c4 100644
--- a/mpn/generic/gcdext_lehmer.c
+++ b/mpn/generic/gcdext_lehmer.c
@@ -24,102 +24,95 @@ along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
 
 /* Here, d is the index of the cofactor to update. FIXME: Could use qn
    = 0 for the common case q = 1. */
-static void
-gcdext_update (void *p, mp_srcptr qp, mp_size_t qn, unsigned d)
+void
+mpn_gcdext_hook (void *p, mp_srcptr gp, mp_size_t gn,
+		 mp_srcptr qp, mp_size_t qn, int d)
 {
   struct gcdext_ctx *ctx = (struct gcdext_ctx *) p;
-  mp_ptr u0;
-  mp_ptr u1;
   mp_size_t un = ctx->un;
-  mp_limb_t cy;
 
-  u0 = ctx->u0;  
-  u1 = ctx->u1;
-  if (d)
-    MP_PTR_SWAP (u0, u1);
+  if (gp)
+    {
+      mp_srcptr up;
 
-  qn -= (qp[qn-1] == 0);
+      MPN_COPY (ctx->gp, gp, gn);
+      ctx->gn = gn;
+      
+      if (d < 0)
+	{
+	  int c;
+      
+	  /* Must return the smallest cofactor, +u1 or -u0 */
+	  MPN_CMP (c, ctx->u0, ctx->u1, un);
+	  ASSERT (c != 0 || (un == 1 && ctx->u0[0] == 1 && ctx->u1[0] == 1));
 
-  /* Update u0 += q  * u1 */
-  if (qn == 1)
-    {
-      mp_limb_t q = qp[0];
+	  d = c < 0;
+	}
 
-      if (q == 1)
-	/* A common case. */
-	cy = mpn_add_n (u0, u0, u1, un);
-      else
-	cy = mpn_addmul_1 (u0, u1, un, q);
+      up = d ? ctx->u0 : ctx->u1;
+
+      MPN_NORMALIZE (up, un);
+      MPN_COPY (ctx->up, up, un);
+
+      *ctx->usize = d ? -un : un;
     }
   else
     {
-      mp_size_t u1n;
-      mp_ptr tp;
-
-      u1n = un;
-      MPN_NORMALIZE (u1, u1n);
-
-      if (u1n == 0)
-	return;
+      mp_limb_t cy;
+      mp_ptr u0 = ctx->u0;
+      mp_ptr u1 = ctx->u1;
 
-      tp = ctx->tp;
+      ASSERT (d >= 0);
+  
+      if (d)
+	MP_PTR_SWAP (u0, u1);
       
-      if (qn > u1n)
-	mpn_mul (tp, qp, qn, u1, u1n);
-      else
-	mpn_mul (tp, u1, u1n, qp, qn);
+      qn -= (qp[qn-1] == 0);
 
-      u1n += qn;
-      u1n -= tp[u1n-1] == 0;
+      /* Update u0 += q  * u1 */
+      if (qn == 1)
+	{
+	  mp_limb_t q = qp[0];
 
-      if (u1n > un)
-	cy = mpn_add (u0, tp, u1n, u0, un);
+	  if (q == 1)
+	    /* A common case. */
+	    cy = mpn_add_n (u0, u0, u1, un);
+	  else
+	    cy = mpn_addmul_1 (u0, u1, un, q);
+	}
       else
-	cy = mpn_add (u0, u0, un, tp, u1n);
-      
-      un = u1n;
-    }
-  u0[un] = cy;
-  ctx->un = un + (cy > 0);
-}
-
-static void
-gcdext_done (void *p, mp_srcptr gp, mp_size_t gn, unsigned d)
-{
-  struct gcdext_ctx *ctx = (struct gcdext_ctx *) p;
-  mp_size_t un;
-  mp_srcptr up;
+	{
+	  mp_size_t u1n;
+	  mp_ptr tp;
 
-  MPN_COPY (ctx->gp, gp, gn);
-  ctx->gn = gn;
+	  u1n = un;
+	  MPN_NORMALIZE (u1, u1n);
 
-  un = ctx->un;
+	  if (u1n == 0)
+	    return;
 
-  if (d == 2)
-    {
-      int c;
+	  tp = ctx->tp;
       
-      /* Must return the smallest cofactor, +u1 or -u0 */
-      MPN_CMP (c, ctx->u0, ctx->u1, un);
-      ASSERT (c != 0 || (un == 1 && ctx->u0[0] == 1 && ctx->u1[0] == 1));
-
-      d = c < 0;
+	  if (qn > u1n)
+	    mpn_mul (tp, qp, qn, u1, u1n);
+	  else
+	    mpn_mul (tp, u1, u1n, qp, qn);
+
+	  u1n += qn;
+	  u1n -= tp[u1n-1] == 0;
+
+	  if (u1n > un)
+	    cy = mpn_add (u0, tp, u1n, u0, un);
+	  else
+	    cy = mpn_add (u0, u0, un, tp, u1n);
+      
+	  un = u1n;
+	}
+      u0[un] = cy;
+      ctx->un = un + (cy > 0);
     }
-
-  up = d ? ctx->u0 : ctx->u1;
-
-  MPN_NORMALIZE (up, un);
-  MPN_COPY (ctx->up, up, un);
-
-  *ctx->usize = d ? -un : un;
 }
 
-const struct gcd_subdiv_step_hook
-gcdext_hook = {
-  gcdext_update,
-  gcdext_done
-};
-
 /* Temporary storage: 3*(n+1) for u. n+1 for the matrix-vector
    multiplications (if hgcd2 succeeds). If hgcd fails, n+1 limbs are
    needed for the division, with most n for the quotient, and n+1 for
@@ -223,7 +216,7 @@ mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize,
 
 	  /* Temporary storage n for the quotient and ualloc for the
 	     new cofactor. */
-	  n = mpn_gcd_subdiv_step (ap, bp, n, &gcdext_hook, &ctx, tp);
+	  n = mpn_gcd_subdiv_step (ap, bp, n, mpn_gcdext_hook, &ctx, tp);
 	  if (n == 0)
 	    return ctx.gn;