[libmath-prime-util-perl] 09/40: Add PP code and simple test for frobenius_underwood probable prime test

[Partha P. Mukherjee]

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 ea22a6b79afaa38c4cd6fe59e0dee45b1d6fbb8d
Author: Dana Jacobsen <dana at acm.org>
Date:   Mon Jun 24 15:16:27 2013 -0700

    Add PP code and simple test for frobenius_underwood probable prime test
---
 TODO                      |  5 -----
 lib/Math/Prime/Util/PP.pm | 45 +++++++++++++++++++++++++++++++++++++++++++++
 t/17-pseudoprime.t        | 24 +++++++++++++++++++++---
 3 files changed, 66 insertions(+), 8 deletions(-)

diff --git a/TODO b/TODO
index f9c19c0..33f6bc5 100644
--- a/TODO
+++ b/TODO
@@ -48,8 +48,3 @@
 
 - Write a standalone function that demonstrates the memory leak with MULTICALL,
   so we can use MULTICALL.
-
-- Tests for:
-   - F-U pseudoprime
-
-- PP F-U pseudoprime
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 435e5df..8aab410 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1090,6 +1090,51 @@ sub is_extra_strong_lucas_pseudoprime {
   return 0;
 }
 
+sub is_frobenius_underwood_pseudoprime {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return 0 if $n < 2;
+  return 1 if $n < 4;
+  return 0 if ($n % 2) == 0;
+  return 0 if _is_perfect_square($n);
+
+  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 $ZERO = $n->copy->bzero;
+  my $a = $ZERO + 1;
+  my $b = $ZERO + 2;
+
+  my ($x, $t, $np1, $len, $na) = (0, -1, $n+1, 1, undef);
+  while ( _jacobi($t, $n) != -1 ) {
+    $x++;
+    $t = $x*$x - 4;
+  }
+  my $result = $x+$x+5;
+  my $multiplier = $x+2;
+  $result %= $n if $result > $n;
+  $multiplier %= $n if $multiplier > $n;
+  { my $v = $np1; $len++ while ($v >>= 1); }
+  foreach my $bit (reverse 0 .. $len-2) {
+    $na = $a * (($a*$x) + ($b+$b));
+    $b = ( ($b + $a) * ($b - $a) ) % $n;
+    $a = $na % $n;
+    if ( ($np1 >> $bit) & 1 ) {
+      $na = $b + ($a * $multiplier);
+      $b += ($b - $a);
+      $a = $na;
+    }
+  }
+  $a->bmod($n);
+  $b->bmod($n);
+  return ($a == 0 && $b == $result) ? 1 : 0;
+}
+
 
 my $_poly_bignum;
 sub _poly_new {
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index 3a7caed..04d8d87 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -4,9 +4,12 @@ use warnings;
 
 use Test::More;
 use Math::Prime::Util qw/is_prime
-                         is_pseudoprime is_strong_pseudoprime
-                         is_lucas_pseudoprime is_strong_lucas_pseudoprime
+                         is_pseudoprime
+                         is_strong_pseudoprime
+                         is_lucas_pseudoprime
+                         is_strong_lucas_pseudoprime
                          is_extra_strong_lucas_pseudoprime
+                         is_frobenius_underwood_pseudoprime
                          lucas_sequence/;
 
 my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -92,6 +95,7 @@ plan tests => 0 + 3
                 + 9  # mr with large bases
                 + scalar @small_lucas_trials
                 + scalar(keys %lucas_sequences)
+                + 1  # frob-underwood
                 + 1*$extra;
 
 ok(!eval { is_strong_pseudoprime(2047); }, "MR with no base fails");
@@ -176,7 +180,7 @@ if ($extra) {
       last;
     }
   }
-  is($mr2fail, 0, "is_strong_pseudoprime bases 2,3 matches is_prime to 1,373,652");
+  is($mr2fail, 0, "is_strong_pseudoprime bases 2,3 matches is_prime");
 }
 
 # Lucas sequences, used for quite a few tests
@@ -186,3 +190,17 @@ sub lucas_sequence_to_native {
 while (my($params, $expect) = each (%lucas_sequences)) {
   is_deeply( [lucas_sequence_to_native(split(' ', $params))], $expect, "Lucas sequence $params" );
 }
+
+{
+  my $fufail = 0;
+  foreach my $i (1 .. 5000) {
+    my $n = int(rand(1000000000)) + 1;
+    my $ispfu = !!is_frobenius_underwood_pseudoprime($n);
+    my $prime = !!is_prime($n);
+    if ($ispfu != $prime) {
+      $fufail = $n;
+      last;
+    }
+  }
+  is($fufail, 0, "is_frobenius_underwood_pseudoprime matches is_prime");
+}

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git