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