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 <d...@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(",",@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

_______________________________________________
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

Reply via email to