[Gmp-commit] /var/hg/gmp: 2 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Fri May 20 20:29:20 CEST 2011
details: /var/hg/gmp/rev/4d399272d8d9
changeset: 14194:4d399272d8d9
user: Niels M?ller <nisse at lysator.liu.se>
date: Fri May 20 20:20:30 2011 +0200
description:
Subquadratic jacobi.
details: /var/hg/gmp/rev/ddd56c4f2b69
changeset: 14195:ddd56c4f2b69
user: Niels M?ller <nisse at lysator.liu.se>
date: Fri May 20 20:23:55 2011 +0200
description:
Ignore gen-jacobitab and mpn/jacobitab.h.
diffstat:
.hgignore | 2 +
ChangeLog | 14 +
configure.in | 5 +-
gmp-impl.h | 13 +-
mpn/generic/jacobi.c | 284 ++++++++++++++
mpn/generic/jacobi_lehmer.c | 859 --------------------------------------------
mpz/jacobi.c | 324 ++++++++++------
7 files changed, 517 insertions(+), 984 deletions(-)
diffs (truncated from 1635 to 300 lines):
diff -r 437e5faee991 -r ddd56c4f2b69 .hgignore
--- a/.hgignore Fri May 20 19:22:15 2011 +0200
+++ b/.hgignore Fri May 20 20:23:55 2011 +0200
@@ -23,6 +23,7 @@
^gen-bases
^gen-fac_ui
^gen-fib
+^gen-jacobitab
^gen-psqr
^gen-trialdivtab
^mp_bases\.h
@@ -30,6 +31,7 @@
^trialdivtab\.h
^mpn/perfsqr\.h
+^mpn/jacobitab\.h
^mpz/fac_ui\.h
^doc/gmp\.info.*
diff -r 437e5faee991 -r ddd56c4f2b69 ChangeLog
--- a/ChangeLog Fri May 20 19:22:15 2011 +0200
+++ b/ChangeLog Fri May 20 20:23:55 2011 +0200
@@ -1,5 +1,19 @@
2011-05-20 Niels Möller <nisse at lysator.liu.se>
+ * gmp-impl.h: Jacobi-related prototypes.
+
+ * configure.in (gmp_mpn_functions): Added jacobi_2, jacobi,
+ hgcd2_jacobi, hgcd_jacobi, and removed jacobi_lehmer.
+
+ * mpz/jacobi.c (STRIP_TWOS): Deleted macro.
+ (mpz_jacobi): Partially rewritten, to no longer makes the A
+ operand odd. Use new mpn_jacobi_n.
+
+ * mpn/generic/jacobi_lehmer.c: Deleted file.
+
+ * mpn/generic/jacobi.c (mpn_jacobi_n): New subquadratic jacobi
+ implementation. Supersedes jacobi_lehmer.c.
+
* mpn/generic/hgcd_jacobi.c (mpn_hgcd_jacobi): New file and
function. A copy of mpn_hgcd, using mpn_hgcd2_jacobi, and with calls to
mpn_jacobi_update when appropriate.
diff -r 437e5faee991 -r ddd56c4f2b69 configure.in
--- a/configure.in Fri May 20 19:22:15 2011 +0200
+++ b/configure.in Fri May 20 20:23:55 2011 +0200
@@ -2538,9 +2538,10 @@
perfsqr perfpow \
gcd_1 gcd gcdext_1 gcdext gcd_subdiv_step \
gcdext_lehmer \
- div_q tdiv_qr jacbase jacobi_lehmer get_d \
+ div_q tdiv_qr jacbase jacobi_2 jacobi get_d \
matrix22_mul matrix22_mul1_inverse_vector \
- hgcd_matrix hgcd2 hgcd mullo_n mullo_basecase \
+ hgcd_matrix hgcd2 hgcd hgcd2_jacobi hgcd_jacobi \
+ mullo_n mullo_basecase \
toom22_mul toom32_mul toom42_mul toom52_mul toom62_mul \
toom33_mul toom43_mul toom53_mul toom63_mul \
toom44_mul \
diff -r 437e5faee991 -r ddd56c4f2b69 gmp-impl.h
--- a/gmp-impl.h Fri May 20 19:22:15 2011 +0200
+++ b/gmp-impl.h Fri May 20 20:23:55 2011 +0200
@@ -985,8 +985,11 @@
#define mpn_jacobi_base __MPN(jacobi_base)
__GMP_DECLSPEC int mpn_jacobi_base __GMP_PROTO ((mp_limb_t, mp_limb_t, int)) ATTRIBUTE_CONST;
-#define mpn_jacobi_lehmer __MPN(jacobi_lehmer)
-__GMP_DECLSPEC int mpn_jacobi_lehmer __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, unsigned, mp_ptr));
+#define mpn_jacobi_2 __MPN(jacobi_2)
+__GMP_DECLSPEC int mpn_jacobi_2 __GMP_PROTO ((mp_srcptr, mp_srcptr, unsigned));
+
+#define mpn_jacobi_n __MPN(jacobi_n)
+__GMP_DECLSPEC int mpn_jacobi_n __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, unsigned));
#define mpn_mod_1c __MPN(mod_1c)
__GMP_DECLSPEC mp_limb_t mpn_mod_1c __GMP_PROTO ((mp_srcptr, mp_size_t, mp_limb_t, mp_limb_t)) __GMP_ATTRIBUTE_PURE;
@@ -3914,6 +3917,9 @@
#define mpn_matrix22_mul1_inverse_vector __MPN (matrix22_mul1_inverse_vector)
__GMP_DECLSPEC mp_size_t mpn_matrix22_mul1_inverse_vector __GMP_PROTO ((const struct hgcd_matrix1 *, mp_ptr, mp_srcptr, mp_ptr, mp_size_t));
+#define mpn_hgcd2_jacobi __MPN (hgcd2_jacobi)
+__GMP_DECLSPEC int mpn_hgcd2_jacobi __GMP_PROTO ((mp_limb_t, mp_limb_t, mp_limb_t, mp_limb_t, struct hgcd_matrix1 *, unsigned *));
+
struct hgcd_matrix
{
mp_size_t alloc; /* for sanity checking only */
@@ -3944,6 +3950,9 @@
#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_jacobi __MPN (hgcd_jacobi)
+__GMP_DECLSPEC mp_size_t mpn_hgcd_jacobi __GMP_PROTO ((mp_ptr, mp_ptr, mp_size_t, struct hgcd_matrix *, unsigned *, mp_ptr));
+
typedef void gcd_subdiv_step_hook(void *, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, int);
/* Needs storage for the quotient */
diff -r 437e5faee991 -r ddd56c4f2b69 mpn/generic/jacobi.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/jacobi.c Fri May 20 20:23:55 2011 +0200
@@ -0,0 +1,284 @@
+/* jacobi.c
+
+ 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
+ GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
+
+Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2008, 2010, 2011 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"
+#include "longlong.h"
+
+#ifndef JACOBI_DC_THRESHOLD
+#define JACOBI_DC_THRESHOLD GCD_DC_THRESHOLD
+#endif
+
+/* Schönhage's rules:
+ *
+ * Assume r0 = r1 q1 + r2, with r0 odd, and r1 = q2 r2 + r3
+ *
+ * If r1 is odd, then
+ *
+ * (r1 | r0) = s(r1, r0) (r0 | r1) = s(r1, r0) (r2, r1)
+ *
+ * where s(x,y) = (-1)^{(x-1)(y-1)/4} = (-1)^[x = y = 3 (mod 4)].
+ *
+ * If r1 is even, r2 must be odd. We have
+ *
+ * (r1 | r0) = (r1 - r0 | r0) = (-1)^(r0-1)/2 (r0 - r1 | r0)
+ * = (-1)^(r0-1)/2 s(r0, r0 - r1) (r0 | r0 - r1)
+ * = (-1)^(r0-1)/2 s(r0, r0 - r1) (r1 | r0 - r1)
+ *
+ * Now, if r1 = 0 (mod 4), then the sign factor is +1, and repeating
+ * q1 times gives
+ *
+ * (r1 | r0) = (r1 | r2) = (r3 | r2)
+ *
+ * On the other hand, if r1 = 2 (mod 4), the sign factor is
+ * (-1)^{(r0-1)/2}, and repeating q1 times gives the exponent
+ *
+ * (r0-1)/2 + (r0-r1-1)/2 + ... + (r0 - (q1-1) r1)/2
+ * = q1 (r0-1)/2 + q1 (q1-1)/2
+ *
+ * and we can summarize the even case as
+ *
+ * (r1 | r0) = t(r1, r0, q1) (r3 | r2)
+ *
+ * where t(x,y,q) = (-1)^{[x = 2 (mod 4)] (q(y-1)/2 + y(q-1)/2)}
+ *
+ * What about termination? The remainder sequence ends with (0|1) = 1
+ * (or (0 | r) = 0 if r != 1). What are the possible cases? If r1 is
+ * odd, r2 may be zero. If r1 is even, then r2 = r0 - q1 r1 is odd and
+ * hence non-zero. We may have r3 = r1 - q2 r2 = 0.
+ *
+ * Examples: (11|15) = - (15|11) = - (4|11)
+ * (4|11) = (4| 3) = (1| 3)
+ * (1| 3) = (3|1) = (0|1) = 1
+ *
+ * (2|7) = (2|1) = (0|1) = 1
+ *
+ * Detail: (2|7) = (2-7|7) = (-1|7)(5|7) = -(7|5) = -(2|5)
+ * (2|5) = (2-5|5) = (-1|5)(3|5) = (5|3) = (2|3)
+ * (2|3) = (2-3|3) = (-1|3)(1|3) = -(3|1) = -(2|1)
+ *
+ */
+
+/* In principle, the state consists of four variables: e (one bit), a,
+ b (two bits each), d (one bit). Collected factors are (-1)^e. a and
+ b are the least significant bits of the current remainders. d
+ (denominator) is 0 if we're currently subtracting multiplies of a
+ from b, and 1 if we're subtracting b from a.
+
+ e is stored in the least significant bit, while a, b and d are
+ coded as only 13 distinct values in bits 1-4, according to the
+ following table. For rows not mentioning d, the value is either
+ implied, or it doesn't matter. */
+
+#if WANT_ASSERT
+static const struct
+{
+ unsigned char a;
+ unsigned char b;
+} decode_table[13] = {
+ /* 0 */ { 0, 1 },
+ /* 1 */ { 0, 3 },
+ /* 2 */ { 1, 1 },
+ /* 3 */ { 1, 3 },
+ /* 4 */ { 2, 1 },
+ /* 5 */ { 2, 3 },
+ /* 6 */ { 3, 1 },
+ /* 7 */ { 3, 3 }, /* d = 1 */
+ /* 8 */ { 1, 0 },
+ /* 9 */ { 1, 2 },
+ /* 10 */ { 3, 0 },
+ /* 11 */ { 3, 2 },
+ /* 12 */ { 3, 3 }, /* d = 0 */
+};
+#define JACOBI_A(bits) (decode_table[(bits)>>1].a)
+#define JACOBI_B(bits) (decode_table[(bits)>>1].b)
+#endif /* WANT_ASSERT */
+
+const unsigned char jacobi_table[208] = {
+#include "jacobitab.h"
+};
+
+#define BITS_FAIL 31
+
+static void
+jacobi_hook (void *p, mp_srcptr gp, mp_size_t gn,
+ mp_srcptr qp, mp_size_t qn, int d)
+{
+ unsigned *bitsp = (unsigned *) p;
+
+ if (gp)
+ {
+ ASSERT (gn > 0);
+ if (gn != 1 || gp[0] != 1)
+ {
+ *bitsp = BITS_FAIL;
+ return;
+ }
+ }
+
+ if (qp)
+ {
+ ASSERT (qn > 0);
+ ASSERT (d >= 0);
+ *bitsp = mpn_jacobi_update (*bitsp, d, qp[0] & 3);
+ }
+}
+
+#define CHOOSE_P(n) (2*(n) / 3)
+
+int
+mpn_jacobi_n (mp_ptr ap, mp_ptr bp, mp_size_t n, unsigned bits)
+{
+ mp_size_t scratch;
+ mp_size_t matrix_scratch;
+ mp_ptr tp;
+
+ TMP_DECL;
+
+ ASSERT (n > 0);
+ ASSERT ( (ap[n-1] | bp[n-1]) > 0);
+ ASSERT ( (bp[0] | ap[0]) & 1);
+
+ /* FIXME: Check for small sizes first, before setting up temporary
+ storage etc. */
+ scratch = MPN_GCD_SUBDIV_STEP_ITCH(n);
+
+ if (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
+ {
+ mp_size_t hgcd_scratch;
+ mp_size_t update_scratch;
+ mp_size_t p = CHOOSE_P (n);
+ mp_size_t dc_scratch;
+
+ matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
+ hgcd_scratch = mpn_hgcd_itch (n - p);
+ update_scratch = p + n - 1;
+
+ dc_scratch = matrix_scratch + MAX(hgcd_scratch, update_scratch);
+ if (dc_scratch > scratch)
+ scratch = dc_scratch;
+ }
+
+ TMP_MARK;
+ tp = TMP_ALLOC_LIMBS(scratch);
+
+ while (ABOVE_THRESHOLD (n, JACOBI_DC_THRESHOLD))
+ {
+ struct hgcd_matrix M;
+ mp_size_t p = 2*n/3;
+ mp_size_t matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
+ mp_size_t nn;
+ mpn_hgcd_matrix_init (&M, n - p, tp);
+
+ nn = mpn_hgcd_jacobi (ap + p, bp + p, n - p, &M, &bits,
+ tp + matrix_scratch);
+ if (nn > 0)
+ {
+ ASSERT (M.n <= (n - p - 1)/2);
+ ASSERT (M.n + p <= (p + n - 1) / 2);
+ /* Temporary storage 2 (p + M->n) <= p + n - 1. */
+ n = mpn_hgcd_matrix_adjust (&M, p + nn, ap, bp, p, tp + matrix_scratch);
More information about the gmp-commit
mailing list