mini-gmp: mpz_probab_prime_p

Niels Möller nisse at lysator.liu.se
Wed Mar 5 21:47:10 UTC 2014


bodrato at mail.dm.unipi.it writes:

> Then I suggest to simplify the first few lines of prob_prime with:

Below is a revised implementation, including some test cases (similar
testcases could be added to tests/mpz/t-pprime_p.c too). I think I'll
check this in within a day or two.

Regards,
/Niels


diff -r 967e17283e82 mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c	Mon Mar 03 20:33:06 2014 +0100
+++ b/mini-gmp/mini-gmp.c	Wed Mar 05 22:43:10 2014 +0100
@@ -3358,6 +3358,104 @@
 }
 
 

+/* Primality testing */
+static int
+gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
+		 const mpz_t q, mp_bitcnt_t k)
+{
+  mp_bitcnt_t i;
+
+  /* Caller must initialize y to the base. */
+  mpz_powm (y, y, q, n);
+
+  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
+    return 1;
+
+  for (i = 1; i < k; i++)
+    {
+      mpz_powm_ui (y, y, 2, n);
+      if (mpz_cmp (y, nm1) == 0)
+	return 1;
+      if (mpz_cmp_ui (y, 1) == 0)
+	return 0;
+    }
+  return 0;
+}
+
+/* This product is 0xc0cfd797, and fits in 32 bits. */
+#define GMP_PRIME_PRODUCT \
+  (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
+
+/* Bit (p+1)/2 is set, for each odd prime <= 61 */
+#define GMP_PRIME_MASK 0xc96996dcUL
+
+int
+mpz_probab_prime_p (const mpz_t n_in, int reps)
+{
+  mpz_t n;
+  mpz_t nm1;
+  mpz_t q;
+  mpz_t y;
+  mp_bitcnt_t k;
+  int is_prime;
+  int j;
+
+  /* Use absolute value only, for compatibility with the real GMP. */
+  mpz_roinit_n (n, n_in->_mp_d, GMP_ABS (n_in->_mp_size));
+
+  if (mpz_even_p (n))
+    return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
+
+  /* Above test excludes n == 0 */
+  assert (n->_mp_size > 0);
+  
+  if (mpz_cmpabs_ui (n, 64) < 0)
+    return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
+
+  if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
+    return 0;
+
+  /* All prime factors are >= 31. */
+  if (mpz_cmpabs_ui (n, 31*31) < 0)
+    return 2;
+  
+  /* 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.  */
+  mpz_sub_ui (nm1, n, 1);
+  k = mpz_scan1 (nm1, 0);
+  mpz_tdiv_q_2exp (q, nm1, k);
+  
+  for (j = 0, is_prime = 1; j < reps; j++)
+    {
+      mpz_set_ui (y, (unsigned long) j*j+j+41);
+      if (mpz_cmp (y, nm1) >= 0)
+	{
+	  /* Don't try any further bases. */
+	  assert (j >= 30);
+	  break;
+	}
+      if (!gmp_millerrabin (n, nm1, y, q, k))
+	{
+	  is_prime = 0;
+	  break;
+	}
+    }
+  mpz_clear (nm1);
+  mpz_clear (q);
+  mpz_clear (y);
+    
+  return is_prime;
+}
+
+

 /* Logical operations and bit manipulation. */
 
 /* Numbers are treated as if represented in two's complement (and
diff -r 967e17283e82 mini-gmp/mini-gmp.h
--- a/mini-gmp/mini-gmp.h	Mon Mar 03 20:33:06 2014 +0100
+++ b/mini-gmp/mini-gmp.h	Wed Mar 05 22:43:10 2014 +0100
@@ -215,6 +215,9 @@
 void mpz_fac_ui (mpz_t, unsigned long);
 void mpz_bin_uiui (mpz_t, unsigned long, unsigned long);
 
+int
+mpz_probab_prime_p (const mpz_t, int);
+
 int mpz_tstbit (const mpz_t, mp_bitcnt_t);
 void mpz_setbit (mpz_t, mp_bitcnt_t);
 void mpz_clrbit (mpz_t, mp_bitcnt_t);
diff -r 967e17283e82 mini-gmp/tests/Makefile
--- a/mini-gmp/tests/Makefile	Mon Mar 03 20:33:06 2014 +0100
+++ b/mini-gmp/tests/Makefile	Wed Mar 05 22:43:10 2014 +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-reuse t-aorsmul t-limbs t-cong t-pprime_p
 
 MISC_OBJS = hex-random.o mini-random.o testutils.o
 
diff -r 967e17283e82 mini-gmp/tests/t-pprime_p.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mini-gmp/tests/t-pprime_p.c	Wed Mar 05 22:43:10 2014 +0100
@@ -0,0 +1,181 @@
+/* test mpz_probab_prime_p
+
+Copyright 2001, 2002, 2004, 2011, 2012, 2014 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library test suite.
+
+The GNU MP Library test suite is free software; you can redistribute it
+and/or modify it under the terms of the GNU 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 test suite 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 General
+Public License for more details.
+
+You should have received a copy of the GNU General Public License along with
+the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
+
+#include "testutils.h"
+
+static int
+isprime (unsigned long int t)
+{
+  unsigned long int q, r, d;
+
+  if (t < 32)
+    return (0xa08a28acUL >> t) & 1;
+  if ((t & 1) == 0)
+    return 0;
+
+  if (t % 3 == 0)
+    return 0;
+  if (t % 5 == 0)
+    return 0;
+  if (t % 7 == 0)
+    return 0;
+
+  for (d = 11;;)
+    {
+      q = t / d;
+      r = t - q * d;
+      if (q < d)
+	return 1;
+      if (r == 0)
+	break;
+      d += 2;
+      q = t / d;
+      r = t - q * d;
+      if (q < d)
+	return 1;
+      if (r == 0)
+	break;
+      d += 4;
+    }
+  return 0;
+}
+
+static void
+check_one (mpz_srcptr n, int want)
+{
+  int  got;
+
+  got = mpz_probab_prime_p (n, 25);
+
+  /* "definitely prime" is fine if we only wanted "probably prime" */
+  if (got == 2 && want == 1)
+    want = 2;
+
+  if (got != want)
+    {
+      printf ("mpz_probab_prime_p\n");
+      dump   ("  n    ", n);
+      printf ("  got =%d", got);
+      printf ("  want=%d\n", want);
+      abort ();
+    }
+}
+
+static void
+check_pn (mpz_ptr n, int want)
+{
+  check_one (n, want);
+  mpz_neg (n, n);
+  check_one (n, want);
+}
+
+static void
+check_small (void)
+{
+  mpz_t  n;
+  long   i;
+
+  mpz_init (n);
+
+  for (i = 0; i < 300; i++)
+    {
+      mpz_set_si (n, i);
+      check_pn (n, isprime (i));
+    }
+
+  mpz_clear (n);
+}
+
+void
+check_composites (void)
+{
+  int i;
+  int reps = 1000;
+  mpz_t a, b, n, bs;
+  unsigned long size_range, size;
+
+  mpz_init (a);
+  mpz_init (b);
+  mpz_init (n);
+  mpz_init (bs);
+
+  for (i = 0; i < reps; i++)
+    {
+      mini_urandomb (bs, 32);
+      size_range = mpz_get_ui (bs) % 12 + 1; /* 0..4096 bit operands */
+
+      mini_urandomb (bs, size_range);
+      size = mpz_get_ui (bs);
+      mini_rrandomb (a, size);
+
+      mini_urandomb (bs, 32);
+      size_range = mpz_get_ui (bs) % 12 + 1; /* 0..4096 bit operands */
+
+      /* Exclude trivial factors */
+      if (mpz_cmp_ui (a, 1) == 0)
+	mpz_set_ui (a, 2);
+      if (mpz_cmp_ui (b, 1) == 0)
+	mpz_set_ui (b, 2);
+
+      mpz_mul (n, a, b);
+      
+      check_pn (n, 0);
+    }
+  mpz_clear (a);
+  mpz_clear (b);
+  mpz_clear (n);
+  mpz_clear (bs);
+}
+
+static void
+check_primes (void)
+{
+  static const char * const primes[] = {
+    "2", "17", "65537",
+    /* diffie-hellman-group1-sha1, also "Well known group 2" in RFC
+       2412, 2^1024 - 2^960 - 1 + 2^64 * { [2^894 pi] + 129093 } */  
+    "0xFFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD1"
+    "29024E088A67CC74020BBEA63B139B22514A08798E3404DD"
+    "EF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245"
+    "E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7ED"
+    "EE386BFB5A899FA5AE9F24117C4B1FE649286651ECE65381"
+    "FFFFFFFFFFFFFFFF",
+    NULL
+  };
+
+  mpz_t n;
+  int i;
+
+  mpz_init (n);
+
+  for (i = 0; primes[i]; i++)
+    {
+      mpz_set_str_or_abort (n, primes[i], 0);
+      check_one (n, 1);
+    }
+  mpz_clear (n);  
+}
+
+void
+testmain (int argc, char *argv[])
+{
+  check_small ();
+  check_composites ();
+  check_primes ();
+}

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.


More information about the gmp-devel mailing list