[libmath-prime-util-perl] 07/15: Add two more probable prime tests
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:48:46 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.29
in repository libmath-prime-util-perl.
commit 61ad7e51e631c1818a41cec3a5566f6ccd4ec3c1
Author: Dana Jacobsen <dana at acm.org>
Date: Wed May 29 15:30:22 2013 -0700
Add two more probable prime tests
---
Changes | 13 +++-
TODO | 6 +-
XS.xs | 2 +-
factor.c | 31 ++++++---
lib/Math/Prime/Util.pm | 83 ++++++++++++++++++-----
lib/Math/Prime/Util/PP.pm | 167 +++++++++++++++++++++++++++++++++++-----------
t/17-pseudoprime.t | 16 ++++-
7 files changed, 245 insertions(+), 73 deletions(-)
diff --git a/Changes b/Changes
index 38bbf9f..1f7eec3 100644
--- a/Changes
+++ b/Changes
@@ -2,13 +2,20 @@ Revision history for Perl extension Math::Prime::Util.
0.29 xx May 2013
- - Fix a signed vs. unsigned char issue in ranged moebius.
+ - Fix a signed vs. unsigned char issue in ranged moebius. Thanks to the
+ Debian testers for finding this.
+
+ - XS is_prob_prime / is_prime now use a BPSW-style test (SPRP2 plus
+ extra strong Lucas test) for values over 2^32. This results in up
+ to 2.5x faster performance for large 64-bit values on most machines.
+ All PSP2s have been verified with Jan Feitsma's database.
- forprimes now uses a segmented sieve. This (1) allows arbitrary 64-bit
ranges with good memory use, and (2) allows nesting on threaded perls.
- - Added XS strong Lucas test, but not using Selfridge parameters.
- Hence we only use it internally.
+ - Added is_pseudoprime (for Fermat probable prime test) and
+ is_extra_strong_lucas_pseudoprime (for the Mo/Jones/Grantham extra
+ strong Lucas test).
0.28 23 May 2013
diff --git a/TODO b/TODO
index aaa6311..9be763e 100644
--- a/TODO
+++ b/TODO
@@ -41,9 +41,11 @@
- LMO prime count
- QS factoring
-- Make forprimes use a segment sieve
-
- Figure out a way to make the internal FOR_EACH_PRIME macros use a segmented
sieve.
- prime_count_approx on a 100-digit bigint is really slow without MPFR.
+
+- Add is_lucas_pseudoprime to complete our set
+
+- re-verify with Feitsma database using extra strong test
diff --git a/XS.xs b/XS.xs
index 2d23cd5..9342e4a 100644
--- a/XS.xs
+++ b/XS.xs
@@ -429,7 +429,7 @@ _XS_miller_rabin(IN UV n, ...)
RETVAL
int
-_XS_is_strong_lucas_pseudoprime(IN UV n)
+_XS_is_extra_strong_lucas_pseudoprime(IN UV n)
int
_XS_is_frobenius_underwood_pseudoprime(IN UV n)
diff --git a/factor.c b/factor.c
index db8bfba..fce1268 100644
--- a/factor.c
+++ b/factor.c
@@ -314,7 +314,8 @@ int _XS_is_pseudoprime(UV n, UV a)
UV x;
UV const nm1 = n-1;
- if (n <= 3) return 0;
+ if (n == 2 || n == 3) return 1;
+ if (n < 5) return 0;
if (a < 2) croak("Base %"UVuf" is invalid", a);
if (a >= n) {
a %= n;
@@ -427,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_strong_lucas_pseudoprime(n);
+ prob_prime = _SPRP2(n) && _XS_is_extra_strong_lucas_pseudoprime(n);
return 2*prob_prime;
}
@@ -484,12 +485,22 @@ int _XS_is_prob_prime(UV n)
#endif
}
-/* Strong Lucas with non-Selfridge params. Basically the same speed as
- * the standard test above. The parameter choice leads to almost 2x
- * as many pseudoprimes as the Selfridge params. This runs about 7x
- * faster than the GMP version, and about 2x slower than our M-R tests.
- * The goal: no false results when combined with the SPRP-2 test. */
-int _XS_is_strong_lucas_pseudoprime(UV n)
+/* Extra Strong Lucas test.
+ *
+ * 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.
+ *
+ * This runs about 7x faster than the GMP strong test, and about 2x slower
+ * than our M-R tests.
+ */
+int _XS_is_extra_strong_lucas_pseudoprime(UV n)
{
UV P, D, Q, U, V, t, t2, d, s, b;
int const _verbose = _XS_get_verbose();
@@ -507,7 +518,7 @@ int _XS_is_strong_lucas_pseudoprime(UV n)
break;
/* Perhaps n is a perfect square? */
if (P == 21 && is_perfect_square(n, 0)) return 0;
- P += 2;
+ 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");
@@ -534,7 +545,7 @@ int _XS_is_strong_lucas_pseudoprime(UV n)
}
if (_verbose>3) printf("U=%lu V=%lu\n", U, V);
}
- if (U == 0 || V == 0)
+ if ( (U == 0 && (V == 2 || V == (n-2))) || (V == 0) )
return 1;
while (s--) {
V = muladdmod(V, V, n-2, n);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index ce7f1c5..134e924 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -17,7 +17,8 @@ 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_strong_pseudoprime is_strong_lucas_pseudoprime
+ is_pseudoprime is_strong_pseudoprime
+ is_strong_lucas_pseudoprime is_extra_strong_lucas_pseudoprime
is_aks_prime
miller_rabin
primes
@@ -53,6 +54,7 @@ sub _import_nobigint {
undef *factor; *factor = \&_XS_factor;
#undef *prime_count; *prime_count = \&_XS_prime_count;
undef *nth_prime; *nth_prime = \&_XS_nth_prime;
+ undef *is_pseudoprime; *is_pseudoprime = \&_XS_is_pseudoprime;
undef *is_strong_pseudoprime; *is_strong_pseudoprime = \&_XS_miller_rabin;
undef *miller_rabin; *miller_rabin = \&_XS_miller_rabin;
undef *moebius; *moebius = \&_XS_moebius;
@@ -169,7 +171,7 @@ our $_Neg_Infinity = -$_Infinity;
# make $n into a bigint. This is debatable, but they *did* hand us a
# string with a big integer in it. The big gotcha here is that
# is_strong_lucas_pseudoprime does bigint computations, so it will load
-# up bigint and there is no way to unload it.
+# up Math::BigInt and there is no way to unload it.
#
# 3) if (ref($n) =~ /^Math::Big/)
# $n is a big int, float, or rat. We probably want this as an int.
@@ -1585,6 +1587,17 @@ sub factor {
return Math::Prime::Util::PP::factor($n);
}
+sub is_pseudoprime {
+ my($n, $a) = @_;
+ _validate_num($n) || _validate_positive_integer($n);
+ _validate_num($a) || _validate_positive_integer($a);
+ return _XS_is_pseudoprime($n, $a)
+ if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
+ return Math::Prime::Util::GMP::is_pseudoprime($n, $a)
+ if $_HAVE_GMP && defined &Math::Prime::Util::GMP::is_pseudoprime;
+ return Math::Prime::Util::PP::is_pseudoprime($n, $a);
+}
+
sub is_strong_pseudoprime {
my($n) = shift;
_validate_num($n) || _validate_positive_integer($n);
@@ -1594,13 +1607,27 @@ sub is_strong_pseudoprime {
return Math::Prime::Util::PP::miller_rabin($n, @_);
}
+# TODO: standard lucas pseudoprime
+
sub is_strong_lucas_pseudoprime {
my($n) = shift;
_validate_num($n) || _validate_positive_integer($n);
- return Math::Prime::Util::GMP::is_strong_lucas_pseudoprime("$n") if $_HAVE_GMP;
+ return Math::Prime::Util::GMP::is_strong_lucas_pseudoprime("$n")
+ if $_HAVE_GMP;
return Math::Prime::Util::PP::is_strong_lucas_pseudoprime($n);
}
+sub is_extra_strong_lucas_pseudoprime {
+ my($n) = shift;
+ _validate_num($n) || _validate_positive_integer($n);
+ 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
+ && defined &Math::Prime::Util::GMP::is_extra_strong_lucas_pseudoprime;
+ return Math::Prime::Util::PP::is_extra_strong_lucas_pseudoprime($n);
+}
+
sub miller_rabin {
#warn "miller_rabin() is deprecated. Use is_strong_pseudoprime instead.";
return is_strong_pseudoprime(@_);
@@ -1746,7 +1773,10 @@ sub verify_prime {
return 0;
}
my $bpsw = 0;
- if ($_HAVE_GMP) {
+ if ($n <= $_XS_MAXVAL) {
+ $bpsw = _XS_miller_rabin($n, 2)
+ && _XS_is_extra_strong_lucas_pseudoprime($n);
+ } elsif ($_HAVE_GMP) {
$bpsw = Math::Prime::Util::GMP::is_prob_prime($n);
} else {
$bpsw = Math::Prime::Util::PP::miller_rabin($n, 2)
@@ -2411,8 +2441,9 @@ Version 0.28
# Strong pseudoprime test with multiple bases, using Miller-Rabin
say "$n is a prime or 2/7/61-psp" if is_strong_pseudoprime($n, 2, 7, 61);
- # Strong Lucas-Selfridge test
+ # Strong Lucas-Selfridge and extra strong Lucas tests
say "$n is a prime or slpsp" if is_strong_lucas_pseudoprime($n);
+ say "$n is a prime or eslpsp" if is_extra_strong_lucas_pseudoprime($n);
# step to the next prime (returns 0 if not using bigints and we'd overflow)
$n = next_prime($n);
@@ -2556,6 +2587,7 @@ it will make it run much faster.
Some of the functions, including:
factor
+ is_pseudoprime
is_strong_pseudoprime
nth_prime
moebius
@@ -2848,14 +2880,21 @@ generate any primes. Uses the Cipolla 1902 approximation with two
polynomials, plus a correction for small values to reduce the error.
+=head2 is_pseudoprime
+
+Takes a positive number C<n> and a base C<a> as input, and returns 1 if
+C<n> is a probable prime to base C<a>. This is the simple Fermat primality
+test. Removing primes, given base 2 this produces the sequence
+L<OEIS A001567|http://oeis.org/A001567>.
+
=head2 is_strong_pseudoprime
my $maybe_prime = is_strong_pseudoprime($n, 2);
my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17);
Takes a positive number as input and one or more bases. The bases must be
-greater than C<1>. Returns 1 if the input is a prime or a strong
-pseudoprime to all of the bases, and 0 if not.
+greater than C<1>. Returns 1 if the input is a strong probable prime
+to all of the bases, and 0 if not.
If 0 is returned, then the number really is a composite. If 1 is returned,
then it is either a prime or a strong pseudoprime to all the given bases.
@@ -2869,7 +2908,7 @@ selected bases are required to give correct primality test results for any
32-bit number). Given the small chances of passing multiple bases, there
are some math packages that just use multiple MR tests for primality testing.
-Even numbers other than 2 will always return 0 (composite). While the
+Even inputs other than 2 will always return 0 (composite). While the
algorithm does run with even input, most sources define it only on odd input.
Returning composite for all non-2 even input makes the function match most
other implementations including L<Math::Primality>'s C<is_strong_pseudoprime>
@@ -2883,10 +2922,21 @@ An alias for C<is_strong_pseudoprime>. This name is being deprecated.
=head2 is_strong_lucas_pseudoprime
Takes a positive number as input, and returns 1 if the input is a strong
-Lucas pseudoprime using the Selfridge method of choosing D, P, and Q (some
+Lucas probable prime using the Selfridge method of choosing D, P, and Q (some
sources call this a strong Lucas-Selfridge pseudoprime). This is one half
of the BPSW primality test (the Miller-Rabin strong pseudoprime test with
-base 2 being the other half).
+base 2 being the other half). Removing primes, this produces the sequence
+L<OEIS A217255|http://oeis.org/A217255>.
+
+=head2 is_extra_strong_lucas_pseudoprime
+
+Takes a positive number as input, and returns 1 if the input passes the extra
+strong Lucas test (as defined in
+L<Grantham 2000|http://www.ams.org/mathscinet-getitem?mr=1680879>).
+This has slightly more restrictive conditions than the strong Lucas test,
+but uses different starting parameters so is not directly comparable.
+Removing primes, this produces the sequence
+L<OEIS A217719|http://oeis.org/A217719>.
=head2 is_prob_prime
@@ -2897,11 +2947,12 @@ base 2 being the other half).
Takes a positive number as input and returns back either 0 (composite),
2 (definitely prime), or 1 (probably prime).
-For 64-bit input (native or bignum), this uses a tuned set of Miller-Rabin
-tests such that the result will be deterministic. 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 change the internals, but the results will be identical.
+For 64-bit input (native or bignum), this uses either a deterministic set of
+Miller-Rabin tests (1, 2, or 3 tests) or a strong BPSW test consisting of a
+single base-2 strong probable prime test followed by a strong Lucas test.
+This has been verified with Jan Feitsma's 2-PSP database to produce no false
+results for 64-bit inputs. Hence the result will always be 0 (composite) or
+2 (prime).
For inputs larger than C<2^64>, a strong Baillie-PSW primality test is
performed (also called BPSW or BSW). This is a probabilistic test, so only
@@ -2926,7 +2977,7 @@ and do not return false positives.
Using the L<Math::Prime::Util::GMP> module is B<highly recommended> for doing
primality proofs, as it is much, much faster. The pure Perl code is just not
-very fast for this type of operation, nor does it have the best algorithms.
+fast for this type of operation, nor does it have the best algorithms.
It should suffice fine for proofs of up to 40 digit primes, while the latest
MPU::GMP works for primes of hundreds of digits.
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index e413d0b..5708dd5 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -748,6 +748,24 @@ sub _log2 {
+sub is_pseudoprime {
+ my($n, $base) = @_;
+ _validate_positive_integer($n);
+ _validate_positive_integer($base);
+
+ return 1 if $n == 2 || $n == 3;
+ return 0 if $n < 5;
+ croak "Base $base is invalid" if $base < 2;
+ if ($base >= $n) {
+ $base = $base % $n;
+ return 1 if $base <= 1 || $base == $n-1;
+ }
+ my $x = (ref($n) eq 'Math::BigInt')
+ ? $n->copy->bzero->badd($base)->bmodpow($n-1,$n)
+ : _native_powmod($base, $n-1, $n);
+ return ($x == 1);
+}
+
sub miller_rabin {
my($n, @bases) = @_;
_validate_positive_integer($n);
@@ -821,6 +839,26 @@ sub miller_rabin {
1;
}
+sub _is_perfect_square {
+ my($n) = @_;
+
+ if (ref($n) eq 'Math::BigInt') {
+ my $mc = int(($n & 31)->bstr);
+ if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
+ my $sq = $n->copy->bsqrt->bfloor;
+ $sq->bmul($sq);
+ return 1 if $sq == $n;
+ }
+ } else {
+ my $mc = $n & 31;
+ if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
+ my $sq = int(sqrt($n));
+ return 1 if ($sq*$sq) == $n;
+ }
+ }
+ 0;
+}
+
# Calculate Jacobi symbol (M|N)
sub _jacobi {
my($n, $m) = @_;
@@ -880,37 +918,9 @@ sub is_strong_lucas_pseudoprime {
my($n) = @_;
_validate_positive_integer($n);
- # We're trying to limit the bignum calculations as much as possible.
- # It's also important to try to match whatever they passed in. For instance
- # if they use a GMP or Pari object, we must do the same. Hence instead of:
- # my $U = Math::BigInt->bone;
- # we do
- # my $U = $n->copy->bone;
- # so U is the same class as n. If they passed in a string or a small value,
- # then we just make it up.
-
return 1 if $n == 2;
return 0 if $n < 2 || ($n % 2) == 0;
-
- # References:
- # http://www.trnicely.net/misc/bpsw.html
- # Math::Primality
-
- # Check for perfect square
- if (ref($n) eq 'Math::BigInt') {
- my $mc = int(($n & 31)->bstr);
- if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
- my $sq = $n->copy->bsqrt->bfloor;
- $sq->bmul($sq);
- return 0 if $sq == $n;
- }
- } else {
- my $mc = $n & 31;
- if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
- my $sq = int(sqrt($n));
- return 0 if ($sq*$sq) == $n;
- }
- }
+ return 0 if _is_perfect_square($n);
# Determine Selfridge D, P, and Q parameters
my $D = _find_jacobi_d_sequence($n);
@@ -944,7 +954,7 @@ sub is_strong_lucas_pseudoprime {
my $V = $ZERO + $P;
my $Qk = $ZERO + $Q;
my $bpos = 0;
- while ($bpos++ < length($dstr)) {
+ while (++$bpos < length($dstr)) {
$U = ($U * $V) % $n;
$V = ($V * $V - 2*$Qk) % $n;
$Qk = ($Qk * $Qk) % $n;
@@ -974,6 +984,71 @@ sub is_strong_lucas_pseudoprime {
return 0;
}
+sub is_extra_strong_lucas_pseudoprime {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ return 1 if $n == 2;
+ return 0 if $n < 2 || ($n % 2) == 0;
+ return 0 if _is_perfect_square($n);
+
+ my ($P, $Q, $D) = (3, 1, 5);
+ while (1) {
+ my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($D, $n)
+ : _gcd_ui($D, $n);
+ return 0 if $gcd > 1 && $gcd != $n; # Found divisor $d
+ last if _jacobi($D, $n) == -1;
+ $P++;
+ $D = $P*$P - 4;
+ }
+ die "Lucas incorrect DPQ: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);
+
+ if (ref($n) ne 'Math::BigInt') {
+ if (!defined $Math::BigInt::VERSION) {
+ eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
+ or do { croak "Cannot load Math::BigInt "; }
+ }
+ $n = Math::BigInt->new("$n");
+ }
+
+ my $m = $n->copy->badd(1);
+ # Traditional d,s:
+ # my $d=$m->copy; my $s=0; while ($d->is_even) { $s++; $d->brsft(1); }
+ # die "Invalid $m, $d, $s\n" unless $m == $d * 2**$s;
+ my $dstr = substr($m->as_bin, 2);
+ $dstr =~ s/(0*)$//;
+ my $s = length($1);
+
+ my $ZERO = $n->copy->bzero;
+ my $U = $ZERO + 1;
+ my $V = $ZERO + $P;
+ my $bpos = 0;
+ while (++$bpos < length($dstr)) {
+ $U->bmul($V)->bmod($n);
+ $V->bmul($V)->bsub(2)->bmod($n);
+ if (substr($dstr,$bpos,1)) {
+ my $U_times_D = $U->copy->bmul($D);
+ $U->bmul($P)->badd($V);
+ $U->badd($n) if $U->is_odd;
+ $U->brsft(1);
+
+ $V->bmul($P)->badd($U_times_D);
+ $V->badd($n) if $V->is_odd;
+ $V->brsft(1);
+ }
+ }
+ $U->bmod($n);
+ $V->bmod($n);
+ return 1 if $U->is_zero && ($V == 2 || $V == ($n-2));
+ return 1 if $V->is_zero;
+
+ foreach my $r (1 .. $s-1) {
+ $V->bmul($V)->bsub(2)->bmod($n);
+ return 1 if $V->is_zero;
+ }
+ return 0;
+}
+
my $_poly_bignum;
sub _poly_new {
@@ -1485,7 +1560,7 @@ sub pminus1_factor {
while (1) {
$pc_end = $B1 if $pc_end > $B1;
@bprimes = @{ primes($pc_beg, $pc_end) };
- foreach my $q (@bprimes) {
+ foreach my $q (@bprimes) {
my($k, $kmin) = ($q, int($B1 / $q));
while ($k <= $kmin) { $k *= $q; }
$t *= $k; # accumulate powers for a
@@ -2653,6 +2728,13 @@ methods take quite a bit of time and space. Think abut whether a bound or
approximation would be acceptable instead.
+=head2 is_pseudoprime
+
+Takes a positive number C<n> and a base C<a> as input, and returns 1 if
+C<n> is a probable prime to base C<a>. This is the simple Fermat primality
+test. Removing primes, given base 2 this produces the sequence
+L<OEIS A001567|http://oeis.org/A001567>.
+
=head2 is_strong_pseudoprime
=head2 miller_rabin
@@ -2661,16 +2743,17 @@ approximation would be acceptable instead.
my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17);
Takes a positive number as input and one or more bases. The bases must be
-greater than 1. 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, there is a good chance the number is prime
-(depending on the input and the bases), but we cannot be sure.
+greater than C<1>. Returns 1 if the input is a strong probable prime
+to all of the bases, and 0 if not.
+
+If 0 is returned, then the number really is a composite. If 1 is returned,
+then it is either a prime or a strong pseudoprime to all the given bases.
+Given enough distinct bases, the chances become very, very strong that the
+number is actually prime.
This is usually used in combination with other tests to make either stronger
tests (e.g. the strong BPSW test) or deterministic results for numbers less
-than some verified limit (such as the C<is_prob_prime> function in this module).
+than some verified limit.
=head2 is_strong_lucas_pseudoprime
@@ -2681,6 +2764,14 @@ sources call this a strong Lucas-Selfridge pseudoprime). This is one half
of the BPSW primality test (the Miller-Rabin strong pseudoprime test with
base 2 being the other half).
+=head2 is_extra_strong_lucas_pseudoprime
+
+Takes a positive number as input, and returns 1 if the input passes the extra
+strong Lucas test (as defined in
+L<Grantham 2000|http://www.ams.org/mathscinet-getitem?mr=1680879>).
+This has slightly more restrictive conditions than the strong Lucas test,
+but uses different starting parameters so is not directly comparable.
+
=head2 is_aks_prime
Takes a positive number as input, and returns 1 if the input can be proven
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index 126a106..4acd817 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -3,7 +3,9 @@ use strict;
use warnings;
use Test::More;
-use Math::Prime::Util qw/is_prime is_strong_pseudoprime is_strong_lucas_pseudoprime/;
+use Math::Prime::Util qw/is_prime
+ is_pseudoprime is_strong_pseudoprime
+ is_strong_lucas_pseudoprime is_extra_strong_lucas_pseudoprime/;
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -57,7 +59,10 @@ my %pseudoprimes = (
1795265022 => [ qw/1795265023 1797174457 1797741901 1804469753 1807751977 1808043283 1808205701 1813675681 1816462201 1817936371 1819050257/ ],
3046413974 => [ qw/3046413975 3048698683 3051199817 3068572849 3069705673 3070556233 3079010071 3089940811 3090723901 3109299161 3110951251 3113625601/ ],
3613982119 => [ qw/3626488471 3630467017 3643480501 3651840727 3653628247 3654142177 3672033223 3672036061 3675774019 3687246109 3690036017 3720856369/ ],
- lucas => [ qw/5459 5777 10877 16109 18971 22499 24569 25199 40309 58519 75077 97439/ ],
+ psp2 => [ qw/341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 8481 8911 10261 10585 11305 12801 13741 13747 13981 14491 15709 15841 16705 18705 18721 19951 23001 23377 25761 29341/ ],
+ psp3 => [ qw/91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 7381 8401 8911 10585 11011 12403 14383 15203 15457 15841 16471 16531 18721 19345 23521 24046 24661 24727 28009 29161/ ],
+ slucas => [ qw/5459 5777 10877 16109 18971 22499 24569 25199 40309 58519 75077 97439 100127 113573 115639 130139/ ],
+ eslucas => [ qw/989 3239 5777 10877 27971 29681 30739 31631 39059 72389 73919 75077 100127 113573 125249 137549 137801/ ],
);
my $num_pseudoprimes = 0;
foreach my $ppref (values %pseudoprimes) {
@@ -87,8 +92,13 @@ is( is_strong_pseudoprime(3, 2), 1, "MR with 3 shortcut prime");
# Check that each strong pseudoprime base b makes it through MR with that base
while (my($base, $ppref) = each (%pseudoprimes)) {
foreach my $p (@$ppref) {
- if ($base eq 'lucas') {
+ if ($base eq 'eslucas') {
+ ok(is_extra_strong_lucas_pseudoprime($p), "$p is an extra strong Lucas pseudoprime");
+ } elsif ($base eq 'slucas') {
ok(is_strong_lucas_pseudoprime($p), "$p is a strong Lucas-Selfridge pseudoprime");
+ } elsif ($base =~ /^psp(\d+)/) {
+ my $base = $1;
+ ok(is_pseudoprime($p, $base), "$p is a pseudoprime to base $base");
} else {
ok(is_strong_pseudoprime($p, $base), "Pseudoprime (base $base) $p passes MR");
}
--
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