[Gmp-commit] /home/hgfiles/gmp: New test with large quotients.
mercurial at gmplib.org
mercurial at gmplib.org
Fri Feb 26 14:03:09 CET 2010
details: /home/hgfiles/gmp/rev/72bd403bf8b8
changeset: 13469:72bd403bf8b8
user: Niels M?ller <nisse at lysator.liu.se>
date: Fri Feb 26 14:01:09 2010 +0100
description:
New test with large quotients.
diffstat:
ChangeLog | 6 +
tests/mpz/t-jac.c | 222 +++++++++++++++++++++++++++++++++++++++++++++++++++++-
2 files changed, 227 insertions(+), 1 deletions(-)
diffs (259 lines):
diff -r 0ea6d9d8163d -r 72bd403bf8b8 ChangeLog
--- a/ChangeLog Fri Feb 26 10:32:50 2010 +0100
+++ b/ChangeLog Fri Feb 26 14:01:09 2010 +0100
@@ -1,3 +1,9 @@
+2010-02-26 Niels Möller <nisse at lysator.liu.se>
+
+ * tests/mpz/t-jac.c (check_large_quotients): New test. Currently
+ disabled, since it's quite slow.
+ (mpz_nextprime_step): New function.
+
2010-02-26 Torbjorn Granlund <tege at gmplib.org>
* mpn/pa64/aors_n.asm: Fix typo in last change.
diff -r 0ea6d9d8163d -r 72bd403bf8b8 tests/mpz/t-jac.c
--- a/tests/mpz/t-jac.c Fri Feb 26 10:32:50 2010 +0100
+++ b/tests/mpz/t-jac.c Fri Feb 26 14:01:09 2010 +0100
@@ -41,6 +41,8 @@
#include "gmp-impl.h"
#include "tests.h"
+/* For count_leading_zeros in mpz_nextprime_step. */
+#include "longlong.h"
#ifdef _LONG_LONG_LIMB
#define LL(l,ll) ll
@@ -814,6 +816,222 @@
#undef PRIME_MAX_B_SIZE
}
+static const unsigned char primegap[] =
+{
+ 2,2,4,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,14,4,6,
+ 2,10,2,6,6,4,6,6,2,10,2,4,2,12,12,4,2,4,6,2,10,6,6,6,2,6,4,2,10,14,4,2,
+ 4,14,6,10,2,4,6,8,6,6,4,6,8,4,8,10,2,10,2,6,4,6,8,4,2,4,12,8,4,8,4,6,
+ 12,2,18,6,10,6,6,2,6,10,6,6,2,6,6,4,2,12,10,2,4,6,6,2,12,4,6,8,10,8,10,8,
+ 6,6,4,8,6,4,8,4,14,10,12,2,10,2,4,2,10,14,4,2,4,14,4,2,4,20,4,8,10,8,4,6,
+ 6,14,4,6,6,8,6,12
+};
+
+#define NUMBER_OF_PRIMES 167
+
+/* Similar to mpz_nextprime, finds the first (odd) prime of the form n
+ + k * step, with k >= 1. If n and step has a cmmon factor, it never
+ temrinates... */
+static void
+mpz_nextprime_step (mpz_ptr p, mpz_srcptr n, mpz_srcptr step_in)
+{
+ unsigned short *moduli;
+ unsigned short *step_moduli;
+ unsigned long difference;
+ int i;
+ unsigned prime_limit;
+ unsigned long prime;
+ int cnt;
+ mp_size_t pn;
+ mp_bitcnt_t nbits;
+ unsigned incr;
+ mpz_t step;
+ TMP_SDECL;
+
+ mpz_init (step);
+
+ switch ( (mpz_odd_p (n) << 1) + mpz_odd_p (step_in))
+ {
+ default:
+ case 0:
+ /* Both even. */
+ abort ();
+ case 1:
+ /* n even, step odd. Use odd k. */
+ mpz_mul_2exp (step, step_in, 1);
+ mpz_add (p, n, step_in);
+ break;
+ case 2:
+ /* n odd, step even. All k > 0 give odd result. */
+ mpz_set (step, step_in);
+ mpz_add (p, n, step_in);
+ break;
+ case 3:
+ /* Both n and step odd. Use even k. */
+ mpz_mul_2exp (step, step_in, 1);
+ mpz_add (p, n, step);
+ break;
+ }
+
+ ASSERT_ALWAYS (mpz_odd_p (p));
+ ASSERT_ALWAYS (mpz_even_p (step));
+
+ if (mpz_cmp_ui (p, 7) <= 0)
+ {
+ mpz_clear (step);
+ return;
+ }
+
+ pn = SIZ(p);
+ count_leading_zeros (cnt, PTR(p)[pn - 1]);
+ nbits = pn * GMP_NUMB_BITS - (cnt - GMP_NAIL_BITS);
+ if (nbits / 2 >= NUMBER_OF_PRIMES)
+ prime_limit = NUMBER_OF_PRIMES - 1;
+ else
+ prime_limit = nbits / 2;
+
+ TMP_SMARK;
+
+ /* Compute residues modulo small odd primes */
+ moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short);
+ step_moduli = TMP_SALLOC_TYPE (prime_limit * sizeof step_moduli[0], unsigned short);
+
+ for (;;)
+ {
+ ASSERT_ALWAYS (mpz_odd_p (p));
+
+ /* FIXME: Compute lazily? */
+ prime = 3;
+ for (i = 0; i < prime_limit; i++)
+ {
+ moduli[i] = mpz_fdiv_ui (p, prime);
+ step_moduli[i] = mpz_fdiv_ui (step, prime);
+ prime += primegap[i];
+ }
+
+ /* INCR_LIMIT * (max_prime - 1) must fit in an unsigned. */
+#define INCR_LIMIT 0x10000
+
+ for (difference = incr = 0; incr < INCR_LIMIT; difference ++)
+ {
+ /* First check residues */
+ prime = 3;
+ for (i = 0; i < prime_limit; i++)
+ {
+ unsigned r;
+ /* FIXME: Reduce moduli + incr and store back, to allow
+ for division-free reductions. Alternatively, table
+ primes[]'s inverses. */
+ r = (moduli[i] + incr*step_moduli[i]) % prime;
+ prime += primegap[i];
+
+ if (r == 0)
+ goto next;
+ }
+
+ mpz_addmul_ui (p, step, difference);
+ difference = 0;
+
+ ASSERT_ALWAYS (mpz_odd_p (p));
+
+ /* Miller-Rabin test */
+ if (mpz_millerrabin (p, 10))
+ goto done;
+ next:;
+ incr ++;
+ }
+
+ mpz_addmul_ui (p, step, difference);
+ difference = 0;
+ }
+ done:
+ mpz_clear (step);
+ TMP_SFREE;
+}
+
+void
+check_large_quotients (void)
+{
+#define COUNT 3
+#define MAX_THRESHOLD 30
+
+ gmp_randstate_ptr rands = RANDS;
+ unsigned i;
+ mpz_t op1, op2, temp1, temp2, bs;
+
+ mpz_init (op1);
+ mpz_init (op2);
+ mpz_init (temp1);
+ mpz_init (temp2);
+ mpz_init (bs);
+
+ for (i = 0; i < COUNT; i++)
+ {
+ unsigned j;
+ unsigned chain_len;
+ int answer;
+
+ /* Code copied from t-gcd.c */
+ mpz_set_ui (op1, 0);
+ mpz_set_ui (op2, 1);
+
+ mpz_urandomb (bs, rands, 32);
+ chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * MAX_THRESHOLD / 256);
+
+ for (j = 0; j < chain_len; j++)
+ {
+ mpz_urandomb (bs, rands, 32);
+ mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
+ mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
+ mpz_add_ui (temp2, temp2, 1);
+ mpz_mul (temp1, op2, temp2);
+ mpz_add (op1, op1, temp1);
+
+ /* Don't generate overly huge operands. */
+ if (SIZ (op1) > 3 * MAX_THRESHOLD)
+ {
+ mpz_swap (op1, op2);
+ break;
+ }
+
+ mpz_urandomb (bs, rands, 32);
+ mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
+ mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
+ mpz_add_ui (temp2, temp2, 1);
+ mpz_mul (temp1, op1, temp2);
+ mpz_add (op2, op2, temp1);
+
+ /* Don't generate overly huge operands. */
+ if (SIZ (op2) > 3 * MAX_THRESHOLD)
+ break;
+ }
+ ASSERT_ALWAYS (mpz_cmp (op1, op2) < 0);
+
+ if (mpz_odd_p (op1) && mpz_probab_prime_p (op1, 5))
+ {
+ answer = refmpz_legendre (op2, op1);
+ try_all (op2, op1, answer);
+ }
+ else if (mpz_odd_p (op2) && mpz_probab_prime_p (op2, 5))
+ {
+ answer = refmpz_legendre (op1, op2);
+ try_all (op1, op2, answer);
+ }
+ else
+ {
+ mpz_nextprime_step (op1, op2, op1);
+ answer = refmpz_legendre (op2, op1);
+ try_all (op2, op1, answer);
+ }
+ }
+ mpz_clear (op1);
+ mpz_clear (op2);
+ mpz_clear (temp1);
+ mpz_clear (temp2);
+ mpz_clear (bs);
+#undef COUNT
+#undef MAX_THRESHOLD
+}
+
int
main (int argc, char *argv[])
{
@@ -837,7 +1055,9 @@
check_squares_zi ();
check_a_zero ();
check_jacobi_factored ();
-
+#if 0
+ check_large_quotients ();
+#endif
tests_end ();
exit (0);
}
More information about the gmp-commit
mailing list