[Gmp-commit] /home/hgfiles/gmp: 2 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Tue May 11 16:22:55 CEST 2010
details: /home/hgfiles/gmp/rev/49abac073223
changeset: 13623:49abac073223
user: Niels M?ller <nisse at lysator.liu.se>
date: Tue May 11 16:08:19 2010 +0200
description:
Minor comment fixes.
details: /home/hgfiles/gmp/rev/5b0baf874eb9
changeset: 13624:5b0baf874eb9
user: Niels M?ller <nisse at lysator.liu.se>
date: Tue May 11 16:09:10 2010 +0200
description:
Enabled new Jacobi code. Reorganized mpz_jacobi to handle small inputs more efficiently.
diffstat:
ChangeLog | 9 +
gmp-impl.h | 8 +-
mpz/jacobi.c | 477 ++++++++++++++++-------------------------------------
tests/mpz/t-jac.c | 36 ++-
4 files changed, 181 insertions(+), 349 deletions(-)
diffs (truncated from 694 to 300 lines):
diff -r 563e3c3cf26a -r 5b0baf874eb9 ChangeLog
--- a/ChangeLog Mon May 10 11:22:02 2010 +0200
+++ b/ChangeLog Tue May 11 16:09:10 2010 +0200
@@ -1,3 +1,12 @@
+2010-05-11 Niels Möller <nisse at lysator.liu.se>
+
+ * mpz/jacobi.c (mpz_jacobi): Deleted old implementation.
+ Reorganized new implementation, to handle small inputs effciently.
+
+ * tests/mpz/t-jac.c (check_large_quotients): Reduced test sizes.
+ (check_data): One more input pair related to a fixed bug.
+ (main): Enable check_large_quotients.
+
2010-05-10 Torbjorn Granlund <tege at gmplib.org>
* mpn/x86_64/aorrlsh2_n.asm: Fix typo.
diff -r 563e3c3cf26a -r 5b0baf874eb9 gmp-impl.h
--- a/gmp-impl.h Mon May 10 11:22:02 2010 +0200
+++ b/gmp-impl.h Tue May 11 16:09:10 2010 +0200
@@ -3572,9 +3572,6 @@
#define PP_FIRST_OMITTED 3
#endif
-
-/* FIXME: Macros for old jacobi code, review what's still needed. */
-
/* BIT1 means a result value in bit 1 (second least significant bit), with a
zero bit representing +1 and a one bit representing -1. Bits other than
bit 1 are garbage. These are meant to be kept in "int"s, and casts are
@@ -3592,6 +3589,9 @@
/* (a/0), with a unsigned; is 1 if a=+/-1, 0 otherwise */
#define JACOBI_U0(a) ((a) == 1)
+/* FIXME: JACOBI_LS0 and JACOBI_0LS are the same, so delete one and
+ come up with a better name. */
+
/* (a/0), with a given by low and size;
is 1 if a=+/-1, 0 otherwise */
#define JACOBI_LS0(alow,asize) \
@@ -3719,7 +3719,7 @@
} \
} while (0)
-/* For new jacobi code. */
+/* State for the Jacobi computation using Lehmer. */
#define jacobi_table __gmp_jacobi_table
__GMP_DECLSPEC extern const unsigned char jacobi_table[208];
diff -r 563e3c3cf26a -r 5b0baf874eb9 mpz/jacobi.c
--- a/mpz/jacobi.c Mon May 10 11:22:02 2010 +0200
+++ b/mpz/jacobi.c Tue May 11 16:09:10 2010 +0200
@@ -1,6 +1,6 @@
/* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
-Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
+Copyright 2000, 2001, 2002, 2005, 2010 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -23,22 +23,6 @@
#include "longlong.h"
-/* Change this to "#define TRACE(x) x" for some traces. */
-#define TRACE(x)
-
-
-#define MPN_RSHIFT_OR_COPY(dst,src,size,shift) \
- do { \
- if ((shift) != 0) \
- { \
- ASSERT_NOCARRY (mpn_rshift (dst, src, size, shift)); \
- (size) -= ((dst)[(size)-1] == 0); \
- } \
- else \
- MPN_COPY (dst, src, size); \
- } while (0)
-
-
/* This code does triple duty as mpz_jacobi, mpz_legendre and
mpz_kronecker. For ABI compatibility, the link symbol is
__gmpz_jacobi, not __gmpz_kronecker, even though the latter would
@@ -54,10 +38,6 @@
multiple of b), but the checking for that takes little time compared to
other operations.
- FIXME: The old code had special cases for a or b fitting in one
- limb, let mod_1 or modexact_1 get used, without any copying, and end
- up just as efficient as the mixed precision mpz_kronecker_ui etc.
-
Enhancements:
mpn_bdiv_qr should be used instead of mpn_tdiv_qr.
@@ -66,18 +46,35 @@
Lehmer. It could stright-forwardly be made subquadratic by using
hgcd in the same way as mpn_gcd. */
-#if 0
+/* (a/2) = -1 iff a = 3 or a = -3 (mod 8), and (2/b) = -1 iff b = 3 or
+ b = - 3 (mod 8). Note that when this is used, we have already
+ excluded the case that a and both have a common factor of two. */
+
+#define STRIP_TWOS(bit1, twos, other_low, p, n, low) \
+ do { \
+ JACOBI_STRIP_LOW_ZEROS ((bit1), (other_low), (p), (n), (low)); \
+ count_trailing_zeros ((twos), (low)); \
+ (bit1) ^= JACOBI_TWOS_U_BIT1((twos), (other_low)); \
+ (low) >>= (twos); \
+ if ((n) > 1 && (twos) > 0) \
+ { \
+ mp_limb_t __second = (p)[1]; \
+ (low) |= __second << (GMP_NUMB_BITS - (twos)); \
+ if ((n) == 2 && (__second >> (twos) == 0)) \
+ n = 1; \
+ } \
+ } while(0)
+
int
mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
{
mp_srcptr asrcp, bsrcp;
- mp_size_t asize, bsize;
+ mp_size_t asize, bsize, itch;
mp_limb_t alow, blow;
- mp_ptr ap, bp;
- int shift;
- unsigned bits;
- unsigned a3;
- int res;
+ mp_ptr ap, bp, scratch;
+ unsigned atwos, btwos;
+ int result_bit1;
+ int res;
TMP_DECL;
asize = SIZ(a);
@@ -88,339 +85,157 @@
bsrcp = PTR(b);
blow = bsrcp[0];
- /* The MPN jacobi function needs positive a and b, and b odd. So we
- * need to handle the cases of a or b zero, then signs, and then the
- * case of even b. */
- if (bsize == 0)
- /* (a/0) = [ a = 1 or a = -1 ] */
- return ABS (asize) == 1 && alow == 1;
+ /* The MPN jacobi functions requies positive a and b, and b odd. So
+ we must to handle the cases of a or b zero, then signs, and then
+ the case of even b.
+
+ In addition, to reduce the number of cases, we arrange so that a
+ is odd, and asize >= bsize. */
if ( (((alow | blow) & 1) == 0))
/* Common factor of 2 ==> (a/b) = 0 */
return 0;
+ if (bsize == 0)
+ /* (a/0) = [ a = 1 or a = -1 ] */
+ return JACOBI_LS0 (alow, asize);
+
if (asize == 0)
/* (0/b) = [ b = 1 or b = - 1 ] */
- return ABS (bsize) == 1 && blow == 1;
+ return JACOBI_0LS (blow, bsize);
- /* (a/-1) = -1 if a < 0, +1 if a >= 0 */
- bits = (asize < 0) && (bsize < 0);
+ if (bsize < 0)
+ {
+ /* (a/-1) = -1 if a < 0, +1 if a >= 0 */
+ result_bit1 = (asize < 0) << 1;
+ bsize = -bsize;
+ }
+ else
+ result_bit1 = 0;
- bsize = ABS (bsize);
+ STRIP_TWOS (result_bit1, btwos, alow, bsrcp, bsize, blow);
- /* (a/2) = -1 iff a = 3 or a = -3 mod 8 */
- a3 = (alow >> 1) ^ (alow >> 2);
+ if (asize < 0)
+ {
+ /* (-1/b) = -1 iff b = 3 (mod 4) */
+ result_bit1 ^= JACOBI_N1B_BIT1(blow);
+ asize = -asize;
+ }
+
+ STRIP_TWOS(result_bit1, atwos, blow, asrcp, asize, alow);
- while (blow == 0)
+ /* Both numbers odd, so arrange so that asize >= bsize */
+ if (asize < bsize)
{
- if (GMP_NUMB_BITS & 1)
- bits ^= a3;
- bsize--;
- blow = *++bsrcp;
+ unsigned t;
+ MPN_SRCPTR_SWAP (asrcp, asize, bsrcp, bsize);
+ MP_LIMB_T_SWAP (alow, blow);
+
+ t = atwos;
+ atwos = btwos;
+ btwos = t;
+
+ result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
+ }
+
+ if (bsize == 1)
+ {
+ if (blow == 1)
+ return JACOBI_BIT1_TO_PN (result_bit1);
+
+ if (asize > 1)
+ {
+ /* We work with {asrcp, asize} mod b, hence throw away the
+ old alow and undo the shift right by atwos. */
+ result_bit1 ^= JACOBI_TWOS_U_BIT1 (atwos, blow);
+
+ JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
+ }
+
+ return mpn_jacobi_base (alow, blow, result_bit1);
}
TMP_MARK;
- bp = TMP_ALLOC_LIMBS (bsize);
- if (blow & 1)
- MPN_COPY (bp, bsrcp, bsize);
- else
- {
- count_trailing_zeros (shift, blow);
- ASSERT_NOCARRY (mpn_rshift (bp, bsrcp, bsize, shift));
- bsize -= (bp[bsize-1] == 0);
- blow = bp[0];
-
- bits ^= shift & a3;
- }
-
- if (asize < 0)
- /* (-1/b) = -1 iff b = 3 mod 4 */
- bits ^= (blow >> 1) & 1;
-
- asize = ABS(asize);
-
- /* FIXME: Take out powers of two in a? */
-
- bits = mpn_jacobi_init (alow, blow, bits & 1);
+ itch = 3*bsize;
if (asize > bsize)
{
- mp_size_t itch = bsize;
- mp_ptr scratch;
+ if (btwos > 0)
+ {
+ if (asize >= 2 * bsize)
+ itch = asize + bsize + 1;
+ }
+ else if (atwos > 0)
+ {
+ if (asize >= 2*bsize)
+ itch = 2*asize - bsize + 1;
+ }
+ else
+ {
+ if (asize >= 3*bsize)
+ itch = asize + 1;
+ }
+ }
- if (asize >= 2*bsize)
- itch = asize - bsize + 1;
+ ap = TMP_ALLOC_LIMBS (itch);
+ bp = ap + bsize;
+ scratch = bp + bsize;
- scratch = TMP_ALLOC_LIMBS (itch);
- ap = TMP_ALLOC_LIMBS (bsize);
+ if (asize > bsize)
+ {
+ /* Do an initial divide. */
+ if (btwos > 0)
+ {
+ /* Result size: 2*bsize, extra: asize - bsize + 1 for
+ quotient, total: asize + bsize + 1 */
+ ASSERT (atwos == 0);
- mpn_tdiv_qr (scratch, ap, 0, asrcp, asize, bp, bsize);
- bits = mpn_jacobi_update (bits, 1, scratch[0] & 3);
- res = mpn_jacobi_lehmer (ap, bp, bsize, bits, scratch);
- }
- else if (asize < bsize)
- {
- mp_size_t itch = 2*asize;
- mp_ptr scratch;
+ ASSERT_NOCARRY (mpn_rshift (bp, bsrcp, bsize, btwos));
More information about the gmp-commit
mailing list