[Gmp-commit] /home/hgfiles/gmp: Deleted function mpn_hgcd_lehmer.

mercurial at gmplib.org mercurial at gmplib.org
Sat Feb 20 20:44:31 CET 2010


details:   /home/hgfiles/gmp/rev/d0223fefc57e
changeset: 13433:d0223fefc57e
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Sat Feb 20 20:44:12 2010 +0100
description:
Deleted function mpn_hgcd_lehmer.

diffstat:

 ChangeLog          |   17 ++++++
 gmp-impl.h         |    5 -
 mpn/generic/hgcd.c |  139 +++++++++++++++++++++-------------------------------
 tune/Makefile.am   |    2 +-
 tune/gcd_accel.c   |   28 ----------
 tune/hgcd_lehmer.c |   33 ++++++++++++
 tune/speed.h       |    3 +
 7 files changed, 110 insertions(+), 117 deletions(-)

diffs (truncated from 324 to 300 lines):

diff -r 5bb4e6bea5be -r d0223fefc57e ChangeLog
--- a/ChangeLog	Fri Feb 19 15:29:32 2010 +0100
+++ b/ChangeLog	Sat Feb 20 20:44:12 2010 +0100
@@ -1,3 +1,20 @@
+2010-02-20  Niels Möller  <nisse at lysator.liu.se>
+
+	* tune/speed.h (mpn_gcd_accel): Deleted prototype.
+	(mpn_hgcd_lehmer): New prototype.
+	(MPN_HGCD_LEHMER_ITCH): New macro (previously in gmp-impl.h).
+
+	* tune/Makefile.am (libspeed_la_SOURCES): Added hgcd_lehmer.c.
+	* tune/hgcd_lehmer.c: New file.
+	* tune/gcd_accel.c: Deleted obsolete file.
+
+	* gmp-impl.h (MPN_HGCD_LEHMER_ITCH): Deleted macro.
+
+	* mpn/generic/hgcd.c (mpn_hgcd_lehmer): Deleted function,
+	(mpn_hgcd): Don't call mpn_hgcd_lehmer, instead use inlined loop
+	around hgcd_step.
+	(mpn_hgcd_itch): Substitute n for MPN_HGCD_LEHMER_ITCH (n).
+
 2010-02-19  Niels Möller  <nisse at lysator.liu.se>
 
 	* Makefile.am (mpn/jacobitab.h): Added the rules needed to
diff -r 5bb4e6bea5be -r d0223fefc57e gmp-impl.h
--- a/gmp-impl.h	Fri Feb 19 15:29:32 2010 +0100
+++ b/gmp-impl.h	Sat Feb 20 20:44:12 2010 +0100
@@ -3791,11 +3791,6 @@
 #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));
 
-#define MPN_HGCD_LEHMER_ITCH(n) (n)
-
-#define mpn_hgcd_lehmer __MPN (hgcd_lehmer)
-__GMP_DECLSPEC mp_size_t mpn_hgcd_lehmer __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, struct hgcd_matrix *, mp_ptr));
-
 /* Needs storage for the quotient */
 #define MPN_GCD_SUBDIV_STEP_ITCH(n) (n)
 
diff -r 5bb4e6bea5be -r d0223fefc57e mpn/generic/hgcd.c
--- a/mpn/generic/hgcd.c	Fri Feb 19 15:29:32 2010 +0100
+++ b/mpn/generic/hgcd.c	Sat Feb 20 20:44:12 2010 +0100
@@ -356,32 +356,6 @@
   return bn;
 }
 
-/* Reduces a,b until |a-b| fits in n/2 + 1 limbs. Constructs matrix M
-   with elements of size at most (n+1)/2 - 1. Returns new size of a,
-   b, or zero if no reduction is possible. */
-mp_size_t
-mpn_hgcd_lehmer (mp_ptr ap, mp_ptr bp, mp_size_t n,
-		 struct hgcd_matrix *M, mp_ptr tp)
-{
-  mp_size_t s = n/2 + 1;
-  mp_size_t nn;
-
-  ASSERT (n > s);
-  ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
-
-  nn = hgcd_step (n, ap, bp, s, M, tp);
-  if (!nn)
-    return 0;
-
-  for (;;)
-    {
-      n = nn;
-      ASSERT (n > s);
-      nn = hgcd_step (n, ap, bp, s, M, tp);
-      if (!nn )
-	return n;
-    }
-}
 
 /* Multiply M by M1 from the right. Needs 3*(M->n + M1->n) + 5 limbs
    of temporary storage (see mpn_matrix22_mul_itch). */
