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 <[email protected]> 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 _______________________________________________ Pkg-perl-cvs-commits mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
