[libmath-prime-util-perl] 07/29: Convenient primality proof random 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.27
in repository libmath-prime-util-perl.

commit de8f5aa09db4c7b31667930b3ef251153a44ce8a
Author: Dana Jacobsen <dana at acm.org>
Date:   Wed May 8 18:31:02 2013 -0700

    Convenient primality proof random test
---
 lib/Math/Prime/Util/PP.pm | 29 ++++++++++++++-----
 xt/primality-proofs.pl    | 72 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 94 insertions(+), 7 deletions(-)

diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index a1349f7..de34252 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1875,12 +1875,22 @@ sub primality_proof_bls75 {
   my $B = $nm1->copy;         # unfactored part
   my @factors = (2);
   croak "BLS75 error: n-1 not even" unless $nm1->is_even();
-  while ($B->is_even) { $B /= 2; $A *= 2; }
-  foreach my $f (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79) {
-    last if $f*$f > $B;
-    if (($B % $f) == 0) {
-      push @factors, $f;
-      do { $B /= $f;  $A *= $f; } while (($B % $f) == 0);
+  {
+    while ($B->is_even) { $B /= 2; $A *= 2; }
+    my $tlim = $_primes_small[-1];
+    if ($B > $tlim*$tlim) {
+      my @small_primes = @_primes_small[2 .. $#_primes_small];
+      foreach my $f (@small_primes) {
+        if (($B % $f) == 0) {
+          push @factors, $f;
+          do { $B /= $f;  $A *= $f; } while (($B % $f) == 0);
+          last if $B <= $tlim*$tlim;
+        }
+      }
+    }
+    if ($B <= $tlim*$tlim) {
+      push @factors, factor($B);
+      $A *= $B;  $B /= $B;
     }
   }
   my @nstack;
@@ -1893,10 +1903,15 @@ sub primality_proof_bls75 {
   } else {
     push @nstack, $B;
   }
+  my $using_theorem_7 = 0;
   while (@nstack) {
     my ($s,$r) = $B->copy->bdiv($A->copy->bmul(2));
     my $fpart = ($A+1) * (2*$A*$A + ($r-1) * $A + 1);
     last if $n < $fpart;
+    # Theorem 7....
+    # my $tlim = $_primes_small[-1];
+    # $fpart = ($A*$tlim+1) * (2*$A*$A + ($r-$tlim) * $A + 1);
+    # if ($n < $fpart) { $using_theorem_7 = 1; last; }
 
     my $m = pop @nstack;
     # Don't use bignum if it has gotten small enough.
@@ -1905,8 +1920,8 @@ sub primality_proof_bls75 {
     my @ftry;
     $_holf_r = 1;
     foreach my $sub (@_fsublist) {
-      last if scalar @ftry >= 2;
       @ftry = $sub->($m);
+      last if scalar @ftry >= 2;
     }
     # If we couldn't find a factor, skip it.
     next unless scalar @ftry > 1;
diff --git a/xt/primality-proofs.pl b/xt/primality-proofs.pl
new file mode 100644
index 0000000..580be8d
--- /dev/null
+++ b/xt/primality-proofs.pl
@@ -0,0 +1,72 @@
+#!/usr/bin/perl
+use warnings;
+use strict;
+use Math::Prime::Util ':all';
+use Math::BigInt lib=>"GMP";
+#use Crypt::Primes 'maurer';
+use Math::Pari qw/isprime/;
+use Crypt::Random 'makerandom';
+use Data::Dump 'dump';
+
+$|++;
+my $num = 71;
+my $size = 300;
+my $prime_method = 'pari';   # mpu, pari, or cpmaurer
+
+my @ns;
+print "Generate ";
+die "invalid size, must be > 4" unless $size > 4;
+foreach my $i (1..$num) {
+  my $bits = int(rand($size-4)) + 4;
+  my $n;
+
+  # How do we get random primes?
+  #   MPU is the fastest, but it's our own code with identical primality tests.
+  #   Pari + Crypt::Random works pretty well if you have them.
+  #   Crypt::Primes::maurer will sometimes output composites (!!!).
+  if      ($prime_method eq 'cpmaurer') {
+    $n = Crypt::Primes::maurer(Size=>$bits); 
+  } elsif ($prime_method eq 'pari') {
+    do { $n = makerandom(Size=>$bits,Strength=>0); } while !isprime($n);
+  } elsif ($prime_method eq 'mpu') {
+    $n = random_nbit_prime($bits);
+  } else {
+    die "Unknown random prime generation method\n";
+  }
+
+  push @ns, Math::BigInt->new("$n");
+  # print a number corresponding to hundreds of bits
+  print int(3.322*length("$n")/100);
+}
+print "\n";
+
+my @certs;
+print "Prove    ";
+foreach my $n (@ns) {
+  my ($isp,$cert) = is_provable_prime_with_cert($n);
+  die "$n is reported as $isp\n" unless $isp == 2;
+  push @certs, [$n, $cert];
+  print proof_mark($cert);
+}
+print "\n";
+
+print "Verify   ";
+foreach my $certn (@certs) {
+  my $v = verify_prime($certn->[1]);
+  print proof_mark($certn->[1]);
+  next if $v;
+  print "\n\n$certn->[0] didn't verify!\n\n";
+  print dump($certn->[1]);
+  die;
+}
+print "\n";
+
+sub proof_mark {
+  my $cert = shift;
+  my $type = (scalar @$cert == 1) ? "bpsw" : $cert->[1];
+  if (!defined $type) { die "\nNo cert:\n\n", dump($cert); }
+  if    ($type =~ /n-1/i)       { return '1'; }
+  elsif ($type =~ /bpsw/i)      { return '.'; }
+  elsif ($type =~ /ecpp|agkm/i) { return 'E'; }
+  return '?';
+}

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