@@ -530,15 +504,14 @@
   mp_size_t nscaled;
 
   if (BELOW_THRESHOLD (n, HGCD_THRESHOLD))
-    return MPN_HGCD_LEHMER_ITCH (n);
+    return n;
 
   /* Get the recursion depth. */
   nscaled = (n - 1) / (HGCD_THRESHOLD - 1);
   count_leading_zeros (count, nscaled);
   k = GMP_LIMB_BITS - count;
 
-  return 20 * ((n+3) / 4) + 22 * k
-    + MPN_HGCD_LEHMER_ITCH (HGCD_THRESHOLD);
+  return 20 * ((n+3) / 4) + 22 * k + HGCD_THRESHOLD;
 }
 
 /* Reduces a,b until |a-b| fits in n/2 + 1 limbs. Constructs matrix M
@@ -550,9 +523,8 @@
 	  struct hgcd_matrix *M, mp_ptr tp)
 {
   mp_size_t s = n/2 + 1;
-  mp_size_t n2 = (3*n)/4 + 1;
 
-  mp_size_t p, nn;
+  mp_size_t nn;
   int success = 0;
 
   if (n <= s)
@@ -564,72 +536,73 @@
 
   ASSERT ((n+1)/2 - 1 < M->alloc);
 
-  if (BELOW_THRESHOLD (n, HGCD_THRESHOLD))
-    return mpn_hgcd_lehmer (ap, bp, n, M, tp);
+  if (ABOVE_THRESHOLD (n, HGCD_THRESHOLD))
+    {
+      mp_size_t n2 = (3*n)/4 + 1;
+      mp_size_t p = n/2;
 
-  p = n/2;
-  nn = mpn_hgcd (ap + p, bp + p, n - p, M, tp);
-  if (nn > 0)
-    {
-      /* Needs 2*(p + M->n) <= 2*(floor(n/2) + ceil(n/2) - 1)
-	 = 2 (n - 1) */
-      n = mpn_hgcd_matrix_adjust (M, p + nn, ap, bp, p, tp);
-      success = 1;
-    }
-  while (n > n2)
-    {
-      /* Needs n + 1 storage */
-      nn = hgcd_step (n, ap, bp, s, M, tp);
-      if (!nn)
-	return success ? n : 0;
-      n = nn;
-      success = 1;
-    }
-
-  if (n > s + 2)
-    {
-      struct hgcd_matrix M1;
-      mp_size_t scratch;
-
-      p = 2*s - n + 1;
-      scratch = MPN_HGCD_MATRIX_INIT_ITCH (n-p);
-
-      mpn_hgcd_matrix_init(&M1, n - p, tp);
-      nn = mpn_hgcd (ap + p, bp + p, n - p, &M1, tp + scratch);
+      nn = mpn_hgcd (ap + p, bp + p, n - p, M, tp);
       if (nn > 0)
 	{
-	  /* We always have max(M) > 2^{-(GMP_NUMB_BITS + 1)} max(M1) */
-	  ASSERT (M->n + 2 >= M1.n);
+	  /* Needs 2*(p + M->n) <= 2*(floor(n/2) + ceil(n/2) - 1)
+	     = 2 (n - 1) */
+	  n = mpn_hgcd_matrix_adjust (M, p + nn, ap, bp, p, tp);
+	  success = 1;
+	}
+      while (n > n2)
+	{
+	  /* Needs n + 1 storage */
+	  nn = hgcd_step (n, ap, bp, s, M, tp);
+	  if (!nn)
+	    return success ? n : 0;
+	  n = nn;
+	  success = 1;
+	}
 
-	  /* Furthermore, assume M ends with a quotient (1, q; 0, 1),
-	     then either q or q + 1 is a correct quotient, and M1 will
-	     start with either (1, 0; 1, 1) or (2, 1; 1, 1). This
-	     rules out the case that the size of M * M1 is much
-	     smaller than the expected M->n + M1->n. */
+      if (n > s + 2)
+	{
+	  struct hgcd_matrix M1;
+	  mp_size_t scratch;
 
-	  ASSERT (M->n + M1.n < M->alloc);
+	  p = 2*s - n + 1;
+	  scratch = MPN_HGCD_MATRIX_INIT_ITCH (n-p);
 
-	  /* Needs 2 (p + M->n) <= 2 (2*s - n2 + 1 + n2 - s - 1)
-	     = 2*s <= 2*(floor(n/2) + 1) <= n + 2. */
-	  n = mpn_hgcd_matrix_adjust (&M1, p + nn, ap, bp, p, tp + scratch);
+	  mpn_hgcd_matrix_init(&M1, n - p, tp);
+	  nn = mpn_hgcd (ap + p, bp + p, n - p, &M1, tp + scratch);
+	  if (nn > 0)
+	    {
+	      /* We always have max(M) > 2^{-(GMP_NUMB_BITS + 1)} max(M1) */
+	      ASSERT (M->n + 2 >= M1.n);
 
-	  /* We need a bound for of M->n + M1.n. Let n be the original
-	     input size. Then
+	      /* Furthermore, assume M ends with a quotient (1, q; 0, 1),
+		 then either q or q + 1 is a correct quotient, and M1 will
+		 start with either (1, 0; 1, 1) or (2, 1; 1, 1). This
+		 rules out the case that the size of M * M1 is much
+		 smaller than the expected M->n + M1->n. */
 
-	       ceil(n/2) - 1 >= size of product >= M.n + M1.n - 2
+	      ASSERT (M->n + M1.n < M->alloc);
 
-	     and it follows that
+	      /* Needs 2 (p + M->n) <= 2 (2*s - n2 + 1 + n2 - s - 1)
+		 = 2*s <= 2*(floor(n/2) + 1) <= n + 2. */
+	      n = mpn_hgcd_matrix_adjust (&M1, p + nn, ap, bp, p, tp + scratch);
 
-	       M.n + M1.n <= ceil(n/2) + 1
+	      /* We need a bound for of M->n + M1.n. Let n be the original
+		 input size. Then
 
-	     Then 3*(M.n + M1.n) + 5 <= 3 * ceil(n/2) + 8 is the
-	     amount of needed scratch space. */
-	  mpn_hgcd_matrix_mul (M, &M1, tp + scratch);
-	  success = 1;
+		 ceil(n/2) - 1 >= size of product >= M.n + M1.n - 2
+
+		 and it follows that
+
+		 M.n + M1.n <= ceil(n/2) + 1
+
+		 Then 3*(M.n + M1.n) + 5 <= 3 * ceil(n/2) + 8 is the
+		 amount of needed scratch space. */
+	      mpn_hgcd_matrix_mul (M, &M1, tp + scratch);
+	      success = 1;
+	    }
 	}
     }
 
