[libmath-prime-util-perl] 12/18: Primality proof updates
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:48:03 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.26
in repository libmath-prime-util-perl.
commit aa43b9b989ab0d56daad9eeb0ba4e84d768f3861
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Apr 17 19:30:40 2013 -0700
Primality proof updates
---
lib/Math/Prime/Util.pm | 177 ++++++++++++++++++++-----------
lib/Math/Prime/Util/ECProjectivePoint.pm | 25 +++--
lib/Math/Prime/Util/PP.pm | 58 +++++++---
3 files changed, 179 insertions(+), 81 deletions(-)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index e023cb2..e7e4d47 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -15,7 +15,7 @@ use base qw( Exporter );
our @EXPORT_OK =
qw( prime_get_config prime_set_config
prime_precalc prime_memfree
- is_prime is_prob_prime is_provable_prime
+ is_prime is_prob_prime is_provable_prime is_provable_prime_with_cert
prime_certificate verify_prime
is_strong_pseudoprime is_strong_lucas_pseudoprime
is_aks_prime
@@ -26,7 +26,7 @@ our @EXPORT_OK =
prime_count_lower prime_count_upper prime_count_approx
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime random_nbit_prime
- random_strong_prime random_maurer_prime
+ random_strong_prime random_maurer_prime random_maurer_prime_with_cert
primorial pn_primorial consecutive_integer_lcm
factor all_factors
moebius mertens euler_phi jordan_totient exp_mangoldt
@@ -834,10 +834,21 @@ 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);
+ return $n;
+ }
+
+ sub random_maurer_prime_with_cert {
my($k) = @_;
_validate_positive_integer($k, 2);
+ my @cert;
if ($] < 5.008 && $_Config{'maxbits'} > 32) {
- return random_nbit_prime($k) if $k <= 49;
+ if ($k <= 49) {
+ my $n = random_nbit_prime($k);
+ return ($n, [$n]);
+ }
croak "Random Maurer not supported on old Perl";
}
@@ -846,7 +857,12 @@ sub primes {
# returning 2. This should be reasonably fast to ~128 bits with MPU::GMP.
my $p0 = $_Config{'maxbits'};
- return random_nbit_prime($k) if $k <= $p0;
+ if ($k <= $p0) {
+ my $n = random_nbit_prime($k);
+ my ($isp, $cert) = is_provable_prime_with_cert($n);
+ croak "small nbit prime could not be proven" if $isp != 2;
+ return ($n, $cert);
+ }
if (!defined $Math::BigInt::VERSION) {
eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
@@ -877,7 +893,7 @@ sub primes {
}
# I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1.
- my $q = random_maurer_prime( ($r * $k)->bfloor + 1 );
+ my ($q, $certref) = 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;
@@ -955,7 +971,15 @@ sub primes {
# 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);
- return $n;
+ # 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);
+ }
}
# Didn't pass the selected a values. Try another R.
}
@@ -1600,46 +1624,51 @@ sub is_prob_prime {
return ($n <= 18446744073709551615) ? 2 : 1;
}
-
+# Return just the non-cert portion.
sub is_provable_prime {
- my($n, $ref_proof) = @_;
+ my($n) = @_;
return 0 if defined $n && $n < 2;
_validate_positive_integer($n);
- if (defined $ref_proof) {
- croak "Second argument must be an array ref" if ref($ref_proof) ne 'ARRAY';
- @$ref_proof = ();
- }
+ return _XS_is_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
+ return Math::Prime::Util::GMP::is_provable_prime($n)
+ if $_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime;
+
+ my ($is_prime, $cert) = is_provable_prime_with_cert($n);
+ return $is_prime;
+}
+
+# Return just the cert portion.
+sub prime_certificate {
+ my($n) = @_;
+ my ($is_prime, $cert) = is_provable_prime_with_cert($n);
+ return @$cert;
+}
+
+
+sub is_provable_prime_with_cert {
+ my($n) = @_;
+ return 0 if defined $n && $n < 2;
+ _validate_positive_integer($n);
# Set to 0 if you want the proof to go down to 11.
if (1) {
- if (defined $ref_proof) {
- if (ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL) {
- my $isp = _XS_is_prime($n);
- @$ref_proof = ($n) if $isp == 2;
- return $isp;
- }
- if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime_with_cert) {
- my($isp, $pref) = Math::Prime::Util::GMP::is_provable_prime_with_cert($n);
- @$ref_proof = @$pref;
- return $isp;
- }
- } else {
- return _XS_is_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
- return Math::Prime::Util::GMP::is_provable_prime($n)
- if $_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime;
+ 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);
}
- my $is_prob_prime = is_prob_prime($n);
- if ($is_prob_prime != 1) {
- @$ref_proof = ($n) if defined $ref_proof && $is_prob_prime == 2;
- return $is_prob_prime;
+ 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) {
- @$ref_proof = ($n) if defined $ref_proof;
- return 2;
+ return (2, [$n]);
}
return 0;
}
@@ -1648,32 +1677,23 @@ sub is_provable_prime {
# Choice of methods for proof:
# ECPP needs a fair bit of programming work
# APRCL needs a lot of programming work
- # BLS75 Requires factoring n-1 to (n/2)^1/3
- # Pocklington Requires factoring n-1 to n^1/2
- # Lucas Easy, required complete factoring of n-1
+ # BLS75 combo Corollary 11 of BLS75. Trial factor n-1 and n+1 to B, find
+ # factors F1 of n-1 and F2 of n+1. Quit when:
+ # B > (N/(F1*F1*(F2/2)))^1/3 or B > (N/((F1/2)*F2*F2))^1/3
+ # BLS75 n+1 Requires factoring n+1 to (n/2)^1/3 (theorem 19)
+ # BLS75 n-1 Requires factoring n-1 to (n/2)^1/3 (theorem 5 or 7)
+ # Pocklington Requires factoring n-1 to n^1/2 (BLS75 theorem 4)
+ # Lucas Easy, requires factoring of n-1 (BLS75 theorem 1)
# 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);
-
- @$ref_proof = @$pref if defined $ref_proof;
carp "proved $n is not prime\n" if !$isp;
- return $isp;
+ return ($isp, $pref);
}
-sub prime_certificate {
- my($n) = @_;
- return () if defined $n && $n < 2;
- _validate_positive_integer($n);
-
- my @cert;
- my $is_prime = is_provable_prime($n, \@cert);
- return () unless $is_prime == 2;
- return @cert;
-}
-
sub verify_prime {
my @pdata = @_;
my $verbose = $_Config{'verbose'};
@@ -1681,6 +1701,9 @@ sub verify_prime {
# Empty input = no certificate = not prime
return 0 if scalar @pdata == 0;
+ # Handle case of being handed a reference to the certificate.
+ @pdata = @{$pdata[0]} if scalar @pdata == 1 && ref($pdata[0]) eq 'ARRAY';
+
my $n = shift @pdata;
if (length($n) == 1) {
return 1 if $n =~ /^[2357]$/;
@@ -1804,7 +1827,7 @@ sub verify_prime {
$f = Math::BigInt->new("$farray->[0]");
return 0 unless verify_prime(@$farray);
} else {
- $f = $farray;
+ $f = Math::BigInt->new("$farray");
return 0 unless verify_prime($f);
}
next if defined $factors_seen{"$f"}; # repeated factors
@@ -1821,7 +1844,13 @@ sub verify_prime {
$factors_seen{"$f"} = 1;
}
croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
- croak "BLS75 error: $A not even" unless $A->is_even();
+
+ # 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;
+ }
# TODO: consider: if B=1 and a single a is given, then Lucas test.
@@ -2777,18 +2806,16 @@ exist, there is a weak conjecture (Martin) that none exist under 10000 digits.
Takes a positive number as input and returns back either 0 (composite),
2 (definitely prime), or 1 (probably prime). This gives it the same return
-values as C<is_prime> and C<is_prob_prime>.
+values as L</is_prime> and L</is_prob_prime>.
-The current implementation uses a Lucas test requiring a complete factorization
-of C<n-1>, which may not be possible in a reasonable amount of time. The GMP
-version uses the BLS (Brillhart-Lehmer-Selfridge) method, requiring C<n-1> to
-be factored to the cube root of C<n>, which is more likely to succeed and will
-usually take less time, but can still fail. Hence you should always test that
-the result is C<2> to ensure the prime is proven.
+The current implementation of both the Perl and GMP proofs is using theorem 5
+of BLS75 (Brillhart-Lehmer-Selfridge), requiring C<n-1> to be factored to
+C<(n/2)^(1/3))>. This takes less time than factoring to C<n^0.5> as required
+by the generalized Pocklington test or C<n-1> for the Lucas test. However it
+is possible a factor cannot be found in a reasonable amount of time, so you
+should always test that the result in C<2> to ensure it was proven.
-An optional second argument may be given, which must be an array reference,
-which will be filled in with a primality certificate. Normally this is used
-via the function L</prime_certificate> below.
+A later implementation will use an ECPP test for larger inputs.
=head2 prime_certificate
@@ -2802,6 +2829,14 @@ This may be examined or given to L</verify_prime> for verification. The latter
function contains the description of the format.
+=head2 is_provable_prime_with_cert
+
+Given a positive integer as input, returns a two element array containing the
+result of L</is_provable_prime> and an array reference containing the primality
+certificate like L</prime_certificate>. The certificate will be an empty
+array reference if the result is not 2 (definitely prime).
+
+
=head2 verify_prime
my @cert = prime_certificate($n);
@@ -3234,6 +3269,7 @@ Similar to L</random_nbit_prime>, the result will be a BigInt if the
number of bits is greater than the native bit size. For better performance
with large bit sizes, install L<Math::Prime::Util::GMP>.
+
=head2 random_maurer_prime
my $bigprime = random_maurer_prime(512);
@@ -3247,6 +3283,23 @@ with large bit sizes, install L<Math::Prime::Util::GMP>.
The differences between this function and that in L<Crypt::Primes> are
described in the L</"SEE ALSO"> section.
+Internally this additionally runs the BPSW probable prime test on every
+partial result, and constructs a primality certificate for the final
+result, which is verified. These add additional checks that the resulting
+value has been properly constructed.
+
+
+=head2 random_maurer_prime_with_cert
+
+ my($n, $cert_ref) = random_maurer_prime_with_cert(512)
+
+As with L</random_maurer_prime>, but returns a two-element array containing
+the n-bit provable prime along with a primality certificate. The certificate
+is the same as produced by L</prime_certificate> or
+L</is_provable_prime_with_cert>, and can be parsed by L</verify_prime> or
+any other software that can parse the certificate (the "n-1" form is described
+in detail in L</verify_prime>).
+
=head1 UTILITY FUNCTIONS
diff --git a/lib/Math/Prime/Util/ECProjectivePoint.pm b/lib/Math/Prime/Util/ECProjectivePoint.pm
index 5af473f..1d062d8 100644
--- a/lib/Math/Prime/Util/ECProjectivePoint.pm
+++ b/lib/Math/Prime/Util/ECProjectivePoint.pm
@@ -187,7 +187,7 @@ sub copy {
__END__
-# ABSTRACT: Elliptic curve operations for affine points
+# ABSTRACT: Elliptic curve operations for projective points
=pod
@@ -196,7 +196,7 @@ __END__
=head1 NAME
-Math::Prime::Util::ECAffinePoint - Elliptic curve operations for affine points
+Math::Prime::Util::ECProjectivePoint - Elliptic curve operations for projective points
=head1 VERSION
@@ -207,7 +207,7 @@ Version 0.26
=head1 SYNOPSIS
# Create a point on a curve (a,b,n) with coordinates 0,1
- my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $n, 0, 1);
+ my $ECP = Math::Prime::Util::ECProjectivePoint->new($a, $b, $n, 0, 1);
# scalar multiplication by k.
$ECP->mul($k)
@@ -228,7 +228,7 @@ To write.
=head2 new
- $point = Math::Prime::Util::EllipticCurve->new(a, b);
+ $point = Math::Prime::Util::ECProjectivePoint->new(a, b);
Returns a new curve defined by a and b.
@@ -242,9 +242,9 @@ Returns the C<a>, C<b>, or C<n> values that describe the curve.
=head2 x
-=head2 y
+=head2 z
-Returns the C<x> or C<y> values that define the point on the curve.
+Returns the C<x> or C<z> values that define the point on the curve.
=head2 f
@@ -254,6 +254,10 @@ Returns a possible factor found during EC multiplication.
Takes another point on the same curve as an argument and adds it this point.
+=head2 double
+
+Double the current point on the curve.
+
=head2 mul
Takes an integer and performs scalar multiplication of the point.
@@ -263,6 +267,15 @@ Takes an integer and performs scalar multiplication of the point.
Returns true if the point is (0,1), which is the point at infinity for
the affine coordinates.
+=head2 copy
+
+Returns a copy of the point.
+
+=head2 normalize
+
+Performs an extended gcd operation to make C<z=1>. If a factor of C<n> is
+found it is put in C<f>.
+
=head1 SEE ALSO
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index b94e684..d8a7bcd 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1815,7 +1815,7 @@ sub primality_proof_lucas {
# Since this can take a very long time with a composite, try some easy cuts
return @composite if !defined $n || $n < 2;
- return ($n,[$n]) if $n < 4;
+ return (2, [$n]) if $n < 4;
return @composite if miller_rabin($n,2,15,325) == 0;
if (!defined $Math::BigInt::VERSION) {
@@ -1843,12 +1843,12 @@ sub primality_proof_lucas {
# Verify each factor and add to proof
my @fac_proofs;
foreach my $f (@factors) {
- my @fproof;
- if (Math::Prime::Util::is_provable_prime($f, \@fproof) != 2) {
+ 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 : [@fproof];
+ push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
}
my @proof = ("$n", "Pratt", [@fac_proofs], $a);
return (2, [@proof]);
@@ -1856,13 +1856,15 @@ sub primality_proof_lucas {
return @composite;
}
+use Data::Dump qw/dump/;
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 ($n,[$n]) if $n < 4;
+ 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) {
@@ -1874,9 +1876,11 @@ sub primality_proof_bls75 {
my $nm1 = $n-1;
my $A = $nm1->copy->bone; # factored part
my $B = $nm1->copy; # unfactored part
- my @factors;
+ 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);
@@ -1884,7 +1888,9 @@ sub primality_proof_bls75 {
}
my @nstack;
# nstack should only hold composites
- if (Math::Prime::Util::is_prob_prime($B)) {
+ 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 {
@@ -1925,8 +1931,7 @@ sub primality_proof_bls75 {
}
}
{ # remove duplicate factors and make a sorted array of bigints
- my %uf;
- undef @uf{@factors};
+ my %uf = map { $_ => 1 } @factors;
@factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
}
# Did we factor enough?
@@ -1955,14 +1960,14 @@ sub primality_proof_bls75 {
last;
}
return @composite unless $success;
- my @fproof;
- if (Math::Prime::Util::is_provable_prime($f, \@fproof) != 2) {
+ 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 : [@fproof];
+ push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
}
- my @proof = ("$n", "n-1", [@fac_proofs], [@as]);
+ my @proof = ($n, "n-1", [@fac_proofs], [@as]);
return (2, [@proof]);
}
@@ -2625,6 +2630,20 @@ 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.
+
=head1 UTILITY FUNCTIONS
@@ -2730,6 +2749,19 @@ given is to run with increasing B1 factors.
This method can rapidly find a factor C<p> of C<n> where C<p-1> is smooth
(i.e. C<p-1> has no large factors).
+=head2 ecm_factor
+
+ my @factors = ecm_factor($n);
+ my @factors = ecm_factor($n, 1000); # B1 = B2 = 1000, curves = 10
+ my @factors = ecm_factor($n, 1000, 5000, 10); # Set B1, B2, and ncurves
+
+Given a positive number input, tries to discover a factor using ECM. The
+resulting array will contain either two factors (it succeeded) or the original
+number (no factor was found). In either case, multiplying @factors yields the
+original input. The B1 and B2 smoothness factors for stage 1 and stage 2
+respectively may be supplied, as can be number of curves to try.
+
+
=head1 MATHEMATICAL FUNCTIONS
--
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