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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits