[libmath-prime-util-perl] 02/40: Add Lucas sequence and have all 3 Lucas tests use it
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:01 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list