This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.30 in repository libmath-prime-util-perl.
commit 0c66a3807c5c37cc1e67b4f338039d7e21f902c4 Author: Dana Jacobsen <d...@acm.org> Date: Thu Jul 4 15:28:45 2013 -0700 Almost extra strong lucas prps --- Changes | 1 + XS.xs | 28 ++++++++--- factor.c | 132 +++++++++++++++++++++++++++++++++++++------------ factor.h | 2 + lib/Math/Prime/Util.pm | 26 +++++++--- 5 files changed, 145 insertions(+), 44 deletions(-) diff --git a/Changes b/Changes index 2e7b980..f7ad87e 100644 --- a/Changes +++ b/Changes @@ -11,6 +11,7 @@ Revision history for Perl extension Math::Prime::Util. - Added: is_frobenius_underwood_pseudoprime + is_almost_extra_strong_lucas_pseudoprime lucas_sequence pplus1_factor diff --git a/XS.xs b/XS.xs index fe7851a..bd36861 100644 --- a/XS.xs +++ b/XS.xs @@ -466,10 +466,25 @@ _XS_lucas_sequence(IN UV n, IN IV P, IN IV Q, IN UV k) XPUSHs(sv_2mortal(newSVuv( Qk ))); int -_XS_is_lucas_pseudoprime(IN UV n, int strength) - -int -_XS_is_frobenius_underwood_pseudoprime(IN UV n) +_XS_is_lucas_pseudoprime(IN UV n) + ALIAS: + _XS_is_strong_lucas_pseudoprime = 1 + _XS_is_extra_strong_lucas_pseudoprime = 2 + _XS_is_almost_extra_strong_lucas_pseudoprime = 3 + _XS_is_pari_lucas_pseudoprime = 4 + _XS_is_frobenius_underwood_pseudoprime = 5 + CODE: + switch (ix) { + case 0: RETVAL = _XS_is_lucas_pseudoprime(n, 0); break; + case 1: RETVAL = _XS_is_lucas_pseudoprime(n, 1); break; + case 2: RETVAL = _XS_is_lucas_pseudoprime(n, 2); break; + case 3: RETVAL = _XS_is_almost_extra_strong_lucas_pseudoprime(n,1);break; + case 4: RETVAL = _XS_is_almost_extra_strong_lucas_pseudoprime(n,2);break; + case 5: RETVAL = _XS_is_frobenius_underwood_pseudoprime(n); break; + default:RETVAL = 0; break; + } + OUTPUT: + RETVAL int _XS_is_prob_prime(IN UV n) @@ -734,9 +749,8 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0) if (beg <= end) { ctx = start_segment_primes(beg, end, &segment); while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) { - /* I'm getting a memory leak in the MULTICALL and I'm having no luck - * finding out why it is happening. Don't use this for now. */ - if (0 && !CvISXSUB(cv)) { + /* TODO: Make sure this doesn't leak memory in 5.16 and previous */ + if (!CvISXSUB(cv)) { dMULTICALL; I32 gimme = G_VOID; PUSH_MULTICALL(cv); diff --git a/factor.c b/factor.c index 92b3234..f97bc9a 100644 --- a/factor.c +++ b/factor.c @@ -405,12 +405,14 @@ int _XS_is_prob_prime(UV n) int nbases; int prob_prime; - if (n == 2 || n == 3 || n == 5) return 2; - if (n < 2 || !(n%2) || !(n%3) || !(n%5)) return 0; - if (n < 49) /* 7*7 */ return 2; - if (!(n% 7) || !(n%11) || !(n%13) || !(n%17) || !(n%19) || !(n%23)) return 0; - if (!(n%29) || !(n%31) || !(n%37) || !(n%41) || !(n%43) || !(n%47)) return 0; - if (n < 2809) /* 53*53 */ return 2; + if (n == 2 || n == 3 || n == 5 || n == 7) return 2; + if (n < 2 || !(n%2) || !(n%3) || !(n%5) || !(n%7)) return 0; + if (n < 121) /* 11*11 */ return 2; + if (!(n%11) || !(n%13) || !(n%17) || !(n%19) || !(n%23) || !(n%29) || + !(n%31) || !(n%37) || !(n%41) || !(n%43) || !(n%47) || !(n%53) || + !(n%59) || !(n%61) || !(n%67) || !(n%71) || !(n%73) || !(n%79) || + !(n%83) || !(n%89) || !(n%97)) return 0; + if (n < 10201) /* 101*101 */ return 2; #if BITS_PER_WORD == 32 @@ -428,7 +430,7 @@ int _XS_is_prob_prime(UV n) /* Verified with Feitsma database. No counterexamples below 2^64. * This is faster than multiple M-R routines once we're over 32-bit */ if (n >= UVCONST(4294967295)) { - prob_prime = _SPRP2(n) && _XS_is_lucas_pseudoprime(n,2); + prob_prime = _SPRP2(n) && _XS_is_almost_extra_strong_lucas_pseudoprime(n,1); return 2*prob_prime; } @@ -489,7 +491,6 @@ int _XS_is_prob_prime(UV n) void lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, IV P, IV Q, UV k) { UV U, V, b, Dmod, Qmod, Pmod, Qk; - IV D; if (k == 0) { *Uret = 0; @@ -498,8 +499,8 @@ void lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, IV P, IV Q, UV k) return; } - Qmod = (Q < 0) ? Q + n : Q; - Pmod = (P < 0) ? P + n : P; + Qmod = (Q < 0) ? (UV)(Q + (IV)n) : (UV)Q; + Pmod = (P < 0) ? (UV)(P + (IV)n) : (UV)P; Dmod = submod( mulmod(Pmod, Pmod, n), mulmod(4, Qmod, n), n); MPUassert(Dmod != 0, "lucas_seq: D is 0"); U = 1; @@ -650,6 +651,70 @@ int _XS_is_lucas_pseudoprime(UV n, int strength) return 0; } +/* A generalization of Pari's shortcut to the extra-strong Lucas test. + * I've added a gcd check at the top, which needs to be done and also results + * in fewer pseudoprimes. Pari always does trial division to 100 first so + * is unlikely to come up there. This only calculate V, which can be done + * faster, but that means we have more pseudoprimes than the standard + * extra-strong test. + * + * increment: 1 for Baillie OEIS, 2 for Pari. + * + * With increment = 1, these results will be a subset of the extra-strong + * Lucas pseudoprimes. With increment = 2, we produce Pari's results. + */ +int _XS_is_almost_extra_strong_lucas_pseudoprime(UV n, UV increment) +{ + UV P, V, d, s; + + if (n == 2 || n == 3 || n == 5) return 1; + if (n < 7 || (n%2) == 0) return 0; + if (n == UV_MAX) return 0; + + P = 3; + while (1) { + UV D = P*P - 4; + d = gcd_ui(D, n); + if (d > 1 && d < n) + return 0; + if (jacobi_iu(D, n) == -1) + break; + if (P == (3+20*increment) && is_perfect_square(n, 0)) return 0; + P += increment; + } + + d = n+1; + s = 0; + while ( (d & 1) == 0 ) { s++; d >>= 1; } + + { + UV W, b; + V = P; + W = mulsubmod(P, P, 2, n); + { UV v = d; b = 1; while (v >>= 1) b++; } + while (b-- > 1) { + if ( (d >> (b-1)) & UVCONST(1) ) { + V = mulsubmod(V, W, P, n); + W = mulsubmod(W, W, 2, n); + } else { + W = mulsubmod(V, W, P, n); + V = mulsubmod(V, V, 2, n); + } + } + } + + if (V == 2 || V == (n-2)) + return 1; + while (s-- > 1) { + if (V == 0) + return 1; + V = mulsubmod(V, V, 2, n); + if (V == 2) + return 0; + } + return 0; +} + UV _XS_divisor_sum(UV n) { @@ -1348,38 +1413,43 @@ int racing_squfof_factor(UV n, UV *factors, UV rounds) * test from section 9 of "Quadratic Composite Tests" by Paul Underwood. * * Given the script: - * time mpu 'forprimes { Math::Prime::Util::_XS_is_frobenius_underwood_pseudoprime($_); Math::Prime::Util::_XS_is_frobenius_underwood_pseudoprime($_+2); } 100_000_000' + * time mpu 'forprimes { Math::Prime::Util::_XS_is_frobenius_underwood_pseudoprime($_); Math::Prime::Util::_XS_is_frobenius_underwood_pseudoprime($_+2); } 500_000_000' * and replacing the tests appropriately, I get these times: * - * 0.57 $_ (cost of empty loop) - * 6.89 _XS_is_pseudoprime($_,2) - * 6.82 _XS_miller_rabin($_,2) - * 17.42 _XS_is_lucas_pseudoprime($_,0) (standard) - * 17.15 _XS_is_lucas_pseudoprime($_,1) (strong) - * 11.81 _XS_is_lucas_pseudoprime($_,2) (extra strong) - * 13.07 _XS_is_frobenius_underwood_pseudoprime($_) - * 7.87 _XS_is_prob_prime($_) - * 8.74 _XS_is_prime($_) + * 0.87 $_ (cost of empty loop) + * 21.37 _XS_is_pseudoprime($_,2) + * 22.42 _XS_miller_rabin($_,2) + * 44.53 _XS_is_lucas_pseudoprime($_) + * 43.95 _XS_is_strong_lucas_pseudoprime($_) + * 40.09 _XS_is_extra_strong_lucas_pseudoprime($_) + * 25.86 _XS_is_almost_extra_strong_lucas_pseudoprime($_) + * 42.40 _XS_is_frobenius_underwood_pseudoprime($_) + * 27.02 _XS_is_prob_prime($_) + * 27.24 _XS_is_prime($_) * * At these sizes is_prob_prime is doing 1-2 M-R tests. The input validation * is adding a noticeable overhead to is_prime. * - * With a set of 10k 64-bit random primes; 'do { die unless ... } for 1..500' + * With a set of 100k 64-bit random primes; 'do { die unless ... } for 1..50' * - * 0.36 empty loop - * 12.38 _XS_is_pseudoprime($_,2) - * 12.05 _XS_miller_rabin($_,2) - * 24.95 _XS_is_extra_strong_lucas_pseudoprime($_) - * 22.35 _XS_is_frobenius_underwood_pseudoprime($_) - * 36.67 _XS_is_prob_prime($_) - * 37.24 _XS_is_prime($_) + * 0.32 empty loop + * 10.25 _XS_is_pseudoprime($_,2) + * 10.06 _XS_miller_rabin($_,2) + * 22.02 _XS_is_lucas_pseudoprime($_) + * 21.81 _XS_is_strong_lucas_pseudoprime($_) + * 20.99 _XS_is_extra_strong_lucas_pseudoprime($_) + * 14.01 _XS_is_almost_extra_strong_lucas_pseudoprime($_) + * 18.44 _XS_is_frobenius_underwood_pseudoprime($_) + * 24.11 _XS_is_prob_prime($_) + * 24.06 _XS_is_prime($_) * * At this point is_prob_prime has transitioned to BPSW. * * Calling a powmod a 'Selfridge' unit, then we see: - * 1 Selfridge unit M-R test - * 2 Selfridge units Lucas or Frobenius-Underwood - * 3 Selfridge units BPSW + * 1 Selfridge unit M-R test + * 1.4 Selfridge unit "almost extra strong" Lucas + * 2 Selfridge units Lucas or Frobenius-Underwood + * 3 Selfridge units BPSW (standard, strong, or extra-strong) * * We try to structure the primality test like: * 1) simple divisibility very fast primes and ~10% of composites diff --git a/factor.h b/factor.h index b913345..c7dd3ae 100644 --- a/factor.h +++ b/factor.h @@ -14,6 +14,7 @@ extern int holf_factor(UV n, UV *factors, UV rounds); extern int pbrent_factor(UV n, UV *factors, UV maxrounds, UV a); extern int prho_factor(UV n, UV *factors, UV maxrounds); extern int pminus1_factor(UV n, UV *factors, UV B1, UV B2); +extern int pplus1_factor(UV n, UV *factors, UV B); extern int squfof_factor(UV n, UV *factors, UV rounds); extern int racing_squfof_factor(UV n, UV *factors, UV rounds); @@ -25,6 +26,7 @@ extern int _XS_is_prob_prime(UV n); extern void lucas_seq(UV* U, UV* V, UV* Qk, UV n, IV P, IV Q, UV k); extern int _XS_is_lucas_pseudoprime(UV n, int strength); extern int _XS_is_frobenius_underwood_pseudoprime(UV n); +extern int _XS_is_almost_extra_strong_lucas_pseudoprime(UV n, UV increment); extern UV _XS_divisor_sum(UV n); diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index a06f5b1..182f6f8 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -17,8 +17,11 @@ our @EXPORT_OK = prime_precalc prime_memfree is_prime is_prob_prime is_provable_prime is_provable_prime_with_cert prime_certificate verify_prime - is_pseudoprime is_strong_pseudoprime is_lucas_pseudoprime - is_strong_lucas_pseudoprime is_extra_strong_lucas_pseudoprime + is_pseudoprime is_strong_pseudoprime + is_lucas_pseudoprime + is_strong_lucas_pseudoprime + is_extra_strong_lucas_pseudoprime + is_almost_extra_strong_lucas_pseudoprime is_frobenius_underwood_pseudoprime is_aks_prime miller_rabin @@ -1601,7 +1604,7 @@ sub is_strong_pseudoprime { sub is_lucas_pseudoprime { my($n) = shift; _validate_num($n) || _validate_positive_integer($n); - return _XS_is_lucas_pseudoprime($n, 0) + return _XS_is_lucas_pseudoprime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL; return Math::Prime::Util::GMP::is_lucas_pseudoprime("$n") if $_HAVE_GMP && defined &Math::Prime::Util::GMP::is_lucas_pseudoprime; @@ -1611,7 +1614,7 @@ sub is_lucas_pseudoprime { sub is_strong_lucas_pseudoprime { my($n) = shift; _validate_num($n) || _validate_positive_integer($n); - return _XS_is_lucas_pseudoprime($n, 1) + return _XS_is_strong_lucas_pseudoprime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL; return Math::Prime::Util::GMP::is_strong_lucas_pseudoprime("$n") if $_HAVE_GMP; @@ -1621,7 +1624,7 @@ sub is_strong_lucas_pseudoprime { sub is_extra_strong_lucas_pseudoprime { my($n) = shift; _validate_num($n) || _validate_positive_integer($n); - return _XS_is_lucas_pseudoprime($n, 2) + return _XS_is_extra_strong_lucas_pseudoprime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL; return Math::Prime::Util::GMP::is_extra_strong_lucas_pseudoprime("$n") if $_HAVE_GMP @@ -1629,6 +1632,17 @@ sub is_extra_strong_lucas_pseudoprime { return Math::Prime::Util::PP::is_extra_strong_lucas_pseudoprime($n); } +sub is_almost_extra_strong_lucas_pseudoprime { + my($n) = shift; + _validate_num($n) || _validate_positive_integer($n); + return _XS_is_almost_extra_strong_lucas_pseudoprime($n) + if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL; + return Math::Prime::Util::GMP::is_almost_extra_strong_lucas_pseudoprime("$n") + if $_HAVE_GMP + && defined &Math::Prime::Util::GMP::is_almost_extra_strong_lucas_pseudoprime; + return Math::Prime::Util::PP::is_almost_extra_strong_lucas_pseudoprime($n); +} + sub is_frobenius_underwood_pseudoprime { my($n) = shift; _validate_num($n) || _validate_positive_integer($n); @@ -1806,7 +1820,7 @@ sub verify_prime { my $intn = int($n->bstr); if ($n->bcmp("$intn") == 0 && $intn <= $_XS_MAXVAL) { $bpsw = _XS_miller_rabin($intn, 2) - && _XS_is_lucas_pseudoprime($intn,2); + && _XS_is_extra_strong_lucas_pseudoprime($intn); } elsif ($_HAVE_GMP) { $bpsw = Math::Prime::Util::GMP::is_prob_prime($n); } else { -- 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