[Gmp-commit] /home/hgfiles/gmp: New variant for mpn_jacobi_base.
mercurial at gmplib.org
mercurial at gmplib.org
Wed Mar 10 17:00:08 CET 2010
details: /home/hgfiles/gmp/rev/9ddaaacc501f
changeset: 13476:9ddaaacc501f
user: Niels M?ller <nisse at lysator.liu.se>
date: Wed Mar 10 16:47:01 2010 +0100
description:
New variant for mpn_jacobi_base.
diffstat:
ChangeLog | 12 ++++++++++
mpn/generic/jacbase.c | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++-
tune/Makefile.am | 2 +-
tune/common.c | 5 ++++
tune/jacbase4.c | 27 +++++++++++++++++++++++
tune/speed.c | 1 +
tune/speed.h | 1 +
tune/tuneup.c | 16 +++++++++----
8 files changed, 115 insertions(+), 7 deletions(-)
diffs (213 lines):
diff -r cad9b7b03dcd -r 9ddaaacc501f ChangeLog
--- a/ChangeLog Wed Mar 10 07:58:51 2010 +0100
+++ b/ChangeLog Wed Mar 10 16:47:01 2010 +0100
@@ -1,3 +1,15 @@
+2010-03-10 Niels Möller <nisse at lysator.liu.se>
+
+ * tune/tuneup.c (tune_jacobi_base): Consider mpn_jacobi_base_4.
+ * tune/speed.c (routine): Added mpn_jacobi_base_4.
+ * tune/common.c (speed_mpn_jacobi_base_4): New function.
+ * tune/speed.h (speed_mpn_jacobi_base_4): Declare it.
+ * tune/Makefile.am (libspeed_la_SOURCES): Added jacbase4.c.
+ * tune/jacbase4.c: New file.
+
+ * mpn/generic/jacbase.c (mpn_jacobi_base): New function, for
+ JACOBI_BASE_METHOD 4.
+
2010-03-09 Niels Möller <nisse at lysator.liu.se>
* tests/mpz/t-jac.c (check_large_quotients): Also generate inputs
diff -r cad9b7b03dcd -r 9ddaaacc501f mpn/generic/jacbase.c
--- a/mpn/generic/jacbase.c Wed Mar 10 07:58:51 2010 +0100
+++ b/mpn/generic/jacbase.c Wed Mar 10 16:47:01 2010 +0100
@@ -108,7 +108,7 @@
}
#endif
-
+#if JACOBI_BASE_METHOD < 4
/* Calculate the value of the Jacobi symbol (a/b) of two mp_limb_t's, but
with a restricted range of inputs accepted, namely b>1, b odd, and a<=b.
@@ -166,3 +166,59 @@
done:
return JACOBI_BIT1_TO_PN (result_bit1);
}
+#endif
+
+#if JACOBI_BASE_METHOD == 4
+/* Computes (a/b) for odd b and any a. The initial bit is taken as a
+ * parameter. We have no need for the convention that the sign is in
+ * bit 1, internally we use bit 0. */
+int
+mpn_jacobi_base (mp_limb_t a, mp_limb_t b, int bit)
+{
+ int c;
+
+ ASSERT (b & 1);
+
+ if (a == 0)
+ /* Ok, here we use that the sign is bit 1, after all. */
+ return b == 1 ? (1-(bit & 2)) : 0;
+
+ bit >>= 1;
+
+ /* Below, we represent a and b shifted right so that the least
+ significant one bit is implicit. */
+
+ b >>= 1;
+
+ count_trailing_zeros (c, a);
+ bit ^= c & (b ^ (b >> 1));
+
+ /* We may have c==GMP_LIMB_BITS-1, so we cant use a>>c+1. */
+ a >>= c;
+ a >>= 1;
+
+ while (a != b)
+ {
+ mp_limb_t t = a - b;
+ mp_limb_t bgta = LIMB_HIGHBIT_TO_MASK (t);
+
+ /* If b > a, invoke reciprocity */
+ bit ^= (bgta & a & b);
+
+ /* b <-- min (a, b) */
+ b += (bgta & t);
+
+ /* a <-- |a - b| */
+ a = (t ^ bgta) - bgta;
+
+ /* Number of trailing zeros is the same no matter if we look at
+ * t or a, but using t gives more parallelism. */
+ count_trailing_zeros (c, t);
+ c ++;
+ /* (2/b) = -1 if b = 3 or 5 mod 8 */
+ bit ^= c & (b ^ (b >> 1));
+ a >>= c;
+ }
+ return a == 0 ? 1-2*(bit & 1) : 0;
+}
+#endif /* JACOBI_BASE_METHOD == 4 */
diff -r cad9b7b03dcd -r 9ddaaacc501f tune/Makefile.am
--- a/tune/Makefile.am Wed Mar 10 07:58:51 2010 +0100
+++ b/tune/Makefile.am Wed Mar 10 16:47:01 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 \
- hgcd_lehmer.c jacbase1.c jacbase2.c jacbase3.c \
+ hgcd_lehmer.c jacbase1.c jacbase2.c jacbase3.c jacbase4.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 cad9b7b03dcd -r 9ddaaacc501f tune/common.c
--- a/tune/common.c Wed Mar 10 07:58:51 2010 +0100
+++ b/tune/common.c Wed Mar 10 16:47:01 2010 +0100
@@ -1491,6 +1491,11 @@
{
SPEED_ROUTINE_MPN_JACBASE (mpn_jacobi_base_3);
}
+double
+speed_mpn_jacobi_base_4 (struct speed_params *s)
+{
+ SPEED_ROUTINE_MPN_JACBASE (mpn_jacobi_base_4);
+}
double
diff -r cad9b7b03dcd -r 9ddaaacc501f tune/jacbase4.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tune/jacbase4.c Wed Mar 10 16:47:01 2010 +0100
@@ -0,0 +1,27 @@
+/* mpn/generic/jacbase.c method 4.
+
+Copyright 2002, 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"
+
+#undef JACOBI_BASE_METHOD
+#define JACOBI_BASE_METHOD 4
+#define __gmpn_jacobi_base mpn_jacobi_base_4
+
+#include "mpn/generic/jacbase.c"
diff -r cad9b7b03dcd -r 9ddaaacc501f tune/speed.c
--- a/tune/speed.c Wed Mar 10 07:58:51 2010 +0100
+++ b/tune/speed.c Wed Mar 10 16:47:01 2010 +0100
@@ -283,6 +283,7 @@
{ "mpn_jacobi_base_1", speed_mpn_jacobi_base_1 },
{ "mpn_jacobi_base_2", speed_mpn_jacobi_base_2 },
{ "mpn_jacobi_base_3", speed_mpn_jacobi_base_3 },
+ { "mpn_jacobi_base_4", speed_mpn_jacobi_base_4 },
{ "mpn_mul", speed_mpn_mul, FLAG_R_OPTIONAL },
{ "mpn_mul_basecase", speed_mpn_mul_basecase,FLAG_R_OPTIONAL },
diff -r cad9b7b03dcd -r 9ddaaacc501f tune/speed.h
--- a/tune/speed.h Wed Mar 10 07:58:51 2010 +0100
+++ b/tune/speed.h Wed Mar 10 16:47:01 2010 +0100
@@ -200,6 +200,7 @@
double speed_mpn_jacobi_base_1 __GMP_PROTO ((struct speed_params *s));
double speed_mpn_jacobi_base_2 __GMP_PROTO ((struct speed_params *s));
double speed_mpn_jacobi_base_3 __GMP_PROTO ((struct speed_params *s));
+double speed_mpn_jacobi_base_4 __GMP_PROTO ((struct speed_params *s));
double speed_mpn_lshift __GMP_PROTO ((struct speed_params *s));
double speed_mpn_lshiftc __GMP_PROTO ((struct speed_params *s));
double speed_mpn_mod_1 __GMP_PROTO ((struct speed_params *s));
diff -r cad9b7b03dcd -r 9ddaaacc501f tune/tuneup.c
--- a/tune/tuneup.c Wed Mar 10 07:58:51 2010 +0100
+++ b/tune/tuneup.c Wed Mar 10 16:47:01 2010 +0100
@@ -2091,7 +2091,7 @@
tune_jacobi_base (void)
{
static struct param_t param;
- double t1, t2, t3;
+ double t1, t2, t3, t4;
int method;
s.size = GMP_LIMB_BITS * 3 / 4;
@@ -2108,19 +2108,25 @@
if (option_trace >= 1)
printf ("size=%ld, mpn_jacobi_base_3 %.9f\n", (long) s.size, t3);
- if (t1 == -1.0 || t2 == -1.0 || t3 == -1.0)
+ t4 = tuneup_measure (speed_mpn_jacobi_base_4, ¶m, &s);
+ if (option_trace >= 1)
+ printf ("size=%ld, mpn_jacobi_base_4 %.9f\n", (long) s.size, t4);
+
+ if (t1 == -1.0 || t2 == -1.0 || t3 == -1.0 || t4 == -1.0)
{
printf ("Oops, can't measure all mpn_jacobi_base methods at %ld\n",
(long) s.size);
abort ();
}
- if (t1 < t2 && t1 < t3)
+ if (t1 < t2 && t1 < t3 && t1 < t4)
method = 1;
- else if (t2 < t3)
+ else if (t2 < t3 && t2 < t4)
method = 2;
+ else if (t3 < t4)
+ method = 3;
else
- method = 3;
+ method = 4;
print_define ("JACOBI_BASE_METHOD", method);
}
More information about the gmp-commit
mailing list