[libmath-prime-util-perl] 23/59: Add BPSW primality test, hooray
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:55 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.10
in repository libmath-prime-util-perl.
commit 61e1c652a3defd14b77d8895a9d6ea70d240091d
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jul 2 21:33:08 2012 -0600
Add BPSW primality test, hooray
---
Changes | 1 +
TODO | 4 +-
XS.xs | 2 +-
factor.c | 2 +-
factor.h | 2 +-
lib/Math/Prime/Util.pm | 140 ++++++++++++++++++-----
lib/Math/Prime/Util/PP.pm | 279 +++++++++++++++++++++++++++++++---------------
t/80-pp.t | 7 +-
t/81-bignum.t | 2 +-
util.c | 4 +-
10 files changed, 316 insertions(+), 127 deletions(-)
diff --git a/Changes b/Changes
index 741dec2..811590c 100644
--- a/Changes
+++ b/Changes
@@ -10,6 +10,7 @@ Revision history for Perl extension Math::Prime::Util.
- factor and all_factor should correctly work with bigints.
- many more bigint/bignum changes
- factor always returns sorted results
+ - added BPSW primality test for large is_prob_prime
0.09 25 June 2012
- Pure Perl code added. Passes all tests. Used only if the XSLoader
diff --git a/TODO b/TODO
index 9736f7f..dcefb57 100644
--- a/TODO
+++ b/TODO
@@ -41,4 +41,6 @@
installed and call that, but I'm not sure how much it will help if we get
the above running.
-- Must add BPSW if we want is_prime to be meaningful for >64-bit.
+- Tune the Lucas sequence code, especially try using more BigInt features.
+
+- Add tests for strong lucas pseudoprime, perhaps export the function
diff --git a/XS.xs b/XS.xs
index c9f93ec..1e8ecfa 100644
--- a/XS.xs
+++ b/XS.xs
@@ -369,7 +369,7 @@ miller_rabin(IN UV n, ...)
RETVAL
int
-is_prob_prime(IN UV n)
+_XS_is_prob_prime(IN UV n)
double
_XS_ExponentialIntegral(double x)
diff --git a/factor.c b/factor.c
index ce55363..8e74e04 100644
--- a/factor.c
+++ b/factor.c
@@ -274,7 +274,7 @@ int miller_rabin(UV n, const UV *bases, int nbases)
return 1;
}
-int is_prob_prime(UV n)
+int _XS_is_prob_prime(UV n)
{
UV bases[12];
int nbases;
diff --git a/factor.h b/factor.h
index d474203..c09ca19 100644
--- a/factor.h
+++ b/factor.h
@@ -20,7 +20,7 @@ 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, int nbases);
-extern int is_prob_prime(UV n);
+extern int _XS_is_prob_prime(UV n);
static UV gcd_ui(UV x, UV y) {
UV t;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 171371a..706f559 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -59,7 +59,6 @@ BEGIN {
# Perhaps these can be moved here?
*miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
- *is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime;
# These probably shouldn't even be exported
*trial_factor = \&Math::Prime::Util::PP::trial_factor;
@@ -87,6 +86,24 @@ if ($_Config{'maxbits'} == 32) {
$_Config{'maxprimeidx'} = 425656284035217743;
}
+# Notes on how we're dealing with big integers:
+#
+# 1) if (ref($n) eq 'Math::BigInt')
+# $n is a bigint, so do bigint stuff
+#
+# 2) if (defined $bigint::VERSION && $n > ~0)
+# 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.
+#
+# 3) if (ref($n) =~ /^Math::Big/)
+# $n is a big int, float, or rat. We probably want this as an int.
+#
+# $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
+# get us out of big math if we can
+
+
sub prime_get_config {
my %config = %_Config;
@@ -104,11 +121,16 @@ sub _validate_positive_integer {
croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
- if ($n > $_Config{'maxparam'}) {
- croak "Parameter '$n' outside of integer range" if !defined $Math::BigInt::VERSION;
- # Make $n a proper object if it isn't one already (or convert from float)
- $_[0] = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
+ if ($n <= $_Config{'maxparam'}) {
+ $_[0] = $n->as_number() if ref($n) eq 'Math::BigFloat';
+ $_[0] = $n->numify() if ref($n) eq 'Math::BigInt';
+ } elsif (ref($n) ne 'Math::BigInt') {
+ croak "Parameter '$n' outside of integer range" if !defined $bigint::VERSION;
+ $_[0] = Math::BigInt->new("$n"); # Make $n a proper bigint object
}
+ # One of these will be true:
+ # 1) $n <= max and $n is not a bigint
+ # 2) $n > max and $n is a bigint
1;
}
@@ -267,7 +289,11 @@ my @_random_ndigit_ranges;
sub random_ndigit_prime {
my($digits) = @_;
- _validate_positive_integer($digits, 1, (defined $bigint::VERSION) ? 10000 : $_Config{'maxdigits'} );
+ if (defined $bigint::VERSION && bigint::in_effect()) {
+ _validate_positive_integer($digits, 1, 10000);
+ } else {
+ _validate_positive_integer($digits, 1, $_Config{'maxdigits'});
+ }
if (!defined $_random_ndigit_ranges[$digits]) {
my $low = ($digits == 1) ? 0 : int(10 ** ($digits-1));
my $high = int(10 ** $digits);
@@ -395,13 +421,69 @@ sub factor {
return _XS_factor($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
- $n = $n->as_number() if ref($n) eq 'Math::BigFloat';
- $n = $n->numify if $n <= ~0;
Math::Prime::Util::PP::factor($n);
}
#############################################################################
+ # Timings for various combinations, given the current possibilities of:
+ # 1) XS MR optimized (either x86-64, 32-bit on 64-bit mach, or half-word)
+ # 2) XS MR non-optimized (big input not on 64-bit machine)
+ # 3) PP MR with small input (non-bigint Perl)
+ # 4) PP MR with large input (using functions for mulmod)
+ # 5) PP MR with full bigints
+ # 6) PP Lucas with small input
+ # 7) PP Lucas with large input
+ # 8) PP Lucas with full bigints
+ #
+ # Time for one test:
+ # 0.5uS XS MR with small input
+ # 0.8uS XS MR with large input
+ # 7uS PP MR with small input
+ # 400uS PP MR with large input
+ # 5000uS PP MR with bigint
+ # 2700uS PP LP with small input
+ # 6100uS PP LP with large input
+ # 7400uS PP LP with bigint
+
+sub is_prob_prime {
+ my($n) = @_;
+ return 0 if defined $n && $n < 2;
+ _validate_positive_integer($n);
+
+ return _XS_is_prob_prime($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+
+ return 2 if $n == 2 || $n == 3 || $n == 5 || $n == 7;
+ return 0 if $n < 11;
+ return 0 if ($n % 2) == 0 || ($n % 3) == 0 || ($n % 5) == 0 || ($n % 7) == 0;
+ foreach my $i (qw/11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71/) {
+ return 2 if $i*$i > $n; return 0 if ($n % $i) == 0;
+ }
+
+ if ($n < 105936894253) { # BPSW seems to be faster after this
+ # Deterministic set of Miller-Rabin tests.
+ my @bases;
+ if ($n < 9080191) { @bases = (31, 73); }
+ elsif ($n < 4759123141) { @bases = (2, 7, 61); }
+ elsif ($n < 105936894253) { @bases = (2, 1005905886, 1340600841); }
+ elsif ($n < 31858317218647) { @bases = (2, 642735, 553174392, 3046413974); }
+ elsif ($n < 3071837692357849) { @bases = (2, 75088, 642735, 203659041, 3613982119); }
+ else { @bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022); }
+ return Math::Prime::Util::PP::miller_rabin($n, @bases) ? 2 : 0;
+ }
+
+ # BPSW probable prime. No composites are known to have passed this test
+ # since it was published in 1980, though we know infinitely many exist.
+ # It has also been verified that no 64-bit composite will return true.
+ # Slow since it's all in PP, but it's the Right Thing To Do.
+
+ return 0 unless Math::Prime::Util::PP::miller_rabin($n, 2);
+ return 0 unless Math::Prime::Util::PP::is_strong_lucas_pseudoprime($n);
+ return ($n <= 18446744073709551615) ? 2 : 1;
+}
+
+#############################################################################
+
sub prime_count_approx {
my($x) = @_;
_validate_positive_integer($x);
@@ -409,7 +491,7 @@ sub prime_count_approx {
return $_prime_count_small[$x] if $x <= $#_prime_count_small;
# Turn on high precision FP if they gave us a big number.
- #do { require Math::BigFloat; Math::BigFloat->import; } if defined $Math::BigInt::VERSION && !defined $Math::BigFloat::VERSION;
+ $x = _upgrade_to_float($x) if ref($x) eq 'Math::BigInt';
# Method 10^10 %error 10^19 %error
# ----------------- ------------ ------------
@@ -433,7 +515,7 @@ sub prime_count_lower {
return $_prime_count_small[$x] if $x <= $#_prime_count_small;
- $x = _upgrade_to_float($x) if defined $Math::BigInt::VERSION && ref($x) ne 'Math::BigFloat';
+ $x = _upgrade_to_float($x) if ref($x) eq 'Math::BigInt';
my $flogx = log($x);
@@ -468,7 +550,7 @@ sub prime_count_upper {
return $_prime_count_small[$x] if $x <= $#_prime_count_small;
- $x = _upgrade_to_float($x) if defined $Math::BigInt::VERSION && ref($x) ne 'Math::BigFloat';
+ $x = _upgrade_to_float($x) if ref($x) eq 'Math::BigInt';
# Chebyshev: 1.25506*x/logx x >= 17
# Rosser & Schoenfeld: x/(logx-3/2) x >= 67
@@ -520,7 +602,7 @@ sub nth_prime_approx {
return $_primes_small[$n] if $n <= $#_primes_small;
- $n = _upgrade_to_float($n) if defined $Math::BigInt::VERSION && ref($n) ne 'Math::BigFloat';
+ $n = _upgrade_to_float($n) if ref($n) eq 'Math::BigInt';
my $flogn = log($n);
my $flog2n = log($flogn);
@@ -572,7 +654,7 @@ sub nth_prime_lower {
return $_primes_small[$n] if $n <= $#_primes_small;
- $n = _upgrade_to_float($n) if defined $Math::BigInt::VERSION && ref($n) ne 'Math::BigFloat';
+ $n = _upgrade_to_float($n) if ref($n) eq 'Math::BigInt';
my $flogn = log($n);
my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n)
@@ -595,7 +677,7 @@ sub nth_prime_upper {
return $_primes_small[$n] if $n <= $#_primes_small;
- $n = _upgrade_to_float($n) if defined $Math::BigInt::VERSION && ref($n) ne 'Math::BigFloat';
+ $n = _upgrade_to_float($n) if ref($n) eq 'Math::BigInt';
my $flogn = log($n);
my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n)
@@ -629,7 +711,7 @@ sub RiemannR {
my($n) = @_;
croak("Invalid input to ReimannR: x must be > 0") if $n <= 0;
- return Math::Prime::Util::PP::RiemannR($n, 1e-30) if defined $Math::BigFloat::VERSION;
+ return Math::Prime::Util::PP::RiemannR($n, 1e-30) if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat';
return Math::Prime::Util::PP::RiemannR($n) if !$_Config{'xs'};
return _XS_RiemannR($n);
@@ -644,7 +726,7 @@ sub ExponentialIntegral {
my($n) = @_;
croak "Invalid input to ExponentialIntegral: x must be != 0" if $n == 0;
- return Math::Prime::Util::PP::ExponentialIntegral($n, 1e-30) if defined $Math::BigFloat::VERSION;
+ return Math::Prime::Util::PP::ExponentialIntegral($n, 1e-30) if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat';
return Math::Prime::Util::PP::ExponentialIntegral($n) if !$_Config{'xs'};
return _XS_ExponentialIntegral($n);
}
@@ -654,7 +736,7 @@ sub LogarithmicIntegral {
return 0 if $n == 0;
croak("Invalid input to LogarithmicIntegral: x must be >= 0") if $n <= 0;
- if (defined $Math::BigFloat::VERSION) {
+ if ( (defined $bignum::VERSION && (!defined &bignum::in_effect || bignum::in_effect())) || (ref($n) eq 'Math::BigFloat') ) {
return Math::BigFloat->binf('-') if $n == 1;
return Math::BigFloat->new('1.045163780117492784844588889194613136522615578151201575832909144075013205210359530172717405626383356306') if $n == 2;
} else {
@@ -797,6 +879,7 @@ your program.
A number of the functions support big numbers, but currently not all. The
ones that do:
+ is_prob_prime
prime_count_lower
prime_count_upper
prime_count_approx
@@ -814,7 +897,6 @@ ones that do:
These still do not:
is_prime
- is_prob_prime
miller_rabin
primes
next_prime
@@ -1019,11 +1101,19 @@ function.
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
-will 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.
+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 inputs larger than C<2^64>, a strong Baillie-PSW primality test is
+performed (aka BPSW or BSW). This is a probabilistic test, so the only times
+a 2 (definitely prime) are returned are when the small trial division succeeds.
+Note that since the test was published in 1980, not a single BPSW pseudoprime
+has been found, so it is extremely likely to be prime. While we know there
+an infinite number of counterexamples exist, there is a weak conjecture that
+none exist under 10000 digits.
=head2 random_prime
@@ -1169,7 +1259,7 @@ process is repeated for each non-prime factor.
While factoring works on bigints, the algorithms are currently set up for
smaller numbers, and bignum support is all in pure Perl. Hence, it will be
-somewhat slow for "easy" numbers and too slow for "hard" numbers.
+somewhat slow for "easy" numbers and very, very slow for "hard" numbers.
=head2 all_factors
@@ -1302,7 +1392,7 @@ Print pseudoprimes base 17:
Print some primes above 64-bit range:
perl -MMath::Prime::Util=:all -Mbigint -E 'my $start=100000000000000000000; say join "\n", @{primes($start,$start+1000)}'
- # Similar behavior:
+ # Similar but much faster:
# perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'
=head1 LIMITATIONS
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 684ecb0..e72ebc7 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -53,6 +53,26 @@ sub _is_positive_int {
((defined $_[0]) && ($_[0] !~ tr/0123456789//c));
}
+sub _validate_positive_integer {
+ my($n, $min, $max) = @_;
+ croak "Parameter must be defined" if !defined $n;
+ croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
+ croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
+ croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
+ if ($n <= ~0) {
+ $_[0] = $n->as_number() if ref($n) eq 'Math::BigFloat';
+ $_[0] = $n->numify() if ref($n) eq 'Math::BigInt';
+ } elsif (ref($n) ne 'Math::BigInt') {
+ croak "Parameter '$n' outside of integer range" if !defined $bigint::VERSION;
+ $_[0] = Math::BigInt->new("$n"); # Make $n a proper bigint object
+ }
+ # One of these will be true:
+ # 1) $n <= max and $n is not a bigint
+ # 2) $n > max and $n is a bigint
+ 1;
+}
+
+
my @_primes_small = (
0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
@@ -81,7 +101,7 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5
return 0 if ($n % $i) == 0;
}
- return is_prob_prime($n) if $n > 10_000_000;
+ return Math::Prime::Util::is_prob_prime($n) if $n > 10_000_000;
my $limit = int(sqrt($n));
my $i = 31;
@@ -110,13 +130,7 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5
sub is_prime {
my($n) = @_;
- croak "Input must be an integer" unless defined $_[0];
- if ($n > ~0) {
- croak "Input out of range" if !defined $Math::BigInt::VERSION;
- $n = Math::BigInt->new($n) unless ref($n) =~ /^Math::Big/;
- } else {
- $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
- }
+ _validate_positive_integer($n);
return 2 if ($n == 2) || ($n == 3) || ($n == 5); # 2, 3, 5 are prime
return 0 if $n < 7; # everything else below 7 is composite
@@ -245,8 +259,8 @@ sub trial_primes {
$high = $low;
$low = 2;
}
- croak "Input must be a positive integer" unless _is_positive_int($low)
- && _is_positive_int($high);
+ _validate_positive_integer($low);
+ _validate_positive_integer($high);
return if $low > $high;
@@ -268,15 +282,13 @@ sub primes {
my $high = shift;
my $sref = [];
- croak "Input must be a positive integer" unless _is_positive_int($low)
- && _is_positive_int($high);
+ _validate_positive_integer($low);
+ _validate_positive_integer($high);
return $sref if ($low > $high) || ($high < 2);
# Ignore method options in this code
- $low = $low->numify if ref($low) =~ /^Math::Big/ && $low <= ~0;
- $high = $high->numify if ref($high) =~ /^Math::Big/ && $high <= ~0;
# At some point even the pretty-fast pure perl sieve is going to be a
# dog, and we should move to trials. This is typical with a small range
# on a large base. More thought on the switchover should be done.
@@ -312,17 +324,10 @@ sub primes {
sub next_prime {
my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
- return 0 if ($n >= ((_PP_prime_maxbits == 32) ? 4294967291 : 18446744073709551557))
- && (!defined $Math::BigInt::VERSION);
+ _validate_positive_integer($n);
+ return 0 if ($n >= ((_PP_prime_maxbits == 32) ? 4294967291 : 18446744073709551557)) && ref($n) ne 'Math::BigInt';
return $_prime_next_small[$n] if $n <= $#_prime_next_small;
- if (ref($n) =~ /^Math::Big/) {
- $n = $n->numify if $n <= ~0;
- } elsif ($n > ~0) {
- $n = Math::BigInt->new($n);
- }
-
# Be careful trying to do:
# my $d = int($n/30);
# my $m = $n - $d*30;
@@ -339,7 +344,7 @@ sub next_prime {
sub prev_prime {
my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
if ($n <= 7) {
return ($n <= 2) ? 0 : ($n <= 3) ? 2 : ($n <= 5) ? 3 : 5;
}
@@ -376,23 +381,8 @@ sub prime_count {
$high = $low;
$low = 2;
}
- croak "Input must be a positive integer" unless _is_positive_int($low)
- && _is_positive_int($high);
-
-if (0) {
- if ($low > ~0) {
- croak "Input out of range" if !defined $Math::BigInt::VERSION;
- $low = Math::BigInt->new($low) unless ref($low) =~ /^Math::Big/;
- } else {
- $low = $low->numify if $low < ~0 && ref($low) =~ /^Math::Big/;
- }
- if ($high > ~0) {
- croak "Input out of range" if !defined $Math::BigInt::VERSION;
- $high = Math::BigInt->new($high) unless ref($high) =~ /^Math::Big/;
- } else {
- $high = $high->numify if $high < ~0 && ref($high) =~ /^Math::Big/;
- }
-}
+ _validate_positive_integer($low);
+ _validate_positive_integer($high);
my $count = 0;
@@ -431,11 +421,11 @@ if (0) {
sub nth_prime {
my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
return $_primes_small[$n] if $n <= $#_primes_small;
- if (!defined $bigint::VERSION) {
+ if (!defined $bigint::VERSION) { # This isn't ideal.
if (_PP_prime_maxbits == 32) {
croak "nth_prime($n) overflow" if $n > 203280221;
} else {
@@ -535,16 +525,9 @@ sub _gcd_ui {
sub miller_rabin {
my($n, @bases) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
croak "No bases given to miller_rabin" unless @bases;
- if ($n > ~0) {
- croak "Input out of range" if !defined $Math::BigInt::VERSION;
- $n = Math::BigInt->new($n) unless ref($n) =~ /^Math::Big/;
- } else {
- $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
- }
-
return 0 if ($n == 0) || ($n == 1);
return 2 if ($n == 2) || ($n == 3);
return 0 if ($n % 2) == 0;
@@ -561,7 +544,7 @@ sub miller_rabin {
$d >>= 1;
}
- if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+ if ( ($n < $_half_word) || (ref($n) eq 'Math::BigInt') ) {
foreach my $a (@bases) {
croak "Base $a is invalid" if $a < 2;
@@ -600,39 +583,157 @@ sub miller_rabin {
1;
}
-sub is_prob_prime {
+# Calculate Jacobi symbol (M|N)
+sub _jacobi {
+ my($n, $m) = @_;
+ return 0 if $m <= 0 || ($m % 2) == 0;
+ my $j = 1;
+ if ($n < 0) {
+ $n = -$n;
+ $j = -$j if ($m % 4) == 3;
+ }
+ while ($n != 0) {
+ while (($n % 2) == 0) {
+ $n >>= 1;
+ $j = -$j if ($m % 8) == 3 || ($m % 8) == 5;
+ }
+ ($n, $m) = ($m, $n);
+ $j = -$j if ($n % 4) == 3 && ($m % 4) == 3;
+ $n = $n % $m;
+ }
+ return ($m == 1) ? $j : 0;
+}
+
+# Find first D in sequence (5,-7,9,-11,13,-15,...) where (D|N) == -1
+sub _find_jacobi_d_sequence {
my($n) = @_;
- if ($n > ~0) {
- croak "Input out of range" if !defined $Math::BigInt::VERSION;
- $n = Math::BigInt->new($n) unless ref($n) =~ /^Math::Big/;
- } else {
- $n = $n->numify if $n < ~0 && ref($n) =~ /^Math::Big/;
+ # D is typically quite small: 67 max for N < 10^19. However, it is
+ # theoretically possible D could grow unreasonably. We are ignoring this.
+ my $d = 5;
+ my $sign = 1;
+ while (1) {
+ my $gcd = _gcd_ui($d, $n);
+ return 0 if $gcd > 1 && $gcd != $n; # Found divisor $d
+ my $j = _jacobi($d * $sign, $n);
+ last if $j == -1;
+ $d += 2;
+ $sign = -$sign;
+ }
+ return ($sign * $d);
+}
+
+
+sub is_strong_lucas_pseudoprime {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ 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
+ {
+ my $mcheck = $n & 127;
+ if (($mcheck*0x8bc40d7d) & ($mcheck*0xa1e2f5d1) & 0x14020a) {
+ # ... 82% of non-squares were rejected by the bloom filter
+ # For bigints we should put in the rest of this filter, which prunes
+ # 99.92% of nonsquares, for the cost of one bigint mod and some
+ # non-bigint operations.
+ my $sq = int(sqrt($n));
+ return 0 if ($sq*$sq) == $n;
+ }
}
- return 2 if ($n == 2) || ($n == 3) || ($n == 5) || ($n == 7);
- return 0 if ($n < 2) || (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0) || (($n % 7) == 0);
- return 2 if $n < 121;
-
- my @bases;
- # We can do this, if MR accepts bases larger than the number
- #if ($n < 316349281) { @bases = (11000544, 31481107); }
- if ($n < 9080191) { @bases = (31, 73); }
- elsif ($n < 4759123141) { @bases = (2, 7, 61); }
- elsif ($n < 105936894253) { @bases = (2, 1005905886, 1340600841); }
- elsif ($n < 31858317218647) { @bases = (2, 642735, 553174392, 3046413974); }
- elsif ($n < 3071837692357849) { @bases = (2, 75088, 642735, 203659041, 3613982119); }
- elsif ($n < 18446744073709551615) { @bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022); }
- # Seriously we need BPSW here.
- else { @bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022, 15, 52, 311); }
-
- # Run our selected set of Miller-Rabin strong probability tests
- my $prob_prime = miller_rabin($n, @bases);
- # We're deterministic for 64-bit numbers
- $prob_prime *= 2 if $n <= 18446744073709551615;
- $prob_prime;
+ # We have to make sure we use bigints
+ my $result = 0;
+ {
+ use bigint;
+ $n = Math::BigInt->new($n) unless ref($n) eq 'Math::BigInt';
+
+ # Determine Selfridge D, P, and Q parameters
+ my $D = _find_jacobi_d_sequence($n);
+ $D = $D->numify if $D <= ~0 && ref($D) eq 'Math::BigInt';
+ return 0 if $D == 0; # We found a divisor in the sequence
+ my $P = 1;
+ my $Q = int( (1 - $D) / 4 );
+ # Verify we've calculated this right
+ die "Selfridge error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);
+ #warn "N: $n D: $D P: $P Q: $Q\n";
+
+ # It's now time to perform the Lucase pseudoprimality test using $D.
+
+ my $m = $n + 1;
+ my $s = 0;
+ my $d = $m;
+ while ( ($d & 1) == 0 ) {
+ $s++;
+ $d >>= 1;
+ }
+ # $m == $d * 2 ** $s
+ #die "Invalid $m, $d, $s\n" unless $m == $d * 2**$s;
+
+ my $U = 1;
+ my $V = $P;
+ my $U_2m = 1;
+ my $V_2m = $P;
+ my $Q_m = $Q;
+ my $Q_m2 = 2 * $Q;
+ my $Qkd = $Q;
+ $d >>= 1;
+ #my $i = 1;
+ while ($d != 0) {
+ #warn "U=$U V=$V Qm=$Q_m Qm2=$Q_m2\n";
+ $U_2m = ($U_2m * $V_2m) % $n;
+ $V_2m = ($V_2m * $V_2m - $Q_m2) % $n;
+ #warn " $i U2m=$U_2m V2m=$V_2m\n"; $i++;
+ $Q_m = ($Q_m * $Q_m) % $n;
+ $Q_m2 = 2 * $Q_m; # no mod
+ if ($d & 1) {
+ my $Uold = $U;
+ $U = $U_2m * $V + $V_2m * $U;
+ $U += $n if $U & 1;
+ $U = int($U / 2) % $n;
+
+ $V = $V_2m * $V + $U_2m * $Uold * $D;
+ $V += $n if $V & 1;
+ $V = int($V / 2) % $n;
+
+ $Qkd = ($Qkd * $Q_m) % $n;
+ }
+ $d >>= 1;
+ }
+ #warn "l0 U=$U V=$V\n";
+ if ($U == 0 || $V == 0) {
+ $result = 1;
+ } else {
+ # Compute powers of V
+ my $Qkd2 = 2 * $Qkd;
+ foreach my $r (1 .. $s-1) {
+ #warn "l$r U=$U V=$V Qkd2=$Qkd2\n";
+ $V = ($V * $V - $Qkd2) % $n;
+ if ($V == 0) {
+ $result = 1;
+ last;
+ }
+ if ($r < ($s-1)) {
+ $Qkd = ($Qkd * $Qkd) % $n;
+ $Qkd2 = 2 * $Qkd;
+ }
+ }
+ }
+ no bigint;
+ }
+ #warn "Math::BigInt is loaded\n" if defined $Math::BigInt::VERSION;
+ #warn "bigint is loaded\n" if defined $bigint::VERSION;
+ #warn "bigint is in effect\n" if bigint::in_effect();
+ return $result;
}
+
sub _basic_factor {
# MODIFIES INPUT SCALAR
return ($_[0]) if $_[0] < 4;
@@ -643,7 +744,7 @@ sub _basic_factor {
while ( ($_[0] % 5) == 0 ) { push @factors, 5; $_[0] = int($_[0] / 5); }
# Stop using bignum if we can
- $_[0] = $_[0]->numify if ref($_[0]) =~ /^Math::Big/ && $_[0] <= ~0;
+ $_[0] = $_[0]->numify if ref($_[0]) eq 'Math::BigInt' && $_[0] <= ~0;
if ( ($_[0] > 1) && _is_prime7($_[0]) ) {
push @factors, $_[0];
@@ -654,7 +755,7 @@ sub _basic_factor {
sub trial_factor {
my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
my @factors = _basic_factor($n);
return @factors if $n < 4;
@@ -677,7 +778,7 @@ sub trial_factor {
sub factor {
my($n) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
return trial_factor($n) if $n < 100000;
@@ -701,7 +802,7 @@ sub factor {
while (@nstack) {
$n = pop @nstack;
# Don't use bignum on $n if it has gotten small enough.
- $n = $n->numify if ref($n) =~ /^Math::Big/ && $n <= ~0;
+ $n = $n->numify if ref($n) eq 'Math::BigInt' && $n <= ~0;
#print "Looking at $n with stack ", join(",", at nstack), "\n";
while ( ($n >= (31*31)) && !_is_prime7($n) ) {
my @ftry;
@@ -749,7 +850,7 @@ sub squfof_factor { trial_factor(@_) }
sub prho_factor {
my($n, $rounds, $a) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
$rounds = 4*1024*1024 unless defined $rounds;
$a = 3 unless defined $a;
@@ -760,7 +861,7 @@ sub prho_factor {
my $U = 7;
my $V = 7;
- if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+ if ( ($n < $_half_word) || (ref($n) eq 'Math::BigInt') ) {
for my $i (1 .. $rounds) {
$U = ($U * $U + $a) % $n;
@@ -807,7 +908,7 @@ sub prho_factor {
sub pbrent_factor {
my($n, $rounds) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
$rounds = 4*1024*1024 unless defined $rounds;
my @factors = _basic_factor($n);
@@ -817,7 +918,7 @@ sub pbrent_factor {
my $Xi = 2;
my $Xm = 2;
- if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+ if ( ($n < $_half_word) || (ref($n) eq 'Math::BigInt') ) {
for my $i (1 .. $rounds) {
$Xi = ($Xi * $Xi + $a) % $n;
@@ -854,7 +955,7 @@ sub pbrent_factor {
sub pminus1_factor {
my($n, $rounds) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
$rounds = 4*1024*1024 unless defined $rounds;
my @factors = _basic_factor($n);
@@ -879,7 +980,7 @@ sub pminus1_factor {
sub holf_factor {
my($n, $rounds, $startrounds) = @_;
- croak "Input must be a positive integer" unless _is_positive_int($n);
+ _validate_positive_integer($n);
$rounds = 64*1024*1024 unless defined $rounds;
$startrounds = 1 unless defined $startrounds;
diff --git a/t/80-pp.t b/t/80-pp.t
index a19dc19..99c5888 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -215,7 +215,7 @@ my %rvals = (
plan tests => 1 +
- 2*(1087 + @primes + @composites) +
+ 1*(1087 + @primes + @composites) +
3 + scalar(keys %small_single) + scalar(keys %small_range) +
2*scalar(keys %primegaps) + 8 + 148 + 148 + 1 +
scalar(keys %pivals_small) + scalar(keys %pi_intervals) +
@@ -239,7 +239,6 @@ require_ok 'Math::Prime::Util::PP';
*prev_prime = \&Math::Prime::Util::PP::prev_prime;
*miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
- *is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime;
*factor = \&Math::Prime::Util::PP::factor;
@@ -252,20 +251,16 @@ require_ok 'Math::Prime::Util::PP';
foreach my $n (0 .. 1086) {
if (defined $small_primes{$n}) {
is( is_prime($n), 2, "$n is prime");
- ok( is_prob_prime($n), "$n is probably prime");
} else {
ok(!is_prime($n), "$n is not prime");
- ok(!is_prob_prime($n), "$n is not probably prime");
}
}
foreach my $n (@primes) {
is( is_prime($n), 2, "$n is prime" );
- ok( is_prob_prime($n), "$n is probably prime");
}
foreach my $n (@composites) {
is( is_prime($n), 0, "$n is not prime" );
- is( is_prob_prime($n), 0, "$n is not probably prime");
}
###############################################################################
diff --git a/t/81-bignum.t b/t/81-bignum.t
index 2ed5b02..25cc484 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -63,6 +63,7 @@ plan tests => 0 +
0;
use Math::Prime::Util qw/
+ is_prob_prime
prime_count_lower
prime_count_upper
prime_count_approx
@@ -88,7 +89,6 @@ use Math::Prime::Util qw/
*prev_prime = \&Math::Prime::Util::PP::prev_prime;
*miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
- *is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime;
###############################################################################
diff --git a/util.c b/util.c
index f57c3d0..b508a33 100644
--- a/util.c
+++ b/util.c
@@ -57,7 +57,7 @@ static int _is_prime7(UV n)
UV limit, i;
if (n > MPU_PROB_PRIME_BEST)
- return is_prob_prime(n); /* We know this works for all 64-bit n */
+ return _XS_is_prob_prime(n); /* We know this works for all 64-bit n */
limit = sqrt(n);
i = 7;
@@ -142,7 +142,7 @@ int is_definitely_prime(UV n)
if (isprime >= 0) return isprime;
if (n > MPU_PROB_PRIME_BEST)
- return (is_prob_prime(n) == 2);
+ return (_XS_is_prob_prime(n) == 2);
return 0;
}
--
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