[libmath-prime-util-perl] 35/40: INCOMPLETE BROKEN : Switching to new text proofs, part 1
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:05 UTC 2015
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 046c5a9024016c614d49abe1cec901c02c3a41b6
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Aug 1 17:33:37 2013 -0700
INCOMPLETE BROKEN : Switching to new text proofs, part 1
---
MANIFEST | 1 +
TODO | 2 +
lib/Math/Prime/Util.pm | 406 +++++-----------------------------------------
lib/Math/Prime/Util/PP.pm | 198 ----------------------
t/23-primality-proofs.t | 32 ++--
t/80-pp.t | 10 +-
xt/primality-proofs.pl | 17 +-
7 files changed, 83 insertions(+), 583 deletions(-)
diff --git a/MANIFEST b/MANIFEST
index 3b39ca2..ab51514 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -7,6 +7,7 @@ lib/Math/Prime/Util/PP.pm
lib/Math/Prime/Util/ZetaBigFloat.pm
lib/Math/Prime/Util/ECAffinePoint.pm
lib/Math/Prime/Util/ECProjectivePoint.pm
+lib/Math/Prime/Util/PrimalityProving.pm
LICENSE
Makefile.PL
MANIFEST
diff --git a/TODO b/TODO
index fe89555..237f548 100644
--- a/TODO
+++ b/TODO
@@ -40,3 +40,5 @@
- Refactor where functions exist in .c. Lots of primality tests in factor.c.
- Switch to new text proofs.
+
+- Add ecm_factor just like the other routines
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index f2e552b..0c2bb5f 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -849,7 +849,7 @@ sub primes {
sub random_maurer_prime {
my ($n, $cert) = random_maurer_prime_with_cert(@_);
croak "maurer prime $n failed certificate verification!"
- unless verify_prime(@$cert);
+ unless verify_prime($cert);
return $n;
}
@@ -1736,7 +1736,7 @@ sub is_provable_prime {
sub prime_certificate {
my($n) = @_;
my ($is_prime, $cert) = is_provable_prime_with_cert($n);
- return @$cert;
+ return $cert;
}
@@ -1744,28 +1744,30 @@ sub is_provable_prime_with_cert {
my($n) = @_;
return 0 if defined $n && $n < 2;
_validate_num($n) || _validate_positive_integer($n);
+ my $header = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
- # Set to 0 if you want the proof to go down to 11.
- if (1) {
- if (ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL) {
- my $isp = _XS_is_prime("$n");
- return ($isp == 2) ? ($isp, [$n]) : ($isp, []);
- }
- if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime_with_cert) {
- return Math::Prime::Util::GMP::is_provable_prime_with_cert($n);
+ if (ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL) {
+ my $isp = _XS_is_prime("$n");
+ return ($isp, '') unless $isp == 2;
+ return (2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\nType Small\nN $n");
+ }
+
+ if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime_with_cert) {
+ my ($isp, $cert) = Math::Prime::Util::GMP::is_provable_prime_with_cert($n);
+ # New version that returns string format.
+ return ($isp, $cert) if ref($cert) ne 'ARRAY';
+ # Old version. Convert.
+ if (!defined $Math::Prime::Util::PrimalityProving::VERSION) {
+ eval { require Math::Prime::Util::PrimalityProving; 1; }
+ or do { croak "Cannot load Math::Prime::Util::PrimalityProving"; };
}
+ return ($isp, Math::Prime::Util::PrimalityProving::convert_array_cert_to_string($cert));
+ }
+ {
my $isp = is_prob_prime($n);
- if ($isp != 1) {
- return ($isp == 2) ? ($isp, [$n]) : ($isp, []);
- }
- } else {
- if ($n <= 10) {
- if ($n==2||$n==3||$n==5||$n==7) {
- return (2, [$n]);
- }
- return 0;
- }
+ return ($isp, '') if $isp == 0;
+ return (2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\nType Small\nN $n") if $isp == 2;
}
# Choice of methods for proof:
@@ -1781,361 +1783,39 @@ sub is_provable_prime_with_cert {
# AKS horribly slow
# See http://primes.utm.edu/prove/merged.html or other sources.
- #my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_lucas($n);
- my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_bls75($n);
+ if (!defined $Math::Prime::Util::PrimalityProving::VERSION) {
+ eval { require Math::Prime::Util::PrimalityProving; 1; }
+ or do { croak "Cannot load Math::Prime::Util::PrimalityProving"; };
+ }
+
+ #my ($isp, $pref) = Math::Prime::Util::PrimalityProving::primality_proof_lucas($n);
+ my ($isp, $pref) = Math::Prime::Util::PrimalityProving::primality_proof_bls75($n);
carp "proved $n is not prime\n" if !$isp;
return ($isp, $pref);
}
sub verify_prime {
- my @pdata = @_;
- my $verbose = $_Config{'verbose'};
-
- # Handle case of being handed a reference to the certificate.
- @pdata = @{$pdata[0]} if scalar @pdata == 1 && ref($pdata[0]) eq 'ARRAY';
-
- # Empty input = no certificate = not prime
- return 0 if scalar @pdata == 0;
-
- my $n = shift @pdata;
- if (length($n) == 1) {
- return 1 if $n =~ /^[2357]$/;
- print "primality fail: $n tiny and not prime\n" if $verbose;
- return 0;
- }
-
- 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") if ref($n) ne 'Math::BigInt';
- if ($n->is_even) {
- print "primality fail: $n even\n" if $verbose;
- return 0;
- }
-
- my $method = (scalar @pdata > 0) ? shift @pdata : 'BPSW';
-
- if ($method eq 'BPSW') {
- if ($n > Math::BigInt->new("18446744073709551615")) {
- print "primality fail: $n too large for BPSW proof\n" if $verbose;
- return 0;
- }
- my $bpsw = 0;
- my $intn = int($n->bstr);
- if ($n->bcmp("$intn") == 0 && $intn <= $_XS_MAXVAL) {
- $bpsw = _XS_miller_rabin($intn, 2)
- && _XS_is_extra_strong_lucas_pseudoprime($intn);
- } elsif ($_HAVE_GMP) {
- $bpsw = Math::Prime::Util::GMP::is_prob_prime($n);
- } else {
- $bpsw = Math::Prime::Util::PP::miller_rabin($n, 2)
- && Math::Prime::Util::PP::is_strong_lucas_pseudoprime($n);
- }
- if (!$bpsw) {
- print "primality fail: BPSW indicated $n is composite\n" if $verbose;
- return 0;
- }
- print "primality success: $n by BPSW\n" if $verbose > 1;
- return 1;
- }
-
- if ($method eq 'Pratt' || $method eq 'Lucas') {
- # Based on Lucas primality test, which requires full n-1 factorization.
- if (scalar @pdata != 2 || (ref($pdata[0]) ne 'ARRAY') || (ref($pdata[1]) eq 'ARRAY')) {
- carp "verify_prime: incorrect Pratt format, must have factors and a value\n";
- return 0;
- }
- my @factors = @{shift @pdata};
- my $a = Math::BigInt->new(shift @pdata);
- my $nm1 = $n - 1;
- my $B = $nm1; # Unfactored part
-
- my @prime_factors;
- my %factors_seen;
- foreach my $farray (@factors) {
- my $f;
- if (ref($farray) eq 'ARRAY') {
- $f = Math::BigInt->new("$farray->[0]");
- return 0 unless verify_prime(@$farray);
- } else {
- $f = $farray;
- return 0 unless verify_prime($f);
- }
- next if defined $factors_seen{"$f"}; # repeated factors
- if (($B % $f) != 0) {
- print "primality fail: given factor $f does not divide $nm1\n" if $verbose;
- return 0;
- }
- while (($B % $f) == 0) {
- $B /= $f;
- }
- push @prime_factors, $f;
- $factors_seen{"$f"} = 1;
- }
- if ($B != 1) {
- print "primality fail: n-1 not completely factored" if $verbose;
- return 0;
- }
-
- # 1. a must be co-prime to n.
- if (Math::BigInt::bgcd($a, $n) != 1) {
- print "primality fail: a and n not coprime\n" if $verbose;
- return 0;
- }
- # 2. n is a psp base a
- if ($a->copy->bmodpow($nm1, $n) != 1) {
- print "primality fail: n is not a psp base a\n" if $verbose;
- return 0;
- }
- # 3. For each factor f of n-1, a^((n-1)/f) != 1 mod n
- foreach my $f (@prime_factors) {
- if ($a->copy->bmodpow(int($nm1/$f),$n) == 1) {
- print "primality fail: factor f fails a^((n-1)/f) != 1 mod n\n" if $verbose;
- return 0;
- }
- }
- print "primality success: $n by Lucas test\n" if $verbose > 1;
- return 1;
- }
-
- if ($method eq 'n-1') {
- # BLS75 or generalized Pocklington
- # http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
- # Pull off optional theorem 7 data
- my ($theorem, $t7_B1, $t7_B, $t7_a) = (5, undef, undef, undef);
- if (scalar @pdata == 3 && ref($pdata[0]) eq 'ARRAY' && $pdata[0]->[0] =~ /^(B|T7|Theorem\s*7)$/i) {
- $theorem = 7;
- (undef, $t7_B1, $t7_B, $t7_a) = @{shift @pdata};
- $t7_B = Math::BigInt->new("$t7_B");
- }
- if (scalar @pdata != 2 || (ref($pdata[0]) ne 'ARRAY') || (ref($pdata[1]) ne 'ARRAY')) {
- carp "verify_prime: incorrect n-1 format, must have factors and a values\n";
- return 0;
- }
- my @factors = @{shift @pdata};
- my @as = @{shift @pdata};
- if ($#factors != $#as) {
- carp "verify_prime: incorrect n-1 format, must have a value for each factor\n";
- return 0;
- }
-
- my $nm1 = $n - 1;
- my $A = $n-$n+1; # Factored part (F_1 in BLS paper)
- my $B = $nm1; # Unfactored part (R_1 in BLS paper)
-
- my @prime_factors;
- my @pfas;
- my %factors_seen;
- foreach my $farray (@factors) {
- my $f;
- my $a = shift @as;
- if (ref($farray) eq 'ARRAY') {
- $f = Math::BigInt->new("$farray->[0]");
- return 0 unless verify_prime(@$farray);
- } else {
- $f = Math::BigInt->new("$farray");
- return 0 unless verify_prime($f);
- }
- next if defined $factors_seen{"$f"}; # repeated factors
- if (($B % $f) != 0) {
- print "primality fail: given factor $f does not divide $nm1\n" if $verbose;
- return 0;
- }
- while (($B % $f) == 0) {
- $B /= $f;
- $A *= $f;
- }
- push @prime_factors, $f;
- push @pfas, $a;
- $factors_seen{"$f"} = 1;
- }
- croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
-
- # The theorems state that A is the even portion, so we are requiring 2 be
- # listed as a factor.
- if ($A->is_odd) {
- print "primality fail: 2 must be included as a factor" if $verbose;
- return 0;
- }
+ my @cdata = @_;
- # TODO: consider: if B=1 and a single a is given, then Lucas test.
-
- if (Math::BigInt::bgcd($A, $B) != 1) {
- print "primality fail: A and B not coprime\n" if $verbose;
- return 0;
- }
- if ($theorem == 7) {
- if ($B != $t7_B) {
- print "primality fail: T7 unfactored != unfactored\n" if $verbose;
- return 0;
- }
- if ($t7_B1 < 1) {
- print "primality fail: T7 B value < 1\n" if $verbose;
- return 0;
- }
- # 1. Check $B has no factors smaller than $t7_B1
- my $no_small_factors = 0;
- if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::trial_factor) {
- my @trial = Math::Prime::Util::GMP::trial_factor($B, $t7_B1);
- $no_small_factors = (scalar @trial == 1);
- } elsif ($B <= ''.~0) {
- my @trial = Math::Prime::Util::PP::trial_factor($B, $t7_B1);
- $no_small_factors = (scalar @trial == 1);
- } else {
- # This is slow when B1 > 1M, but with a bigint B it's faster than
- # doing trial divisions (but much slower with native B).
- $no_small_factors =
- (Math::BigInt::bgcd($B, primorial(Math::BigInt->new($t7_B1))) == 1);
- }
- if (!$no_small_factors) {
- print "primality fail: T7 B has a factor smaller than B1\n" if $verbose;
- return 0;
- }
- # 2. Add $B and $t7_a to lists so they get checked as part of (I).
- push @prime_factors, $B;
- push @pfas, $t7_a;
- }
- { # Theorem 5, m = 1, page 624; Theorem 7, page 626
- my ($s,$r) = $B->copy->bdiv($A->copy->bmul(2));
- my $fpart = ($theorem == 7)
- ? ($A*$t7_B1+1) * (2*$A*$A + ($r-$t7_B1) * $A + 1)
- : ($A+1) * (2*$A*$A + ($r-1) * $A + 1);
- if ($n >= $fpart) {
- print "primality fail: not enough factors\n" if $verbose;
- return 0;
- }
- my $rtest = $r*$r - 8*$s;
- my $rtestroot = $rtest->copy->bsqrt;
- if ($s != 0 && ($rtestroot*$rtestroot) == $rtest) {
- print "primality fail: BLS75 theorem 5: s=$s, r=$r indicates composite\n" if $verbose;
- return 0;
- }
- }
- # Now verify (I), page 623
- foreach my $i (0 .. $#prime_factors) {
- my $f = $prime_factors[$i];
- my $a = Math::BigInt->new("$pfas[$i]");
- if ($a->copy->bmodpow($nm1, $n) != 1 ||
- Math::BigInt::bgcd($a->copy->bmodpow($nm1/$f, $n)->bsub(1), $n) != 1) {
- print "primality fail: BLS75 factor=$f, a=$a failed.\n" if $verbose;
- return 0;
- }
- }
- print "primality success: $n by BLS75 theorem $theorem\n" if $verbose > 1;
- return 1;
+ if (!defined $Math::Prime::Util::PrimalityProving::VERSION) {
+ eval { require Math::Prime::Util::PrimalityProving; 1; }
+ or do { croak "Cannot load Math::Prime::Util::PrimalityProving"; };
}
- if ($method eq 'ECPP' || $method eq 'AGKM') {
- # EC cert: Atkin-Morain etc.
- # Normally we'd have the q values set up recursively, but to follow the
- # standard trend, we have this set up as a list:
- # n, "AGKM", [n,a,b,m,q,P], [n1,a,b,m,q,P], ...
- #
- # Examples:
- # (100000000000000000039, "AGKM", [100000000000000000039, 31484432173069852672, 39553474583282556928, 100000000014867206541, 539348143913549, [39164891430400385024,86449249723524901718]])
- # (677826928624294778921, "AGKM", [677826928624294778921, 404277700094248015180, 599134911995823048257, 677826928656744857936, 104088901820753203, [2293544533, 356794037129589115041]])
- # Ux,Uy should be 600992528322000913770, 206075883056439332684
- # Vx,Vy should be 0, 1
- if (scalar @pdata < 1) {
- carp "verify_prime: incorrect AGKM format\n";
+ my $cert = '';
+ if (scalar @cdata == 1 && ref($cdata[0]) eq '') {
+ $cert = $cdata[0];
+ } else {
+ # We've been given an old array cert
+ $cert = Math::Prime::Util::PrimalityProving::convert_array_cert_to_string(@cdata);
+ if ($cert eq '') {
+ print "primality fail: error converting old certificate" if $_Config{'verbose'};
return 0;
}
- my ($ni, $a, $b, $m, $P);
- my ($qval, $q) = ($n, $n);
- foreach my $block (@pdata) {
- if (ref($block) ne 'ARRAY' || scalar @$block != 6) {
- carp "verify_prime: incorrect AGKM block format\n";
- return 0;
- }
- if (Math::BigInt->new("$block->[0]") != Math::BigInt->new("$q")) {
- carp "verify_prime: incorrect AGKM block format: block n != q\n";
- return 0;
- }
- ($ni, $a, $b, $m, $qval, $P) = @$block;
- $q = ref($qval) eq 'ARRAY' ? $qval->[0] : $qval;
- if (ref($P) ne 'ARRAY' || scalar @$P != 2) {
- carp "verify_prime: incorrect AGKM block point format\n";
- return 0;
- }
- my($Px, $Py) = @$P;
- $ni = $n->copy->bzero->badd("$ni") unless ref($ni) eq 'Math::BigInt';
- $a = $n->copy->bzero->badd("$a") unless ref($a) eq 'Math::BigInt';
- $b = $n->copy->bzero->badd("$b") unless ref($b) eq 'Math::BigInt';
- $m = $n->copy->bzero->badd("$m") unless ref($m) eq 'Math::BigInt';
- $q = $n->copy->bzero->badd("$q") unless ref($q) eq 'Math::BigInt';
- $Px = $n->copy->bzero->badd("$Px") unless ref($Px) eq 'Math::BigInt';
- $Py = $n->copy->bzero->badd("$Py") unless ref($Py) eq 'Math::BigInt';
- if ( $ni <= 0 ) {
- print "primality fail: AGKM block n is 0 or negative\n" if $verbose;
- return 0;
- }
- if (Math::BigInt::bgcd($ni, 6) != 1) {
- print "primality fail: AGKM block n '$ni' is divisible by 2 or 3\n" if $verbose;
- return 0;
- }
- my $c = $a*$a*$a * 4 + $b*$b * 27;
- if (Math::BigInt::bgcd($c, $ni) != 1) {
- print "primality fail: AGKM block gcd 4a^3+27b^2,n incorrect\n" if $verbose;
- return 0;
- }
- if ( ($Py*$Py % $ni) != (($Px*$Px*$Px + $a*$Px + $b) % $ni) ) {
- print "primality fail: AGKM block y^2 != x^3 + ax + b\n" if $verbose;
- return 0;
- }
- if ( $m < ($ni - 2*$ni->copy->bsqrt + 1)) {
- print "primality fail: AGKM block m too small\n" if $verbose;
- return 0;
- }
- if ( $m > ($ni + 2*$ni->copy->bsqrt + 1)) {
- print "primality fail: AGKM block m too large\n" if $verbose;
- return 0;
- }
- if ( $q > $ni || $q <= 0 ) {
- print "primality fail: AGKM block q invalid\n" if $verbose;
- return 0;
- }
- if ( ($m == $q) || ($m % $q) != 0 ) {
- print "primality fail: AGKM block m is not a multiple of q\n" if $verbose;
- return 0;
- }
- if ($q <= $ni->copy->broot(4)->badd(1)->bpow(2)) {
- print "primality fail: AGKM block q is too small\n" if $verbose;
- return 0;
- }
- # Final check, check that we've got a bound above and below (Hasse)
- my $correct_point = 0;
- if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::_validate_ecpp_curve) {
- $correct_point = Math::Prime::Util::GMP::_validate_ecpp_curve($a, $b, $ni, $Px, $Py, $m, $q);
- } else {
- if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) {
- eval { require Math::Prime::Util::ECAffinePoint; 1; }
- or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; };
- }
- my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $ni, $Px, $Py);
- # Compute U = (m/q)P, check U != point at infinity
- $ECP->mul( $m->copy->bdiv($q)->as_int );
- if (!$ECP->is_infinity) {
- # Compute V = qU, check V = point at infinity
- $ECP->mul( $q );
- $correct_point = 1 if $ECP->is_infinity;
- }
- }
- if (!$correct_point) {
- print "primality fail: AGKM point does not multiply correctly.\n" if $verbose;
- return 0;
- }
- }
- # Check primality of last q
- return 0 unless verify_prime($qval);
-
- print "primality success: $n by A-K-G-M elliptic curve\n" if $verbose > 1;
- return 1;
}
-
- carp "verify_prime: Unknown method: '$method'.\n";
- return 0;
+ return 0 if $cert eq '';
+ return Math::Prime::Util::PrimalityProving::verify_cert($cert);
}
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 00e2f2a..faa9f72 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -2070,190 +2070,6 @@ sub ecm_factor {
}
-sub primality_proof_lucas {
- my ($n) = shift;
- my @composite = (0, []);
-
- # Since this can take a very long time with a composite, try some easy cuts
- return @composite if !defined $n || $n < 2;
- return (2, [$n]) if $n < 4;
- return @composite if miller_rabin($n,2,15,325) == 0;
-
- if (!defined $Math::BigInt::VERSION) {
- eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
- or do { croak "Cannot load Math::BigInt"; };
- }
-
- my $nm1 = $n-1;
- my @factors = factor($nm1);
- { # remove duplicate factors and make a sorted array of bigints
- my %uf;
- undef @uf{@factors};
- @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
- }
- for (my $a = 2; $a < $nm1; $a++) {
- my $ap = Math::BigInt->new("$a");
- # 1. a must be coprime to n
- next unless Math::BigInt::bgcd($ap, $n) == 1;
- # 2. a^(n-1) = 1 mod n.
- next unless $ap->copy->bmodpow($nm1, $n) == 1;
- # 3. a^((n-1)/f) != 1 mod n for all f.
- next if (scalar grep { $_ == 1 }
- map { $ap->copy->bmodpow(int($nm1/$_),$n); }
- @factors) > 0;
- # Verify each factor and add to proof
- my @fac_proofs;
- foreach my $f (@factors) {
- my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
- if ($isp != 2) {
- carp "could not prove primality of $n.\n";
- return (1, []);
- }
- push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
- }
- my @proof = ("$n", "Pratt", [@fac_proofs], $a);
- return (2, [@proof]);
- }
- return @composite;
-}
-
-sub primality_proof_bls75 {
- my ($n) = shift;
- my @composite = (0, []);
-
- # Since this can take a very long time with a composite, try some easy cuts
- return @composite if !defined $n || $n < 2;
- return (2, [$n]) if $n < 4;
- return @composite if ($n & 1) == 0;
- return @composite if miller_rabin($n,2,15,325) == 0;
-
- 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") unless ref($n) eq 'Math::BigInt';
- my $nm1 = $n-1;
- my $A = $nm1->copy->bone; # factored part
- my $B = $nm1->copy; # unfactored part
- my @factors = (2);
- croak "BLS75 error: n-1 not even" unless $nm1->is_even();
- my $trial_B = 20000;
- {
- while ($B->is_even) { $B /= 2; $A *= 2; }
- my @tf = trial_factor($B, $trial_B);
- pop @tf if $tf[-1] > $trial_B;
- foreach my $f (@tf) {
- next if $f == $factors[-1];
- push @factors, $f;
- do { $B /= $f; $A *= $f; } while (($B % $f) == 0);
- }
- }
- my @nstack;
- # nstack should only hold composites
- if ($B == 1) {
- # Completely factored. Nothing.
- } elsif (Math::Prime::Util::is_prob_prime($B)) {
- push @factors, $B;
- $A *= $B; $B /= $B; # completely factored already
- } else {
- push @nstack, $B;
- }
- my $theorem = 5;
- 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....
- $fpart = ($A*$trial_B+1) * (2*$A*$A + ($r-$trial_B) * $A + 1);
- if ($n < $fpart) { $theorem = 7; last; }
-
- my $m = pop @nstack;
- # Don't use bignum if it has gotten small enough.
- $m = int($m->bstr) if ref($m) eq 'Math::BigInt' && $m <= ''.~0;
- # Try to find factors of m, using the default set of factor subs.
- my @ftry;
- $_holf_r = 1;
- foreach my $sub (@_fsublist) {
- @ftry = $sub->($m);
- last if scalar @ftry >= 2;
- }
- # If we couldn't find a factor, skip it.
- next unless scalar @ftry > 1;
- # Process each factor
- foreach my $f (@ftry) {
- croak "Invalid factoring: B=$B m=$m f=$f" if $f == 1 || $f == $m || ($B%$f) != 0;
- if (Math::Prime::Util::is_prob_prime($f)) {
- push @factors, $f;
- do { $B /= $f; $A *= $f; } while (($B % $f) == 0);
- } else {
- push @nstack, $f;
- }
- }
- }
- # Just in case:
- foreach my $f (@factors) {
- while (($B % $f) == 0) {
- $B /= $f; $A *= $f;
- }
- }
- { # remove duplicate factors and make a sorted array of bigints
- my %uf = map { $_ => 1 } @factors;
- @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
- }
- # Did we factor enough?
- my ($s,$r) = $B->copy->bdiv($A->copy->bmul(2));
- my $fpart = ($theorem == 5)
- ? ($A+1) * (2*$A*$A + ($r-1) * $A + 1)
- : ($A*$trial_B+1) * (2*$A*$A + ($r-$trial_B) * $A + 1);
- return (1,[]) if $n >= $fpart;
- # Check we didn't mess up
- croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
- croak "BLS75 error: $A not even" unless $A->is_even();
- croak "BLS75 error: A and B not coprime" unless Math::BigInt::bgcd($A, $B)==1;
-
- my $rtest = $r*$r - 8*$s;
- my $rtestroot = $rtest->copy->bsqrt;
- return @composite if $s != 0 && ($rtestroot*$rtestroot) == $rtest;
-
- my @fac_proofs;
- my @as;
- foreach my $f (@factors) {
- my $success = 0;
- foreach my $a (2 .. 10000) {
- my $ap = Math::BigInt->new($a);
- next unless $ap->copy->bmodpow($nm1, $n) == 1;
- next unless Math::BigInt::bgcd($ap->copy->bmodpow($nm1/$f, $n)->bsub(1), $n) == 1;
- push @as, $a;
- $success = 1;
- last;
- }
- return @composite unless $success;
- my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
- if ($isp != 2) {
- carp "could not prove primality of $n.\n";
- return (1, []);
- }
- push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
- }
- # Put n, B back to non-bigints if possible.
- $n = int($n->bstr) if ref($n) eq 'Math::BigInt' && $n <= ''.~0;
- $B = int($B->bstr) if ref($B) eq 'Math::BigInt' && $B <= ''.~0;
- if ($theorem == 5) {
- return (2, [$n, "n-1", [@fac_proofs], [@as]]);
- } else {
- my $t7a = 0;
- my $f = $B;
- foreach my $a (2 .. 10000) {
- my $ap = Math::BigInt->new($a);
- next unless $ap->copy->bmodpow($nm1, $n) == 1;
- next unless Math::BigInt::bgcd($ap->copy->bmodpow($nm1/$f, $n)->bsub(1), $n) == 1;
- return (2, [$n, "n-1", ["B", $trial_B, $B, $a], [@fac_proofs], [@as]]);
- }
- }
- return @composite;
-}
-
my $_const_euler = 0.57721566490153286060651209008240243104215933593992;
my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
@@ -2936,20 +2752,6 @@ Takes a positive number as input, and returns 1 if the input can be proven
prime using the AKS primality test. This code is included for completeness
and as an example, but is impractically slow.
-=head2 primality_proof_lucas
-
-Given a positive number C<n> as input, performs a full factorization of C<n-1>,
-then attempts a Lucas test on the result. A Pratt-style certificate is
-returned. Note that if the input is composite, this will take a B<very> long
-time to return.
-
-=head2 primality_proof_bls75
-
-Given a positive number C<n> as input, performs a partial factorization of
-C<n-1>, then attempts a proof using theorem 5 of Brillhart, Lehmer, and
-Selfridge's 1975 paper. This can take a long time to return if given a
-composite, though it should not be anywhere near as long as the Lucas test.
-
=head2 lucas_sequence
my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k)
diff --git a/t/23-primality-proofs.t b/t/23-primality-proofs.t
index 66a8aed..2195fea 100644
--- a/t/23-primality-proofs.t
+++ b/t/23-primality-proofs.t
@@ -36,10 +36,10 @@ plan tests => 0
+ 2 # is_provable_prime
+ 6 * scalar(@plist)
+ 6 # hand-done proofs
- + 28 # borked up certificates generating warnings
+ + 24 # borked up certificates generating warnings
+ 6 # verification failures (tiny/BPSW)
+ 8 # verification failures (Lucas/Pratt)
- + 12 # verification failures (n-1)
+ + 11 # verification failures (n-1)
+ 7 # verification failures (ECPP)
+ 0;
@@ -54,16 +54,16 @@ foreach my $p (@plist) {
if $broken64 && $p > 2**48;
skip "These take a long time on non-64-bit. Skipping", 5
if !$use64 && !$extra && $p =~ /^(6778|9800)/;
- my($isp, $cert_ref) = is_provable_prime_with_cert($p);
+ my($isp, $cert) = is_provable_prime_with_cert($p);
is( $isp, 2, " is_provable_prime_with_cert returns 2" );
- ok( defined($cert_ref) && ref($cert_ref) eq 'ARRAY' && scalar(@$cert_ref) >= 1,
+ ok( defined($cert) && $cert =~ /^Type/m,
" certificate is non-null" );
- ok( verify_prime($cert_ref), " verification of certificate reference done" );
+ ok( verify_prime($cert), " verification of certificate done" );
# Note, in some cases the two certs could be non-equal (but both must be valid!)
- my @cert = prime_certificate($p);
- ok( scalar(@cert) >= 1, " prime_certificate is also non-null" );
+ my $cert2 = prime_certificate($p);
+ ok( defined($cert2) && $cert2 =~ /^Type/m, " prime_certificate is also non-null" );
# TODO: compare certificates and skip if equal
- ok( verify_prime(@cert), " verification of prime_certificate done" );
+ ok( verify_prime($cert2), " verification of prime_certificate done" );
}
}
@@ -88,10 +88,10 @@ SKIP: {
[ 3,5,3,2,3,3,3,3 ] );
ok( verify_prime(@proof), "2**607-1 primality proof verified" );
}
-{
- my @proof = ('809120722675364249', "n-1", ["B", 20000, '22233477760919', 2], [2, 4549], [3, 2]);
- ok( verify_prime(@proof), "n-1 T7 primality proof of 809120722675364249 verified" );
-}
+#{
+# my @proof = ('809120722675364249', "n-1", ["B", 20000, '22233477760919', 2], [2, 4549], [3, 2]);
+# ok( verify_prime(@proof), "n-1 T7 primality proof of 809120722675364249 verified" );
+#}
{
my @proof = (20907001, "Pratt", [ 2,
3,
@@ -207,10 +207,10 @@ is( verify_prime([1490266103, 'n-1', [13, 19, 1597, 1889], [2, 2, 2, 2]]), 0, "n
is( verify_prime([1490266103, 'n-1', [2, 13, 1889, 30343], [5, 2, 2, 2]]), 0, "n-1 with a non-prime factor" );
is( verify_prime([1490266103, 'n-1', [2, 13, 1889, [30343]], [5, 2, 2, 2]]), 0, "n-1 with a non-prime array factor" );
# I don't know how to make F and R (A and B) to not be coprime
-is( verify_prime(['9848131514359', 'n-1', ["B", 20000, 890588851, 2], [2, 3, 19, 97], [3, 5, 2, 2]]), 1, "n-1 T7 proper" );
-is( verify_prime(['9848131514359', 'n-1', ["B", 20000, 890588951, 2], [2, 3, 19, 97], [3, 5, 2, 2]]), 0, "n-1 T7 with misfactor" );
-is( verify_prime(['9848131514359', 'n-1', ["B", 0, 890588851, 2], [2, 3, 19, 97], [3, 5, 2, 2]]), 0, "n-1 T7 with B < 1" );
-is( verify_prime(['9848131514359', 'n-1', ["B", 20000, 16921188169, 2], [2, 3, 97], [3, 5, 2]]), 0, "n-1 T7 with wrong B" );
+#is( verify_prime(['9848131514359', 'n-1', ["B", 20000, 890588851, 2], [2, 3, 19, 97], [3, 5, 2, 2]]), 1, "n-1 T7 proper" );
+#is( verify_prime(['9848131514359', 'n-1', ["B", 20000, 890588951, 2], [2, 3, 19, 97], [3, 5, 2, 2]]), 0, "n-1 T7 with misfactor" );
+#is( verify_prime(['9848131514359', 'n-1', ["B", 0, 890588851, 2], [2, 3, 19, 97], [3, 5, 2, 2]]), 0, "n-1 T7 with B < 1" );
+#is( verify_prime(['9848131514359', 'n-1', ["B", 20000, 16921188169, 2], [2, 3, 97], [3, 5, 2]]), 0, "n-1 T7 with wrong B" );
is( verify_prime([1490266103, 'n-1', [2, 13], [5, 2]]), 0, "n-1 without enough factors" );
is( verify_prime([914144252447488195, 'n-1', [2, 3, 11, 17, 1531], [2, 2, 2, 2, 2]]), 0, "n-1 with bad BLS75 r/s" );
is( verify_prime([1490266103, 'n-1', [2, 13, 19, 1597, 1889], [3, 2, 2, 2, 2]]), 0, "n-1 with bad a value" );
diff --git a/t/80-pp.t b/t/80-pp.t
index d5a69bb..6c571e3 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -239,7 +239,7 @@ my %ipp = (
36010359 => 0,
);
-plan tests => 1 +
+plan tests => 2 +
3 +
3 + scalar(keys %small_single) + scalar(keys %small_range) +
2*scalar(keys %primegaps) + 8 + 1 + 1 + 1 +
@@ -266,6 +266,8 @@ use Math::Prime::Util qw/primes prime_count_approx prime_count_lower
use Math::BigInt try => 'GMP';
use Math::BigFloat;
require_ok 'Math::Prime::Util::PP';
+require_ok 'Math::Prime::Util::PrimalityProving';
+
# This function skips some setup
undef *primes;
*primes = \&Math::Prime::Util::PP::primes;
@@ -563,7 +565,7 @@ SKIP: {
is( is_aks_prime(74513), 0, "AKS: 74513 is composite (failed anr test)" );
}
-is_deeply( [Math::Prime::Util::PP::primality_proof_lucas(100003)],
+is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_lucas(100003)],
[2, [100003, "Pratt", [2, 3, 7, 2381], 2]],
"primality_proof_lucas(100003)" );
# Had to reduce these to make borked up Perl 5.6.2 work.
@@ -574,10 +576,10 @@ is_deeply( [Math::Prime::Util::PP::primality_proof_lucas(100003)],
# [2, ["809120722675364249", "n-1",
# ["B", 20000, "22233477760919", 2], [2, 4549], [3, 2]]],
# "primality_proof_bls75(809120722675364249)" );
-is_deeply( [Math::Prime::Util::PP::primality_proof_bls75(1490266103)],
+is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_bls75(1490266103)],
[2, [1490266103, 'n-1', [2, 13, 19, 1597, 1889], [5, 2, 2, 2, 2]]],
"primality_proof_bls75(1490266103)" );
-is_deeply( [Math::Prime::Util::PP::primality_proof_bls75(Math::BigInt->new("64020974585221"))],
+is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_bls75(Math::BigInt->new("64020974585221"))],
[2, ["64020974585221", "n-1",
["B", 20000, "5154667841", 2], [2, 3, 5, 23], [2, 2, 2, 2]]],
"primality_proof_bls75(64020974585221)" );
diff --git a/xt/primality-proofs.pl b/xt/primality-proofs.pl
index 4c0cf7b..8d32623 100755
--- a/xt/primality-proofs.pl
+++ b/xt/primality-proofs.pl
@@ -66,10 +66,23 @@ print "\n";
sub proof_mark {
my $cert = shift;
- my $type = (scalar @$cert == 1) ? "bpsw" : $cert->[1];
+ my $type;
+ if (ref($cert) eq 'ARRAY') {
+ $type = (scalar @$cert == 1) ? "bpsw" : $cert->[1];
+ if ($type =~ /n-1/i) {
+ $type = ($cert->[2]->[0] eq 'B') ? 'BLS7' : 'BLS5';
+ }
+ } else {
+ ($type) = $cert =~ /Type (\S+)/;
+ $type = 'ECPP' if $cert =~ /Type\s+ECPP/;
+ }
if (!defined $type) { die "\nNo cert:\n\n", dump_filtered($cert, $bifilter); }
- if ($type =~ /n-1/i) { return ($cert->[2]->[0] eq 'B') ? '7' : '5'; }
+ if ($type =~ /bls5/i) { return '5'; }
+ elsif ($type =~ /bls7/i) { return '7'; }
+ if ($type =~ /bls3/i) { return '-'; }
+ elsif ($type =~ /bls15/i) { return '+'; }
elsif ($type =~ /bpsw/i) { return '.'; }
elsif ($type =~ /ecpp|agkm/i) { return 'E'; }
+ warn "type: $type\n";
return '?';
}
--
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