[Gmp-commit] /var/hg/gmp: 3 new changesets

mercurial at gmplib.org mercurial at gmplib.org
Mon Nov 12 17:37:26 UTC 2018


details:   /var/hg/gmp/rev/1dd99c5d0bf6
changeset: 17691:1dd99c5d0bf6
user:      "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date:      Mon Nov 12 18:15:06 2018 +0100
description:
mini-gmp/mini-gmp.c: Implement BPSW test for primality.

details:   /var/hg/gmp/rev/74a56b8e46f5
changeset: 17692:74a56b8e46f5
user:      "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date:      Mon Nov 12 18:21:26 2018 +0100
description:
mini-gmp/tests/t-lucm.c: New test (and activate it)

details:   /var/hg/gmp/rev/9d9ff751b84a
changeset: 17693:9d9ff751b84a
user:      "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date:      Mon Nov 12 18:33:55 2018 +0100
description:
ChangeLog

diffstat:

 ChangeLog                    |    4 +
 mini-gmp/ChangeLog           |   11 ++
 mini-gmp/mini-gmp.c          |  196 ++++++++++++++++++++++++++++++++++++++++--
 mini-gmp/tests/Makefile      |    2 +-
 mini-gmp/tests/hex-random.c  |   48 ++++++++++
 mini-gmp/tests/hex-random.h  |    3 +
 mini-gmp/tests/mini-random.c |   18 +++
 mini-gmp/tests/mini-random.h |    2 +
 mini-gmp/tests/t-lucm.c      |   95 ++++++++++++++++++++
 9 files changed, 368 insertions(+), 11 deletions(-)

diffs (truncated from 467 to 300 lines):

diff -r 4c54563d50b0 -r 9d9ff751b84a ChangeLog
--- a/ChangeLog	Fri Aug 17 12:29:02 2018 +0200
+++ b/ChangeLog	Mon Nov 12 18:33:55 2018 +0100
@@ -12,6 +12,10 @@
 	* tests/mpn/t-fib2m.c: New file, tests for mpn_fib2m.
 	* tests/mpn/Makefile.am (check_PROGRAMS): Add t-fib2m.
 
+	* mpn/generic/mod_34lsub1.c: Initialise c[012] once.
+	* tests/mpz/t-pprime_p.c (check_primes): Two more primes.
+	* tests/mp?: Use TESTS_REPS in many files.
+
 2018-11-07  Torbjörn Granlund  <tg at gmplib.org>
 
 	* configure.ac (arm): Support a12 and a17.
diff -r 4c54563d50b0 -r 9d9ff751b84a mini-gmp/ChangeLog
--- a/mini-gmp/ChangeLog	Fri Aug 17 12:29:02 2018 +0200
+++ b/mini-gmp/ChangeLog	Mon Nov 12 18:33:55 2018 +0100
@@ -1,3 +1,14 @@
+2018-10-30 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+	* mini-gmp.c (mpz_probab_prime_p): BPSW test for primality.
+
+	* tests/hex-random.c (hex_random_lucm_op): New function.
+	* tests/hex-random.h: Declare it.
+	* tests/mini-random.c (mini_random_lucm_op): New function.
+	* tests/mini-random.h: Declare it.
+	* mini-gmp/tests/t-lucm.c: New test
+	* mini-gmp/tests/Makefile (CHECK_PROGRAMS): Add t-lucm.
+
 2018-09-07  Niels Möller  <nisse at lysator.liu.se>
 
 	* tests/t-div.c (testmain): Add missing const declarations.
diff -r 4c54563d50b0 -r 9d9ff751b84a mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c	Fri Aug 17 12:29:02 2018 +0200
+++ b/mini-gmp/mini-gmp.c	Mon Nov 12 18:33:55 2018 +0100
@@ -3409,6 +3409,177 @@
 
 

 /* Primality testing */
