[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