[libmath-prime-util-perl] 03/11: Miller-Rabin and prob_prime
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:16 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.03
in repository libmath-prime-util-perl.
commit f7deef17d9c02fcbd9eea546025cfc8c9a8342e8
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Jun 6 04:10:36 2012 -0600
Miller-Rabin and prob_prime
---
XS.xs | 41 +++++++++++
factor.c | 190 ++++++++++++++++++++++++++++++++++++++++++++-----
factor.h | 34 +++++++++
lib/Math/Prime/Util.pm | 37 +++++++++-
util.c | 11 ++-
5 files changed, 290 insertions(+), 23 deletions(-)
diff --git a/XS.xs b/XS.xs
index 2be5131..fbeb66a 100644
--- a/XS.xs
+++ b/XS.xs
@@ -257,6 +257,18 @@ factor(IN UV n)
/* trial factorization for each item in the list */
UV f, m, limit;
n = factor_list[i];
+ /* We pulled out everything through this point, so must be prime */
+ if (n < (19*19)) {
+ XPUSHs(sv_2mortal(newSVuv( n )));
+ continue;
+ /* If sufficiently large, run a prob prime test on it. This saves
+ * us a huge amount of work on big primes. We could also look into
+ * some possible shortcuts. What this really needs is to move
+ * the prime test up above SQUFOF. */
+ } else if ( (n > 40000000) && is_prob_prime(n) ) {
+ XPUSHs(sv_2mortal(newSVuv( n )));
+ continue;
+ }
f = startfactor;
m = startfactor % 30;
limit = sqrt((double) n);
@@ -347,3 +359,32 @@ pminus1_factor(IN UV n)
for (i = 0; i < nfactors; i++) {
XPUSHs(sv_2mortal(newSVuv( factors[i] )));
}
+
+int
+miller_rabin(IN UV n, ...)
+ PREINIT:
+ UV bases[64];
+ int prob_prime = 1;
+ int c = 1;
+ CODE:
+ if (items < 2)
+ croak("No bases given to miller_rabin");
+ if ( (n == 0) || (n == 1) ) XSRETURN(0); /* 0 and 1 are composite */
+ if ( (n == 2) || (n == 3) ) XSRETURN(2); /* 2 and 3 are prime */
+ while (c < items) {
+ int b = 0;
+ while (c < items) {
+ bases[b++] = SvUV(ST(c));
+ c++;
+ if (b == 64) break;
+ }
+ prob_prime = miller_rabin(n, bases, b);
+ if (prob_prime != 1)
+ break;
+ }
+ RETVAL = prob_prime;
+ OUTPUT:
+ RETVAL
+
+int
+is_prob_prime(IN UV n)
diff --git a/factor.c b/factor.c
index 18a269a..7497f8a 100644
--- a/factor.c
+++ b/factor.c
@@ -8,6 +8,7 @@
#include "factor.h"
#include "util.h"
#include "sieve.h"
+#include "bitarray.h"
/*
* You need to remember to use UV for unsigned and IV for signed types that
@@ -88,22 +89,176 @@ static UV gcd_ui(UV x, UV y) {
return x;
}
+static UV mulmod(UV a, UV b, UV m) {
+ UV p;
+ UV r = 0;
+ while (b > 0) {
+ if (b & 1) {
+ if (r == 0) {
+ r = a;
+ } else {
+ r = m - r;
+ r = (a >= r) ? a-r : m-r+a;
+ }
+ }
+ a = (a > (m-a)) ? (a-m)+a : a+a;
+ b >>= 1;
+ }
+ return r;
+}
+
/* n^power + a mod m */
-/* This has serious overflow issues, making the programs that use it dubious */
-static UV powmod(UV n, UV power, UV add, UV m) {
+static UV powmodadd(UV n, UV power, UV add, UV m) {
+ UV t = 1;
+ while (power) {
+ if (power & 1)
+ t = mulmod(t, n, m);
+ n = mulmod(n, n, m);
+ power >>= 1;
+ }
+ t = ((m-t) > add) ? t+add : t+add-m; /* (t+a) % m where t < m */
+ return t;
+}
+
+/* n^power mod m */
+static UV powmod(UV n, UV power, UV m) {
UV t = 1;
- n = n % m;
- /* t and n will always be < m from now on */
while (power) {
if (power & 1)
- t = (t * n) % m;
- n = (n * n) % m;
+ t = mulmod(t, n, m);
+ n = mulmod(n, n, m);
power >>= 1;
}
- /* (t+a) % m, noting t is always < m */
- return ( ((m-t) > add) ? (t+add) : (t+add-m) );
+ return t;
+}
+
+
+/* Miller-Rabin probabilistic primality test
+ * Returns 1 if probably prime relative to the bases, 0 if composite.
+ * Bases must be between 2 and n-2
+ */
+int miller_rabin(UV n, const UV *bases, UV nbases)
+{
+ int b;
+ int s = 0;
+ UV d = n-1;
+
+ assert(n > 3);
+
+ while ( (d&1) == 0 ) {
+ s++;
+ d >>= 1;
+ }
+ for (b = 0; b < nbases; b++) {
+ int r;
+ UV a = bases[b];
+ UV x;
+
+ /* Skip invalid bases */
+ if ( (a < 2) || (a > (n-2)) )
+ croak("Base %"UVuf" is invalid for input %"UVuf, a, n);
+
+ x = powmod(a, d, n);
+ if ( (x == 1) || (x == (n-1)) ) continue;
+
+ for (r = 0; r < s; r++) {
+ x = powmod(x, 2, n);
+ if (x == 1) {
+ return 0;
+ } else if (x == (n-1)) {
+ a = 0;
+ break;
+ }
+ }
+ if (a != 0)
+ return 0;
+ }
+ return 1;
+}
+
+int is_prob_prime(UV n)
+{
+ UV bases[12];
+ int nbases;
+ int prob_prime;
+
+ if ( (n == 2) || (n == 3) || (n == 5) || (n == 7) )
+ return 2;
+ if ( (n<2) || ((n% 2)==0) || ((n% 3)==0) || ((n% 5)==0) || ((n% 7)==0) )
+ return 0;
+ if (n < 121)
+ return 2;
+
+#if BITS_PER_WORD == 32
+ if (n < UVCONST(9080191)) {
+ bases[0] = 31; bases[1] = 73; nbases = 2;
+ } else {
+ bases[0] = 2; bases[1] = 7; bases[2] = 61; nbases = 3;
+ }
+#else
+#if 1
+ /* Better basis from: http://miller-rabin.appspot.com/ */
+ if (n < UVCONST(9080191)) {
+ bases[0] = 31; bases[1] = 73; nbases = 2;
+ } else if (n < UVCONST(4759123141)) {
+ bases[0] = 2; bases[1] = 7; bases[2] = 61; nbases = 3;
+ } else if (n < UVCONST(105936894253)) {
+ bases[0] = 2;
+ bases[1] = 1005905886;
+ bases[2] = 1340600841;
+ nbases = 3;
+ } else if (n < UVCONST(31858317218647)) {
+ bases[0] = 2;
+ bases[1] = 642735;
+ bases[2] = 553174392;
+ bases[3] = 3046413974;
+ nbases = 4;
+ } else if (n < UVCONST(3071837692357849)) {
+ bases[0] = 2;
+ bases[1] = 75088;
+ bases[2] = 642735;
+ bases[3] = 203659041;
+ bases[4] = 3613982119;
+ nbases = 5;
+ } else {
+ bases[0] = 2;
+ bases[1] = 325;
+ bases[2] = 9375;
+ bases[3] = 28178;
+ bases[4] = 450775;
+ bases[5] = 9780504;
+ bases[6] = 1795265022;
+ nbases = 7;
+ }
+#else
+ /* More standard bases */
+ if (n < UVCONST(9080191)) {
+ bases[0] = 31; bases[1] = 73; nbases = 2;
+ } else if (n < UVCONST(4759123141)) {
+ bases[0] = 2; bases[1] = 7; bases[2] = 61; nbases = 3;
+ } else if (n < UVCONST(21652684502221)) {
+ bases[0] = 2; bases[1] = 1215; bases[2] = 34862; bases[3] = 574237825;
+ nbases = 4;
+ } else if (n < UVCONST(341550071728321)) {
+ bases[0] = 2; bases[1] = 3; bases[2] = 5; bases[3] = 7; bases[4] = 11;
+ bases[5] = 13; bases[6] = 17; nbases = 7;
+ } else if (n < UVCONST(3825123056546413051)) {
+ bases[0] = 2; bases[1] = 3; bases[2] = 5; bases[3] = 7; bases[4] = 11;
+ bases[5] = 13; bases[6] = 17; bases[7] = 19; bases[8] = 23; nbases = 9;
+ } else {
+ bases[0] = 2; bases[1] = 3; bases[2] = 5; bases[3] = 7; bases[4] = 11;
+ bases[5] = 13; bases[6] = 17; bases[7] = 19; bases[8] = 23; bases[9] = 29;
+ bases[10]= 31; bases[11]= 37;
+ nbases = 12;
+ }
+#endif
+#endif
+ prob_prime = miller_rabin(n, bases, nbases);
+ return 2*prob_prime;
}
+
+
/* Knuth volume 2, algorithm C.
* Very fast for small numbers, grows rapidly.
* SQUFOF is better for numbers nearing the 64-bit limit.
@@ -182,7 +337,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds)
}
for (i = 1; i < rounds; i++) {
- Xi = powmod(Xi, 2, a, n);
+ Xi = powmodadd(Xi, 2, a, n);
f = gcd_ui(Xi - Xm, n);
if ( (f != 1) && (f != n) ) {
factors[nfactors++] = f;
@@ -231,9 +386,9 @@ int prho_factor(UV n, UV *factors, UV rounds)
V = 7;
for (i = 1; i < rounds; i++) {
- U = powmod(U, 2, a, n);
- V = powmod(V, 2, a, n);
- V = powmod(V, 2, a, n);
+ U = powmodadd(U, 2, a, n);
+ V = powmodadd(V, 2, a, n);
+ V = powmodadd(V, 2, a, n);
f = gcd_ui( (U > V) ? U-V : V-U, n);
if ( (f != 1) && (f != n) ) {
@@ -270,12 +425,9 @@ int pminus1_factor(UV n, UV *factors, UV rounds)
return nfactors;
b = 13;
-
for (i = 1; i < rounds; i++) {
- b = powmod(b+1, i, 0, n);
- if (b == 0) b = n;
- b--;
- f = gcd_ui(b, n);
+ b = powmod(b, i, n);
+ f = gcd_ui(b-1, n);
if (f == n) {
factors[nfactors++] = n;
return nfactors;
@@ -435,12 +587,12 @@ int squfof_factor(UV n, UV *factors, UV rounds, UV refactor_above)
factors[nfactors++] = p;
factors[nfactors++] = q;
#else
- if (p >= refactor_above)
+ if ( (p >= refactor_above) && (is_prob_prime(p) == 0) )
nfactors += squfof_factor(p, factors+nfactors, rounds, refactor_above);
else
factors[nfactors++] = p;
- if (q >= refactor_above)
+ if ( (q >= refactor_above) && (is_prob_prime(q) == 0) )
nfactors += squfof_factor(q, factors+nfactors, rounds, refactor_above);
else
factors[nfactors++] = q;
diff --git a/factor.h b/factor.h
index dcd4ee3..3659833 100644
--- a/factor.h
+++ b/factor.h
@@ -18,4 +18,38 @@ extern int prho_factor(UV n, UV *factors, UV maxrounds);
extern int pminus1_factor(UV n, UV *factors, UV maxrounds);
+extern int miller_rabin(UV n, const UV *bases, UV nbases);
+extern int is_prob_prime(UV n);
+
+#if 0
+#include "bitarray.h"
+/* Try to make a quick probable prime test */
+static int is_perhaps_prime(UV n)
+{
+ static const UV bases2[2] = {31, 73};
+ static const UV bases3[3] = {2,7,61};
+ if ( (n == 2) || (n == 3) || (n == 5) || (n == 7) )
+ return 2;
+ if ( (n<2) || ((n% 2)==0) || ((n% 3)==0) || ((n% 5)==0) || ((n% 7)==0) )
+ return 0;
+ if (n < 121)
+ return 2;
+ if (n < UVCONST(9080191))
+ return 2*miller_rabin(n, bases2, 2);
+ else
+ return miller_rabin(n, bases3, 3) * ((n < UVCONST(4759123141)) ? 2 : 1);
+}
+
+static int is_perhaps_prime7(UV n)
+{
+ static const UV bases2[2] = {31, 73};
+ static const UV bases3[3] = {2,7,61};
+ /* n must be >= 73 */
+ if (n < UVCONST(9080191))
+ return 2*miller_rabin(n, bases2, 2);
+ else
+ return miller_rabin(n, bases3, 3) * ((n < UVCONST(4759123141)) ? 2 : 1);
+}
+#endif
+
#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index f83f40f..5ee5404 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -13,7 +13,7 @@ BEGIN {
use base qw( Exporter );
our @EXPORT_OK = qw(
prime_precalc prime_free
- is_prime
+ is_prime is_prob_prime miller_rabin
primes
next_prime prev_prime
prime_count prime_count_lower prime_count_upper prime_count_approx
@@ -215,10 +215,11 @@ including Math::Prime::XS, Math::Prime::FastSieve, and Math::Factor::XS.
=head2 is_prime
-Returns true if the number is prime, false if not.
-
print "$n is prime" if is_prime($n);
+Returns 2 if the number is prime, 0 if not. Also note there are
+probabilistic prime testing functions available.
+
=head2 primes
@@ -351,6 +352,36 @@ generate any primes. Uses the Cipolla 1902 approximation with two
polynomials, plus a correction term for small values to reduce the error.
+=head2 miller_rabin
+
+ my $maybe_prime = miller_rabin($n, 2);
+ my $probably_prime = miller_rabin($n, 2, 3, 5, 7, 11, 13, 17);
+
+Takes a positive number as input and one or more bases. The bases must be
+between C<2> and C<n - 2>. Returns 2 is C<n> is definitely prime, 1 if C<n>
+is probably prime, and 0 if C<n> is definitely composite. Since this is
+just the Miller-Rabin test, a value of 2 is only returned for inputs of
+2 and 3, which are shortcut. If 0 is returned, then the number really is a
+composite. If 1 is returned, we aren't sure.
+
+This is usually used in combination with other tests to make either stronger
+tests or deterministic results for numbers less than some verified limit.
+
+=head2 is_prob_prime
+
+ my $prob_prime = is_prob_prime($n);
+ # Returns 0 (composite), 2 (prime), or 1 (probably prime)
+
+Takes a positive number as input and returns back either 0 (composite),
+2 (definitely prime), or 1 (probably prime).
+
+This is done with a tuned set of Miller-Rabin tests such that the result
+should be deterministic for 64-bit input. Either 2, 3, 4, 5, or 7 Miller-Rabin
+tests are performed (no more than 3 for 32-bit input), and the result will
+then always be 0 (composite) or 2 (prime). A later implementation may switch
+to a BPSW test, depending on speed.
+
+
=head1 UTILITY FUNCTIONS
=head2 prime_precalc
diff --git a/util.c b/util.c
index 169f38b..aa5e545 100644
--- a/util.c
+++ b/util.c
@@ -7,6 +7,7 @@
#include "util.h"
#include "sieve.h"
+#include "factor.h"
#include "bitarray.h"
/*
@@ -59,7 +60,7 @@ static int _is_prime7(UV x)
q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 2;
q = x/i; if (q<i) return 1; if (x==(q*i)) return 0; i += 6;
}
- return 1;
+ return 2;
}
@@ -95,8 +96,16 @@ int is_prime(UV n)
if (n <= get_prime_cache(0, &sieve))
return ((sieve[d] & mtab) == 0);
+#if 0
/* Trial division, mod 30 */
return _is_prime7(n);
+#else
+ if (n < UVCONST(100000000)) {
+ return _is_prime7(n);
+ } else {
+ return is_prob_prime(n); /* We know this works for all 64-bit n */
+ }
+#endif
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list