+
+/* Computes Kronecker (a/b) with odd b, a!=0 and GCD(a,b) = 1 */
+/* Adapted from JACOBI_BASE_METHOD==4 in mpn/generic/jacbase.c */
+static int
+gmp_jacobi_coprime (mp_limb_t a, mp_limb_t b)
+{
+  int c, bit = 0;
+
+  assert (b & 1);
+  assert (a != 0);
+  /* assert (mpn_gcd_11 (a, b) == 1); */
+
+  /* Below, we represent a and b shifted right so that the least
+     significant one bit is implicit. */
+  b >>= 1;
+
+  gmp_ctz(c, a);
+  a >>= 1;
+
+  do
+    {
+      a >>= c;
+      /* (2/b) = -1 if b = 3 or 5 mod 8 */
+      bit ^= c & (b ^ (b >> 1));
+      if (a < b)
+	{
+	  bit ^= a & b;
+	  a = b - a;
+	  b -= a;
+	}
+      else
+	{
+	  a -= b;
+	  assert (a != 0);
+	}
+
+      gmp_ctz(c, a);
+      ++c;
+    }
+  while (b > 0);
+
+  return bit & 1 ? -1 : 1;
+}
+
+static void
+gmp_lucas_step_k_2k (mpz_t V, mpz_t Qk, const mpz_t n)
+{
+  mpz_mod (Qk, Qk, n);
+  /* V_{2k} <- V_k ^ 2 - 2Q^k */
+  mpz_mul (V, V, V);
+  mpz_submul_ui (V, Qk, 2);
+  mpz_tdiv_r (V, V, n);
+  /* Q^{2k} = (Q^k)^2 */
+  mpz_mul (Qk, Qk, Qk);
+}
+
+/* Computes V_k, Q^k (mod n) for the Lucas' sequence */
+/* with P=1, Q=Q; k = (n>>b0)|1. */
+/* Requires an odd n > 4; b0 > 0; -2*Q must not overflow a long */
+/* Returns (U_k == 0) and sets V=V_k and Qk=Q^k. */
+int
+mpz_lucas_mod (mpz_t V, mpz_t Qk, long Q,
+	       mp_bitcnt_t b0, const mpz_t n)
+{
+  mp_bitcnt_t bs;
+  mpz_t U;
+  int res;
+
+  assert (b0 > 0);
+  assert (Q <= (long) (GMP_LIMB_HIGHBIT >> 1));
+  assert (Q > -(long) (GMP_LIMB_HIGHBIT >> 1));
+  assert (mpz_cmp_ui (n, 4) > 0);
+  assert (mpz_odd_p (n));
+
+  mpz_init_set_ui (U, 1); /* U1 = 1 */
+  mpz_set_ui (V, 1); /* V1 = 1 */
+  mpz_set_si (Qk, Q);
+
+  for (bs = mpz_sizeinbase (n, 2) - 1; --bs >= b0;)
+    {
+      /* U_{2k} <- U_k * V_k */
+      mpz_mul (U, U, V);
+      /* V_{2k} <- V_k ^ 2 - 2Q^k */
+      /* Q^{2k} = (Q^k)^2 */
+      gmp_lucas_step_k_2k (V, Qk, n);
+
+      /* A step k->k+1 is performed if the bit in $n$ is 1	*/
+      /* mpz_tstbit(n,bs) or the the bit is 0 in $n$ but	*/
+      /* should be 1 in $n+1$ (bs == b0)			*/
+      if (b0 == bs || mpz_tstbit (n, bs))
+	{
+	  /* Q^{k+1} <- Q^k * Q */
+	  mpz_mul_si (Qk, Qk, Q);
+	  /* U_{k+1} <- (U_k + V_k) / 2 */
+	  mpz_swap (U, V); /* Keep in V the old value of U_k */
+	  mpz_add (U, U, V);
+	  /* We have to compute U/2, so we need an even value, */
+	  /* equivalent (mod n) */
+	  if (mpz_odd_p (U))
+	    mpz_add (U, U, n);
+	  mpz_tdiv_q_2exp (U, U, 1);
+	  /* V_{k+1} <-(D*U_k + V_k) / 2 =
+			U_{k+1} + (D-1)/2*U_k = U_{k+1} - 2Q*U_k */
+	  mpz_mul_si (V, V, -2*Q);
+	  mpz_add (V, U, V);
+	  mpz_tdiv_r (V, V, n);
+	}
+      mpz_tdiv_r (U, U, n);
+    }
+
+  res = U->_mp_size == 0;
+  mpz_clear (U);
+  return res;
+}
+
+/* Performs strong Lucas' test on x, with parameters suggested */
+/* for the BPSW test. Qk is only passed to recycle a variable. */
+/* Requires GCD (x,6) = 1.*/
+static int
+gmp_stronglucas (const mpz_t x, mpz_t Qk)
+{
+  mp_bitcnt_t b0;
+  mpz_t V, n;
+  mp_limb_t maxD, D; /* The absolute value is stored. */
+  long Q;
+  mp_limb_t tl;
+
+  /* Test on the absolute value. */
+  mpz_roinit_normal_n (n, x->_mp_d, GMP_ABS (x->_mp_size));
+
+  assert (mpz_odd_p (n));
+  /* assert (mpz_gcd_ui (NULL, n, 6) == 1); */
+  if (mpz_root (Qk, n, 2))
+    return 0; /* A square is composite. */
+
+  /* Check Ds up to square root (in case, n is prime)
+     or avoid overflows */
+  maxD = (Qk->_mp_size == 1) ? Qk->_mp_d [0] - 1 : GMP_LIMB_MAX;
+
+  D = 3;
+  /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
+  /* For those Ds we have (D/n) = (n/|D|) */
+  do
+    {
+      if (D >= maxD)
+	return 1 + (D != GMP_LIMB_MAX); /* (1 + ! ~ D) */
+      D += 2;
+      tl = mpz_tdiv_ui (n, D);
+      if (tl == 0)
+	return 0;
+    }
+  while (gmp_jacobi_coprime (tl, D) == 1);
+
+  mpz_init (V);
+
+  /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
+  b0 = mpz_scan0 (n, 0);
+
+  /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
+  Q = (D & 2) ? (D >> 2) + 1 : -(D >> 2);
+
+  if (! mpz_lucas_mod (V, Qk, Q, b0, n))	/* If Ud != 0 */
+    while (V->_mp_size != 0 && --b0 != 0)	/* while Vk != 0 */
+      /* V <- V ^ 2 - 2Q^k */
+      /* Q^{2k} = (Q^k)^2 */
+      gmp_lucas_step_k_2k (V, Qk, n);
+
+  mpz_clear (V);
+  return (b0 != 0);
+}
+
 static int
 gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
 		 const mpz_t q, mp_bitcnt_t k)
