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 aa85624de37b697497e8e2106fd6bfd3944325c4 Author: Dana Jacobsen <d...@acm.org> Date: Tue Jun 18 11:42:48 2013 -0700 Add Lucas sequence and have all 3 Lucas tests use it --- Changes | 3 + XS.xs | 2 +- factor.c | 170 +++++++++++++++++++++++++++++++++++-------------- factor.h | 2 +- lib/Math/Prime/Util.pm | 8 ++- 5 files changed, 132 insertions(+), 53 deletions(-) diff --git a/Changes b/Changes index a798d32..75217f4 100644 --- a/Changes +++ b/Changes @@ -1,5 +1,8 @@ Revision history for Perl extension Math::Prime::Util. +0.30 + - Add standard and strong Lucas tests in XS. + 0.29 30 May 2013 - Fix a signed vs. unsigned char issue in ranged moebius. Thanks to the diff --git a/XS.xs b/XS.xs index 8f32f45..dc8a882 100644 --- a/XS.xs +++ b/XS.xs @@ -448,7 +448,7 @@ _XS_miller_rabin(IN UV n, ...) RETVAL int -_XS_is_extra_strong_lucas_pseudoprime(IN UV n) +_XS_is_lucas_pseudoprime(IN UV n, int strength) int _XS_is_frobenius_underwood_pseudoprime(IN UV n) diff --git a/factor.c b/factor.c index 1dfcc4d..5afa18e 100644 --- a/factor.c +++ b/factor.c @@ -428,7 +428,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_extra_strong_lucas_pseudoprime(n); + prob_prime = _SPRP2(n) && _XS_is_lucas_pseudoprime(n,2); return 2*prob_prime; } @@ -485,72 +485,142 @@ int _XS_is_prob_prime(UV n) #endif } -/* Extra Strong Lucas test. +/* Generic Lucas sequence for any appropriate P and Q */ +static 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; + *Vret = 2; + *Qkret = Q; + return; + } + + Qmod = (Q < 0) ? Q + n : Q; + Pmod = (P < 0) ? P + n : P; + Dmod = addmod( mulmod(Pmod, Pmod, n), n - mulmod(4, Qmod, n), n ); + MPUassert(Dmod != 0, "lucas_seq: D is 0"); + U = 1; + V = Pmod; + Qk = Qmod; + { UV v = k; b = 1; while (v >>= 1) b++; } + + if (Q == 1) { + while (b > 1) { + U = mulmod(U, V, n); + V = muladdmod(V, V, n-2, n); + b--; + if ( (k >> (b-1)) & UVCONST(1) ) { + UV t2 = mulmod(U, Dmod, n); + U = muladdmod(U, Pmod, V, n); + if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; } + V = muladdmod(V, P, t2, n); + if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; } + } + } + } else { + while (b > 1) { + U = mulmod(U, V, n); + V = muladdmod(V, V, n - addmod(Qk, Qk, n), n); + Qk = sqrmod(Qk, n); + b--; + if ( (k >> (b-1)) & UVCONST(1) ) { + UV t2 = mulmod(U, Dmod, n); + U = muladdmod(U, Pmod, V, n); + if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; } + V = muladdmod(V, P, t2, n); + if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; } + Qk = mulmod(Qk, Qmod, n); + } + } + } + *Uret = U; + *Vret = V; + *Qkret = Qk; +} + +/* Lucas tests: + * 0: Standard + * 1: Strong + * 2: Extra Strong (Mo/Jones/Grantham) * * Goal: * (1) no false results when combined with the SPRP-2 test. * (2) fast enough to use SPRP-2 + this in place of 3+ M-R tests. * - * Why the extra strong test? If we use Q=1, the code is both simpler and - * faster. But the Selfridge parameters have P==1 Q!=1. Once we decide to - * go with the Q==1, P!=1 method, then we may as well use the extra strong - * test so we can verify results (e.g. OEIS A217719). There is no cost to - * run this vs. the strong or standard test. + * For internal purposes, we typically want to use the extra strong test + * because it is slightly faster (Q = 1 simplies things). None of them have + * any false positives for the BPSW test. * - * This runs about 7x faster than the GMP strong test, and about 2x slower - * than our M-R tests. + * This runs 4-7x faster than MPU::GMP, which means ~10x faster than most GMP + * implementations. It is about 2x slower than a single M-R test. */ -int _XS_is_extra_strong_lucas_pseudoprime(UV n) +int _XS_is_lucas_pseudoprime(UV n, int strength) { - UV P, D, Q, U, V, d, s, b; - int const _verbose = _XS_get_verbose(); + IV P, Q, D; + UV U, V, Qk, d, s, b; if (n == 2 || n == 3) return 1; if (n < 5 || (n%2) == 0) return 0; if (n == UV_MAX) return 0; - P = 3; - Q = 1; - while (1) { - D = P*P - 4; - if (gcd_ui(D, n) > 1 && gcd_ui(D, n) != n) return 0; - if (jacobi_iu(D, n) == -1) - break; - /* Perhaps n is a perfect square? */ - if (P == 21 && is_perfect_square(n, 0)) return 0; - P++; + if (strength < 2) { + UV Du = 5; + IV sign = 1; + while (1) { + D = Du * sign; + if (gcd_ui(Du, n) > 1 && gcd_ui(Du, n) != n) return 0; + if (jacobi_iu(D, n) == -1) + break; + if (Du == 21 && is_perfect_square(n, 0)) return 0; + Du += 2; + sign = -sign; + } + P = 1; + Q = (1 - D) / 4; + } else { + P = 3; + Q = 1; + while (1) { + D = P*P - 4; + if (gcd_ui(D, n) > 1 && gcd_ui(D, n) != n) return 0; + if (jacobi_iu(D, n) == -1) + break; + if (P == 21 && is_perfect_square(n, 0)) return 0; + P++; + } } - if (_verbose>3) printf("N: %lu D: %ld P: %lu Q: %ld\n", n, D, P, Q); - MPUassert( D == ((IV)(P*P)) - 4*Q , "incorrect DPQ"); + MPUassert( D == (P*P - 4*Q) , "is_lucas_pseudoprime: incorrect DPQ"); - U = 1; - V = P; d = n+1; s = 0; - while ( (d & 1) == 0 ) { s++; d >>= 1; } - { UV v = d; b = 1; while (v >>= 1) b++; } - - if (_verbose>3) printf("U=%lu V=%lu\n", U, V); - while (b > 1) { - U = mulmod(U, V, n); - V = muladdmod(V, V, n-2, n); - b--; - if (_verbose>3) printf("U2k=%lu V2k=%lu\n", U, V); - if ( (d >> (b-1)) & UVCONST(1) ) { - UV t2 = mulmod(U, D, n); - U = muladdmod(U, P, V, n); - if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; } - V = muladdmod(V, P, t2, n); - if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; } + if (strength > 0) + while ( (d & 1) == 0 ) { s++; d >>= 1; } + lucas_seq(&U, &V, &Qk, n, P, Q, d); + + if (strength == 0) { + if (U == 0) + return 1; + } else if (strength == 1) { + if (U == 0 || V == 0) + return 1; + while (s--) { + V = muladdmod(V, V, n - addmod(Qk, Qk, n), n); + if (V == 0) + return 1; + if (s) + Qk = sqrmod(Qk, n); } - if (_verbose>3) printf("U=%lu V=%lu\n", U, V); - } - if ( (U == 0 && (V == 2 || V == (n-2))) || (V == 0) ) - return 1; - while (s--) { - V = muladdmod(V, V, n-2, n); - if (V == 0) + } else { + if ( (U == 0 && (V == 2 || V == (n-2))) || (V == 0) ) return 1; + while (s--) { + V = muladdmod(V, V, n-2, n); + if (V == 0) + return 1; + } } return 0; } @@ -1208,7 +1278,9 @@ int racing_squfof_factor(UV n, UV *factors, UV rounds) * 0.57 $_ (cost of empty loop) * 6.89 _XS_is_pseudoprime($_,2) * 6.82 _XS_miller_rabin($_,2) - * 11.81 _XS_is_extra_strong_lucas_pseudoprime($_) + * 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($_) diff --git a/factor.h b/factor.h index d017d49..bce382c 100644 --- a/factor.h +++ b/factor.h @@ -22,7 +22,7 @@ extern int _XS_is_pseudoprime(UV n, UV a); extern int _XS_miller_rabin(UV n, const UV *bases, int nbases); extern int _SPRP2(UV n); extern int _XS_is_prob_prime(UV n); -extern int _XS_is_extra_strong_lucas_pseudoprime(UV n); +extern int _XS_is_lucas_pseudoprime(UV n, int strength); extern int _XS_is_frobenius_underwood_pseudoprime(UV n); extern UV _XS_divisor_sum(UV n); diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index d3c1c98..129dcfe 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -1599,6 +1599,8 @@ 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) + 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; return Math::Prime::Util::PP::is_lucas_pseudoprime($n); @@ -1607,6 +1609,8 @@ 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) + if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL; return Math::Prime::Util::GMP::is_strong_lucas_pseudoprime("$n") if $_HAVE_GMP; return Math::Prime::Util::PP::is_strong_lucas_pseudoprime($n); @@ -1615,7 +1619,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_extra_strong_lucas_pseudoprime($n) + return _XS_is_lucas_pseudoprime($n, 2) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL; return Math::Prime::Util::GMP::is_extra_strong_lucas_pseudoprime("$n") if $_HAVE_GMP @@ -1771,7 +1775,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_extra_strong_lucas_pseudoprime($intn); + && _XS_is_lucas_pseudoprime($intn,2); } 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