[libmath-prime-util-perl] 37/40: Switch to new text proofs, part 2
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:06 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 f282fde046ff4f25a100d67f93f1deaaa0b0055c
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Aug 2 17:45:02 2013 -0700
Switch to new text proofs, part 2
---
Changes | 54 +++++----
TODO | 4 +-
factor.c | 1 +
lib/Math/Prime/Util.pm | 137 ++++++++++++----------
lib/Math/Prime/Util/PrimalityProving.pm | 200 +++++++++++++++++++++-----------
t/23-primality-proofs.t | 14 ++-
t/80-pp.t | 14 +--
7 files changed, 257 insertions(+), 167 deletions(-)
diff --git a/Changes b/Changes
index d16a4f4..ea29786 100644
--- a/Changes
+++ b/Changes
@@ -1,32 +1,40 @@
Revision history for Perl module Math::Prime::Util
0.30
- [Functions Added]
- - is_frobenius_underwood_pseudoprime
- - is_almost_extra_strong_lucas_pseudoprime
- - lucas_sequence
- - pplus1_factor
-
- - Documentation and PP is_prime changed to use extra strong Lucas test
- from the strong test. This matches what the newest MPU::GMP does.
- This has no effect at all for numbers < 2^64. No counter-example is
- known for the standard, strong, extra strong, or almost extra strong
- (increment 1 or 2) tests. The extra strong test is faster than the
- strong test and produces fewer pseudoprimes. It retains the residue
- class properties of the strong Lucas test (where the SPSP-2
- pseudoprimes favor residue class 1 and the Lucas pseudoprimes favor
- residue class -1), hence should retain the BPSW test strength.
-
- - XS code for all 4 Lucas tests.
- - Clean up is_prob_prime, also ~10% faster for n >= 885594169.
+ [API Changes]
- - Fixed a rare refcount / bignum / callback issue in next_prime.
+ - Primality proofs now use the new "MPU Certificate" format, which is
+ text rather than a nested Perl data structure. This is much better
+ for external interaction, especially with non-Perl tools. It is
+ not quite as convenient for all-Perl manipulation.
- - Small mulmod speedup for non-gcc/x86_64 platforms, and for any platform
- with gcc 4.4 or newer.
-
- - Add more conditions to ECPP block verification.
+ [Functions Added]
+ - is_frobenius_underwood_pseudoprime
+ - is_almost_extra_strong_lucas_pseudoprime
+ - lucas_sequence
+ - pplus1_factor
+
+ [Enhancements]
+ - Documentation and PP is_prime changed to use extra strong Lucas test
+ from the strong test. This matches what the newest MPU::GMP does.
+ This has no effect at all for numbers < 2^64. No counter-example is
+ known for the standard, strong, extra strong, or almost extra strong
+ (increment 1 or 2) tests. The extra strong test is faster than the
+ strong test and produces fewer pseudoprimes. It retains the residue
+ class properties of the strong Lucas test (where the SPSP-2
+ pseudoprimes favor residue class 1 and the Lucas pseudoprimes favor
+ residue class -1), hence should retain the BPSW test strength.
+
+ - XS code for all 4 Lucas tests.
+
+ - Clean up is_prob_prime, also ~10% faster for n >= 885594169.
+
+ - Small mulmod speedup for non-gcc/x86_64 platforms, and for any platform
+ with gcc 4.4 or newer.
+
+ [Bug Fixes]
+ - Fixed a rare refcount / bignum / callback issue in next_prime.
0.29 2013-05-30
diff --git a/TODO b/TODO
index 237f548..7662f08 100644
--- a/TODO
+++ b/TODO
@@ -39,6 +39,8 @@
- Refactor where functions exist in .c. Lots of primality tests in factor.c.
+- Add ecm_factor just like the other routines
+
- Switch to new text proofs.
-- Add ecm_factor just like the other routines
+- Rewrite 23-primality-proofs.t for new format (keep some of the old tests?).
diff --git a/factor.c b/factor.c
index 146bfeb..0530b32 100644
--- a/factor.c
+++ b/factor.c
@@ -1030,6 +1030,7 @@ int pplus1_factor(UV n, UV *factors, UV B1)
X1 = 7 % n;
X2 = 11 % n;
+ f = 1;
START_DO_FOR_EACH_PRIME(2, B1) {
UV k = p;
if (p < sqrtB1) {
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 0c2bb5f..c7cf487 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -847,7 +847,13 @@ sub primes {
}
sub random_maurer_prime {
- my ($n, $cert) = random_maurer_prime_with_cert(@_);
+ my $k = shift;
+ _validate_num($k, 2) || _validate_positive_integer($k, 2);
+ if ($k <= $_Config{'maxbits'}) {
+ croak "Random Maurer not supported on old Perl" if $k > 49 && $] < 5.008 && $_Config{'maxbits'} > 32;
+ return random_nbit_prime($k);
+ }
+ my ($n, $cert) = random_maurer_prime_with_cert($k);
croak "maurer prime $n failed certificate verification!"
unless verify_prime($cert);
return $n;
@@ -856,11 +862,10 @@ sub primes {
sub random_maurer_prime_with_cert {
my($k) = @_;
_validate_num($k, 2) || _validate_positive_integer($k, 2);
- my @cert;
if ($] < 5.008 && $_Config{'maxbits'} > 32) {
if ($k <= 49) {
my $n = random_nbit_prime($k);
- return ($n, [$n]);
+ return ($n, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\nType Small\nN $n\n");
}
croak "Random Maurer not supported on old Perl";
}
@@ -890,10 +895,9 @@ sub primes {
my $verbose = $_Config{'verbose'};
local $| = 1 if $verbose > 2;
- my $c = Math::BigFloat->new("0.065"); # higher = more trial divisions
+ # Ignore Maurer's g and c that controls how much trial division is done.
my $r = Math::BigFloat->new("0.5"); # relative size of the prime q
my $m = 20; # makes sure R is big enough
- my $B = ($c * $k * $k)->bfloor;
my $irandf = _get_rand_func();
# Generate a random prime q of size $r*$k, where $r >= 0.5. Try to
@@ -906,10 +910,12 @@ sub primes {
}
# I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1.
- my ($q, $certref) = random_maurer_prime_with_cert( ($r * $k)->bfloor + 1 );
+ my ($q, $qcert) = random_maurer_prime_with_cert( ($r * $k)->bfloor + 1 );
$q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt';
my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int();
- print "B = $B r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3;
+ print "r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3;
+ $qcert = ($q < Math::BigInt->new("18446744073709551615"))
+ ? "" : _strip_proof_header($qcert);
# Big GCD's are hugely fast with GMP or Pari, but super slow with Calc.
if ($_big_gcd_use < 0) {
@@ -927,72 +933,47 @@ sub primes {
my $n = Math::BigInt->new(2)->bmul($R)->bmul($q)->badd(1);
# We constructed a promising looking $n. Now test it.
print "." if $verbose > 2;
-
# Trial divisions, trying to quickly weed out non-primes.
next unless Math::BigInt::bgcd($n, 111546435) == 1;
if ($_big_gcd_use && $n > $_big_gcd_top) {
next unless Math::BigInt::bgcd($n, $_big_gcd[0]) == 1;
next unless Math::BigInt::bgcd($n, $_big_gcd[1]) == 1;
- next unless Math::BigInt::bgcd($n, $_big_gcd[2]) == 1;
- next unless Math::BigInt::bgcd($n, $_big_gcd[3]) == 1;
+ if (!$_HAVE_GMP) {
+ next unless Math::BigInt::bgcd($n, $_big_gcd[2]) == 1;
+ next unless Math::BigInt::bgcd($n, $_big_gcd[3]) == 1;
+ }
}
print "+" if $verbose > 2;
- if ($_HAVE_GMP) {
- next unless Math::Prime::Util::GMP::is_strong_pseudoprime($n, 2);
- }
+ next unless Math::Prime::Util::is_strong_pseudoprime($n, 3);
print "*" if $verbose > 2;
- # Now we do Lemma 1 -- a special case of the Pocklington test.
- # Let F = q where q is prime, and n = 2RF+1.
- # If F > sqrt(n) or F odd and F > R, and a^((n-1)/F)-1 mod n = 1, n prime.
-
- # Based on our construction, this should always be true. Check anyway.
- next unless $q > $R;
-
- # Select random 'a' values. If n is prime, then almost any 'a' value
- # will work, so just try two small ones instead of generating a giant
- # random 'a' between 2 and n-2. This makes the powmods run faster.
- foreach my $try_a (2, 7) {
- # my $a = 2 + $irandf->( $n - 4 );
- my $a = Math::BigInt->new($try_a);
- my $b = $a->copy->bmodpow($n-1, $n);
- next unless $b == 1;
-
- # Now do the one gcd check we need to do.
- $b = $a->copy->bmodpow(2*$R, $n);
- next unless Math::BigInt::bgcd($b-1, $n) == 1;
- if ($verbose) {
- print "", ($verbose == 3) ? "($k)" : "$n passed final gcd\n";
- }
+ # We could set r = 0.3333 earlier, then use BLS75 theorem 5 here.
+ # The chain would be shorter, requiring less overall work for
+ # large inputs. Maurer's paper discusses the idea.
- # Instead of the previous gcd, we could check q >= n**1/3 and also do
- # some tests on x & y from 2R = xq+y (see Lemma 2 from Maurer's paper).
- # Crypt::Primes does the q test but doesn't do the x/y tests.
- # next if ($q <= $n->copy->broot(3));
- # my $x = (2*$R)->bdiv($q)->bfloor;
- # my $y = 2*$R - $x*$q;
- # my $z = $y*$y - 4*$x;
- # next if $z == 0;
- # next if $z->bsqrt->bfloor->bpow(2) == $z; # perfect square
- # Menezes seems to imply only the q test needs to be done, but this
- # doesn't follow from Lemma 2. Also note the entire POINT of going to
- # Lemma 2 is that we now allow r to be 0.3334, making q smaller. If we
- # run this without changing r, then x will typically be 0 and this fails.
-
- # Verify with a BPSW test on the result. This could:
- # 1) save us from accidently outputting a non-prime due to some mistake
- # 2) make history by finding the first known BPSW pseudo-prime
- croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n);
+ # Use BLS75 theorem 3. This is easier and possibly faster than
+ # BLS75 theorem 4 (Pocklington) used by Maurer's paper.
- # Build up cert, knowing n-1 = 2*q*R, q > sqrt(n).
- # We'll need to find the right a value for the factor 2.
- foreach my $f2a (2 .. 200) {
- $a = Math::BigInt->new($f2a);
- next unless $a->copy->bmodpow($n-1, $n) == 1;
- next unless Math::BigInt::bgcd($a->copy->bmodpow(($n-1)/2, $n)->bsub(1), $n) == 1;
- @cert = ("$n", "n-1", [2, [@$certref]], [$f2a, $try_a]);
- return ($n, \@cert);
- }
+ # Check conditions -- these should be redundant.
+ my $m = 2 * $R;
+ if (! ($q->is_odd && $q > 2 && $m > 0 &&
+ $m * $q + 1 == $n && 2*$q+1 > $n->copy->bsqrt()) ) {
+ carp "Maurer prime failed BLS75 theorem 3 conditions. Retry.";
+ next;
+ }
+ # Find a suitable a. Move on if one isn't found quickly.
+ foreach my $trya (2, 3, 5, 7, 11, 13) {
+ my $a = Math::BigInt->new($trya);
+ # m/2 = R (n-1)/2 = (2*R*q)/2 = R*q
+ next unless $a->copy->bmodpow($R,$n) != $n-1;
+ next unless $a->copy->bmodpow($R*$q, $n) == $n-1;
+ print "($k)" if $verbose > 2;
+ croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n);
+ my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" .
+ "Proof for:\nN $n\n\n" .
+ "Type BLS3\nN $n\nQ $q\nA $a\n" .
+ $qcert;
+ return ($n, $cert);
}
# Didn't pass the selected a values. Try another R.
}
@@ -1582,6 +1563,31 @@ sub factor {
return Math::Prime::Util::PP::factor($n);
}
+#sub trial_factor {
+# my($n, $limit) = @_;
+# _validate_num($n) || _validate_positive_integer($n);
+# $limit = 2147483647 unless defined $limit;
+# _validate_num($limit) || _validate_positive_integer($limit);
+# return _XS_trial_factor($n, $limit) if $n <= $_XS_MAXVAL;
+# if ($_HAVE_GMP) {
+# my @factors;
+# while (1) {
+# last if $n <= 1 || is_prob_prime($n);
+# my @f = sort { $a <=> $b }
+# Math::Prime::Util::GMP::trial_factor($n, $limit);
+# pop(@f); # Remove remainder
+# last unless scalar @f > 0;
+# foreach my $f (@f) {
+# push @factors, $f;
+# $n /= $f;
+# }
+# }
+# push @factors, $n if $n > 1;
+# return @factors;
+# }
+# return Math::Prime::Util::PP::trial_factor($n, $limit);
+#}
+
sub is_pseudoprime {
my($n, $a) = @_;
_validate_num($n) || _validate_positive_integer($n);
@@ -1718,6 +1724,13 @@ sub is_aks_prime {
return Math::Prime::Util::PP::is_aks_prime($n);
}
+# For stripping off the header on certificates so they can be combined.
+sub _strip_proof_header {
+ my $proof = shift;
+ $proof =~ s/^\[MPU - Primality Certificate\]\nVersion \S+\n+Proof for:\nN (\d+)\n+//ms;
+ return $proof;
+}
+
# Return just the non-cert portion.
sub is_provable_prime {
my($n) = @_;
diff --git a/lib/Math/Prime/Util/PrimalityProving.pm b/lib/Math/Prime/Util/PrimalityProving.pm
index db16396..8582b49 100644
--- a/lib/Math/Prime/Util/PrimalityProving.pm
+++ b/lib/Math/Prime/Util/PrimalityProving.pm
@@ -19,6 +19,7 @@ BEGIN {
$Math::Prime::Util::PrimalityProving::VERSION = '0.30';
}
+my $_smallval = Math::BigInt->new("18446744073709551615");
###############################################################################
@@ -26,27 +27,43 @@ BEGIN {
###############################################################################
my @_fsublist = (
- sub { Math::Prime::Util::prho_factor (shift, 8*1024) },
- sub { Math::Prime::Util::pminus1_factor(shift, 10_000) },
- sub { Math::Prime::Util::pbrent_factor (shift, 32*1024) },
- sub { Math::Prime::Util::pminus1_factor(shift, 1_000_000) },
- sub { Math::Prime::Util::pbrent_factor (shift, 512*1024) },
- #sub { Math::Prime::Util::ecm_factor (shift, 1_000, 5_000, 10) },
- sub { Math::Prime::Util::pminus1_factor(shift, 4_000_000) },
- #sub { Math::Prime::Util::ecm_factor (shift, 10_000, 50_000, 10) },
- sub { Math::Prime::Util::pminus1_factor(shift,20_000_000) },
- #sub { Math::Prime::Util::ecm_factor (shift, 100_000, 800_000, 10) },
- #sub { Math::Prime::Util::ecm_factor (shift, 1_000_000, 1_000_000, 10) },
- sub { Math::Prime::Util::pminus1_factor(shift, 100_000_000, 500_000_000) },
+ sub { Math::Prime::Util::PP::prho_factor (shift, 8*1024, 3) },
+ sub { Math::Prime::Util::PP::pminus1_factor(shift, 10_000) },
+ sub { Math::Prime::Util::PP::pbrent_factor (shift, 32*1024, 1) },
+ sub { Math::Prime::Util::PP::pminus1_factor(shift, 1_000_000) },
+ sub { Math::Prime::Util::PP::pbrent_factor (shift, 512*1024, 7) },
+ sub { Math::Prime::Util::PP::ecm_factor (shift, 1_000, 5_000, 10) },
+ sub { Math::Prime::Util::PP::pminus1_factor(shift, 4_000_000) },
+ sub { Math::Prime::Util::PP::pbrent_factor (shift, 512*1024, 11) },
+ sub { Math::Prime::Util::PP::ecm_factor (shift, 10_000, 50_000, 10) },
+ sub { Math::Prime::Util::PP::pminus1_factor(shift,20_000_000) },
+ sub { Math::Prime::Util::PP::ecm_factor (shift, 100_000, 800_000, 10) },
+ sub { Math::Prime::Util::PP::pbrent_factor (shift, 2048*1024, 13) },
+ sub { Math::Prime::Util::PP::ecm_factor (shift, 1_000_000, 1_000_000, 20)},
+ sub { Math::Prime::Util::PP::pminus1_factor(shift, 100_000_000, 500_000_000)},
);
+sub _small_cert {
+ my $n = shift;
+ return '' unless is_prob_prime($n);
+ return join "\n", "[MPU - Primality Certificate]",
+ "Version 1.0",
+ "",
+ "Proof for:",
+ "N $n",
+ "",
+ "Type Small",
+ "N $n",
+ "";
+}
+
sub primality_proof_lucas {
my ($n) = shift;
- my @composite = (0, []);
+ 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 (2, _small_cert($n)) if $n < 4;
return @composite if is_strong_pseudoprime($n,2,15,325) == 0;
if (!defined $Math::BigInt::VERSION) {
@@ -61,6 +78,11 @@ sub primality_proof_lucas {
undef @uf{@factors};
@factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
}
+ my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
+ $cert .= "Type Lucas\nN $n\n";
+ foreach my $i (1 .. scalar @factors) {
+ $cert .= "Q[$i] " . $factors[$i-1] . "\n";
+ }
for (my $a = 2; $a < $nm1; $a++) {
my $ap = Math::BigInt->new("$a");
# 1. a must be coprime to n
@@ -77,23 +99,26 @@ sub primality_proof_lucas {
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, []);
+ return (1, '');
}
- push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
+ push @fac_proofs, Math::Prime::Util::_strip_proof_header($fproof) if $f > $_smallval;
+ }
+ $cert .= "A $a\n";
+ foreach my $proof (@fac_proofs) {
+ $cert .= "\n$proof";
}
- my @proof = ("$n", "Pratt", [@fac_proofs], $a);
- return (2, [@proof]);
+ return (2, $cert);
}
return @composite;
}
sub primality_proof_bls75 {
my ($n) = shift;
- my @composite = (0, []);
+ my @composite = (0, '');
# Since this can take a very long time with a composite, try some easy tests
return @composite if !defined $n || $n < 2;
- return (2, [$n]) if $n < 4;
+ return (2, _small_cert($n)) if $n < 4;
return @composite if ($n & 1) == 0;
return @composite if is_strong_pseudoprime($n,2,15,325) == 0;
@@ -106,12 +131,12 @@ sub primality_proof_bls75 {
my $trial_B = 20000;
{
while ($B->is_even) { $B /= 2; $A *= 2; }
- my @tf = Math::Prime::Util::trial_factor($B, $trial_B);
+ my @tf = Math::Prime::Util::PP::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);
+ while (($B % $f) == 0) { $B /= $f; $A *= $f; }
}
}
my @nstack;
@@ -124,14 +149,10 @@ sub primality_proof_bls75 {
} 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.
@@ -167,10 +188,8 @@ sub primality_proof_bls75 {
}
# 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;
+ my $fpart = ($A+1) * (2*$A*$A + ($r-1) * $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();
@@ -180,42 +199,41 @@ sub primality_proof_bls75 {
my $rtestroot = $rtest->copy->bsqrt;
return @composite if $s != 0 && ($rtestroot*$rtestroot) == $rtest;
+ my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
+ $cert .= "Type BLS5\nN $n\n";
+ my $qnum = 0;
+ my $atext = '';
my @fac_proofs;
- my @as;
foreach my $f (@factors) {
my $success = 0;
+ if ($qnum == 0) {
+ die "BLS5 Perl proof: Internal error, first factor not 2" unless $f == 2;
+ } else {
+ $cert .= "Q[$qnum] $f\n";
+ }
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;
+ $atext .= "A[$qnum] $a\n" unless $a == 2;
$success = 1;
last;
}
+ $qnum++;
return @composite unless $success;
my ($isp, $fproof) = is_provable_prime_with_cert($f);
if ($isp != 2) {
carp "could not prove primality of $n.\n";
- return (1, []);
+ return (1, '');
}
- push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
+ push @fac_proofs, Math::Prime::Util::_strip_proof_header($fproof) if $f > $_smallval;
}
- # 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]]);
- }
+ $cert .= $atext;
+ $cert .= "----\n";
+ foreach my $proof (@fac_proofs) {
+ $cert .= "\n$proof";
}
- return @composite;
+ return (2, $cert);
}
###############################################################################
@@ -237,15 +255,16 @@ sub _convert_cert {
my $method = (scalar @$pdata > 0) ? shift @$pdata : 'BPSW';
if ($method eq 'BPSW') {
- return '' if $n > Math::BigInt->new("18446744073709551615");
+ return '' if $n > $_smallval;
return '' if is_prob_prime($n) != 2;
return "Type Small\nN $n";
}
if ($method eq 'Pratt' || $method eq 'Lucas') {
- return '' if scalar @$pdata != 2 ||
- ref($$pdata[0]) ne 'ARRAY' ||
- ref($$pdata[1]) eq 'ARRAY';
+ 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 '';
+ }
my @factors = @{shift @$pdata};
my $a = shift @$pdata;
my $cert = "Type Lucas\nN $n\n";
@@ -266,12 +285,16 @@ sub _convert_cert {
if (scalar @$pdata == 3 && ref($$pdata[0]) eq 'ARRAY' && $$pdata[0]->[0] =~ /^(B|T7|Theorem\s*7)$/i) {
croak "Unsupported BLS7 proof in conversion";
}
- return '' if scalar @$pdata != 2 ||
- ref($$pdata[0]) ne 'ARRAY' ||
- ref($$pdata[1]) ne 'ARRAY';
+ 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 '';
+ }
my @factors = @{shift @$pdata};
my @as = @{shift @$pdata};
- return '' unless scalar @factors == scalar @as;
+ if (scalar @factors != scalar @as) {
+ carp "verify_prime: incorrect n-1 format, must have a value for each factor\n";
+ return '';
+ }
# Make sure 2 is at the top
foreach my $i (1 .. $#factors) {
my $f = (ref($factors[$i]) eq 'ARRAY') ? $factors[$i]->[0] : $factors[$i];
@@ -298,13 +321,27 @@ sub _convert_cert {
return $cert;
}
if ($method eq 'ECPP' || $method eq 'AGKM') {
- return '' if scalar @$pdata < 1;
+ if (scalar @$pdata < 1) {
+ carp "verify_prime: incorrect AGKM format\n";
+ return '';
+ }
my $cert = '';
+ my $q = $n;
foreach my $block (@$pdata) {
- return '' if ref($block) ne 'ARRAY' || scalar @$block != 6;
+ if (ref($block) ne 'ARRAY' || scalar @$block != 6) {
+ carp "verify_prime: incorrect AGKM block format\n";
+ return '';
+ }
my($ni, $a, $b, $m, $qval, $P) = @$block;
- my $q = ref($qval) eq 'ARRAY' ? $qval->[0] : $qval;
- return '' if ref($P) ne 'ARRAY' || scalar @$P != 2;
+ if (Math::BigInt->new("$ni") != Math::BigInt->new("$q")) {
+ carp "verify_prime: incorrect AGKM block format: block n != q\n";
+ return '';
+ }
+ $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 '';
+ }
my ($x, $y) = @{$P};
$cert .= "Type ECPP\nN $ni\nA $a\nB $b\nM $m\nQ $q\nX $x\nY $y\n\n";
if (ref($qval) eq 'ARRAY') {
@@ -313,6 +350,7 @@ sub _convert_cert {
}
return $cert;
}
+ carp "verify_prime: Unknown method: '$method'.\n";
return '';
}
@@ -338,12 +376,12 @@ sub convert_array_cert_to_string {
# Verify certificate
###############################################################################
-sub _primality_error ($) {
+sub _primality_error ($) { ## no critic qw(ProhibitSubroutinePrototypes)
print "primality fail: $_[0]" if prime_get_config->{'verbose'};
return; # error in certificate
}
-sub _pfail ($) {
+sub _pfail ($) { ## no critic qw(ProhibitSubroutinePrototypes)
print "primality fail: $_[0]" if prime_get_config->{'verbose'};
return; # Failed a condition
}
@@ -638,7 +676,7 @@ sub _verify_pock {
sub _verify_small {
my ($n) = @_;
return unless defined $n;
- return _pfail "Small n $n is > 2^64\n" if $n > Math::BigInt->new("18446744073709551615");
+ return _pfail "Small n $n is > 2^64\n" if $n > $_smallval;
return _pfail "Small n $n does not pass BPSW" unless is_prob_prime($n);
($n);
}
@@ -729,7 +767,6 @@ sub verify_cert {
POCKLINGTON => \&_prove_pock, # simple n-1, Primo type 1
LUCAS => \&_prove_lucas, # n-1 completely factored
);
- my $smallval = Math::BigInt->new(2)->bpow(64);
my $base = 10;
my $cert_type = 'Unknown';
my $N;
@@ -775,7 +812,7 @@ sub verify_cert {
my $q = shift @qs;
# Check that this q has a chain
if (!defined $parts{$q}) {
- if ($q > $smallval) {
+ if ($q > $_smallval) {
_primality_error "q value $q has no proof\n";
return 0;
}
@@ -807,7 +844,7 @@ __END__
=head1 NAME
-Math::Prime::Util::PrimalityProving - Support for primality proofs and certs
+Math::Prime::Util::PrimalityProving - Primality proofs and certificates
=head1 VERSION
@@ -838,6 +875,35 @@ 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 convert_array_cert_to_string
+
+Takes as input a Perl structure certificate, used by Math::Prime::Util
+from version 0.26 through 0.29, and converts it to a multi-line text
+certificate starting with "[MPU - Primality Certificate]". This is the
+new format produced and processed by Math::Prime::Util, Math::Prime::Util::GMP,
+and associated tools.
+
+=head2 verify_cert
+
+Takes a MPU primality certificate and verifies that it does prove the
+primality of the number it represents (the N after the "Proof for:" line).
+For backwards compatibility, if given an old-style Perl structure, it will
+be converted then verified.
+
+The return value will be C<0> (failed to verify) or C<1> (verified).
+A result of C<0> does I<not> indicate the number is composite; it only
+indicates the proof given is not sufficient.
+
+If the certificate is malformed, the routine will carp a warning in addition
+to returning 0. If the C<verbose> option is set (see L</prime_set_config>)
+then if the validation fails, the reason for the failure is printed in
+addition to returning 0. If the C<verbose> option is set to 2 or higher, then
+a message indicating success and the certificate type is also printed.
+
+A later release may add support for
+L<Primo|http://www.ellipsa.eu/public/primo/primo.html>
+certificates, as all the method verifications are coded.
+
=head1 SEE ALSO
diff --git a/t/23-primality-proofs.t b/t/23-primality-proofs.t
index 2195fea..6a0fbc8 100644
--- a/t/23-primality-proofs.t
+++ b/t/23-primality-proofs.t
@@ -61,12 +61,20 @@ foreach my $p (@plist) {
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 $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($cert2), " verification of prime_certificate done" );
+ ok( defined($cert2) && $cert2 =~ /^Type/m,
+ " prime_certificate is also non-null" );
+ if ($cert2 eq $cert) {
+ ok(1, " certificate is identical to first");
+ } else {
+ ok( verify_prime($cert2), " different cert, verified" );
+ }
}
}
+# TODO: All these proofs are using the old format. That's ok for now,
+# as verify_prime will convert them for us. But we really should do
+# testing with the new format, including possible errors it could have.
+
# Some hand-done proofs
SKIP: {
skip "Skipping 2**521-1 verification without Math::BigInt::GMP", 1
diff --git a/t/80-pp.t b/t/80-pp.t
index 6c571e3..e4b2c8a 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -251,7 +251,7 @@ plan tests => 2 +
1 + 1 + # factor
10 + 8*3 + # factoring subs
10 + # AKS
- 3 + # Lucas and BLS75 primality proofs
+ 2 + # Lucas and BLS75 primality proofs
3 + # M-R and Lucas on bigint
13 + # Misc util.pm functions
scalar(keys %ipp) +
@@ -566,23 +566,15 @@ SKIP: {
}
is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_lucas(100003)],
- [2, [100003, "Pratt", [2, 3, 7, 2381], 2]],
+ [2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN 100003\n\nType Lucas\nN 100003\nQ[1] 2\nQ[2] 3\nQ[3] 7\nQ[4] 2381\nA 2\n"],
"primality_proof_lucas(100003)" );
# Had to reduce these to make borked up Perl 5.6.2 work.
#is_deeply( [Math::Prime::Util::PP::primality_proof_bls75("210596120454733723")],
# [2, ["210596120454733723", "n-1", [2, 3, 82651, "47185492693"], [2, 2, 2, 2]]],
# "primality_proof_bls75(210596120454733723)" );
-#is_deeply( [Math::Prime::Util::PP::primality_proof_bls75("809120722675364249")],
-# [2, ["809120722675364249", "n-1",
-# ["B", 20000, "22233477760919", 2], [2, 4549], [3, 2]]],
-# "primality_proof_bls75(809120722675364249)" );
is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_bls75(1490266103)],
- [2, [1490266103, 'n-1', [2, 13, 19, 1597, 1889], [5, 2, 2, 2, 2]]],
+ [2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN 1490266103\n\nType BLS5\nN 1490266103\nQ[1] 13\nQ[2] 19\nQ[3] 1597\nQ[4] 1889\nA[0] 5\n----\n"],
"primality_proof_bls75(1490266103)" );
-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)" );
{
my $n = Math::BigInt->new("168790877523676911809192454171451");
--
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