[Gmp-commit] /var/hg/gmp: mini-gmp: New function mpz_probab_prime_p.
mercurial at gmplib.org
mercurial at gmplib.org
Thu Mar 6 17:44:37 UTC 2014
details: /var/hg/gmp/rev/8d71961d8ad2
changeset: 16325:8d71961d8ad2
user: Niels M?ller <nisse at lysator.liu.se>
date: Thu Mar 06 18:44:18 2014 +0100
description:
mini-gmp: New function mpz_probab_prime_p.
diffstat:
ChangeLog | 8 +
mini-gmp/mini-gmp.c | 93 ++++++++++++++++++++++
mini-gmp/mini-gmp.h | 3 +
mini-gmp/tests/Makefile | 2 +-
mini-gmp/tests/t-pprime_p.c | 181 ++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 286 insertions(+), 1 deletions(-)
diffs (truncated from 328 to 300 lines):
diff -r 967e17283e82 -r 8d71961d8ad2 ChangeLog
--- a/ChangeLog Mon Mar 03 20:33:06 2014 +0100
+++ b/ChangeLog Thu Mar 06 18:44:18 2014 +0100
@@ -1,3 +1,11 @@
+2014-03-06 Niels Möller <nisse at lysator.liu.se>
+
+ * mini-gmp/mini-gmp.c (gmp_millerrabin): New internal function.
+ (mpz_probab_prime_p): New function.
+ * mini-gmp/mini-gmp.h (mpz_probab_prime_p): Declare it.
+ * mini-gmp/tests/t-pprime_p.c: New test program.
+ * mini-gmp/tests/Makefile (CHECK_PROGRAMS): Added t-pprime_p.
+
2014-03-03 Niels Möller <nisse at lysator.liu.se>
* mini-gmp/mini-gmp.c (mpz_congruent_p): New function.
diff -r 967e17283e82 -r 8d71961d8ad2 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 Thu Mar 06 18:44:18 2014 +0100
@@ -3358,6 +3358,99 @@
}
+/* 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, int reps)
+{
+ mpz_t nm1;
+ mpz_t q;
+ mpz_t y;
+ mp_bitcnt_t k;
+ int is_prime;
+ int j;
+
+ /* Note that we use the absolute value of n only, for compatibility
+ with the real GMP. */
+ 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_abs (nm1, n);
+ mpz_sub_ui (nm1, nm1, 1);
+ k = mpz_scan1 (nm1, 0);
+ mpz_tdiv_q_2exp (q, nm1, k);
+
+ for (j = 0, is_prime = 1; is_prime && 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;
+ }
+ is_prime &= gmp_millerrabin (n, nm1, y, q, k);
+ }
+ 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 -r 8d71961d8ad2 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 Thu Mar 06 18:44:18 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 -r 8d71961d8ad2 mini-gmp/tests/Makefile
--- a/mini-gmp/tests/Makefile Mon Mar 03 20:33:06 2014 +0100
+++ b/mini-gmp/tests/Makefile Thu Mar 06 18:44:18 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 -r 8d71961d8ad2 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 Thu Mar 06 18:44:18 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"
More information about the gmp-commit
mailing list