-  /* This really is the base case */
   for (;;)
     {
       /* Needs s+3 < n */
diff -r 5bb4e6bea5be -r d0223fefc57e tune/Makefile.am
--- a/tune/Makefile.am	Fri Feb 19 15:29:32 2010 +0100
+++ b/tune/Makefile.am	Sat Feb 20 20:44:12 2010 +0100
@@ -43,7 +43,7 @@
   common.c divrem1div.c divrem1inv.c divrem2div.c divrem2inv.c		\
   freq.c								\
   gcdext_single.c gcdext_double.c gcdextod.c gcdextos.c			\
-  jacbase1.c jacbase2.c jacbase3.c					\
+  hgcd_lehmer.c jacbase1.c jacbase2.c jacbase3.c			\
   mod_1_div.c mod_1_inv.c modlinv.c					\
   noop.c powm_mod.c powm_redc.c pre_divrem_1.c				\
   set_strb.c set_strs.c set_strp.c time.c
diff -r 5bb4e6bea5be -r d0223fefc57e tune/gcd_accel.c
--- a/tune/gcd_accel.c	Fri Feb 19 15:29:32 2010 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,28 +0,0 @@
-/* mpn/generic/gcd.c forced to use the accelerated binary algorithm. */
-
-/*
-Copyright 2003 Free Software Foundation, Inc.
-
-This file is part of the GNU MP Library.
-
-The GNU MP Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-The GNU MP Library is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
-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 "gmp.h"
-#include "gmp-impl.h"
-
-#undef  GCD_SCHOENHAGE_THRESHOLD
-#define GCD_SCHOENHAGE_THRESHOLD  MP_SIZE_T_MAX
-#define __gmpn_gcd  mpn_gcd_accel
-
-#include "../mpn/generic/gcd.c"
diff -r 5bb4e6bea5be -r d0223fefc57e tune/hgcd_lehmer.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tune/hgcd_lehmer.c	Sat Feb 20 20:44:12 2010 +0100
@@ -0,0 +1,33 @@
+/* mpn/generic/hgcd.c forced to use Lehmer's quadratic algorithm. */
+
+/*
+Copyright 2010 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+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 "gmp.h"
+#include "gmp-impl.h"


More information about the gmp-commit mailing list