@@ -3470,21 +3641,26 @@
   if (mpz_cmpabs_ui (n, 31*31) < 0)
     return 2;
 
+  mpz_init (nm1);
+  mpz_init (q);
+
+  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
+  mpz_abs (nm1, n);
+  nm1->_mp_d[0] -= 1;
+  k = mpz_scan1 (nm1, 0);
+  mpz_tdiv_q_2exp (q, nm1, k);
+
+  /* BPSW test */
+  mpz_init_set_ui (y, 2);
+  is_prime = gmp_millerrabin (n, nm1, y, q, k) && gmp_stronglucas (n, y);
+  reps -= 25; /* skip the first 25 repetitions */
+
   /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
      j^2 + j + 41 using Euler's polynomial. We potentially stop early,
      if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
      30 (a[30] == 971 > 31*31 == 961). */
 
-  mpz_init (nm1);
-  mpz_init (q);
-  mpz_init (y);
-
-  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
-  nm1->_mp_size = mpz_abs_sub_ui (nm1, n, 1);
-  k = mpz_scan1 (nm1, 0);
-  mpz_tdiv_q_2exp (q, nm1, k);
-
-  for (j = 0, is_prime = 1; is_prime & (j < reps); j++)
+  for (j = 0; is_prime & (j < reps); j++)
     {
       mpz_set_ui (y, (unsigned long) j*j+j+41);
       if (mpz_cmp (y, nm1) >= 0)
diff -r 4c54563d50b0 -r 9d9ff751b84a mini-gmp/tests/Makefile
--- a/mini-gmp/tests/Makefile	Fri Aug 17 12:29:02 2018 +0200
+++ b/mini-gmp/tests/Makefile	Mon Nov 12 18:33:55 2018 +0100
@@ -30,7 +30,7 @@
 CHECK_PROGRAMS = t-add t-sub t-mul t-invert t-div t-div_2exp \
 	t-double t-cmp_d t-gcd t-lcm t-import t-comb t-signed \
 	t-sqrt t-root t-powm t-logops t-bitops t-scan t-str \
-	t-reuse t-aorsmul t-limbs t-cong t-pprime_p \
+	t-reuse t-aorsmul t-limbs t-cong t-pprime_p t-lucm \
 	t-mpq_addsub t-mpq_muldiv t-mpq_muldiv_2exp
 # Default TESTS to all tests, allowing overriding TESTS for building tests
 # without running them.
diff -r 4c54563d50b0 -r 9d9ff751b84a mini-gmp/tests/hex-random.c
--- a/mini-gmp/tests/hex-random.c	Fri Aug 17 12:29:02 2018 +0200
+++ b/mini-gmp/tests/hex-random.c	Mon Nov 12 18:33:55 2018 +0100
@@ -28,6 +28,8 @@
 #endif
 
 #include "gmp.h"
+/* FIXME: gmp-impl.h included only for mpz_lucas_mod */
+#include "gmp-impl.h"
 
 #include "hex-random.h"
 
@@ -481,3 +483,49 @@
 
   mpz_clear (a);
 }
+
+void hex_random_lucm_op (unsigned long maxbits,
+			 char **vp, char **qp, char **mp,
+			 long *Q, unsigned long *b0, int *res)
+{
+  mpz_t m, v, q, t1, t2;
+  unsigned long mbits;
+
+  mpz_init (m);
+  mpz_init (v);
+  mpz_init (q);
+  mpz_init (t1);
+  mpz_init (t2);
+
+  *Q = gmp_urandomb_ui (state, 14) + 1;
+
+  do
+    {
+      mbits = gmp_urandomb_ui (state, 32) % maxbits + 5;
+
+      mpz_rrandomb (m, state, mbits);
+      *b0 = gmp_urandomb_ui (state, 32) % (mbits - 3) + 2;


More information about the gmp-commit mailing list