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 b83566207a51be1f6a34f575c92ce179e6f80279 Author: Dana Jacobsen <d...@acm.org> Date: Mon Apr 15 00:15:08 2013 -0700 PP: add simple ECM factoring and BLS75 primality proof --- Changes | 10 +- MANIFEST | 1 + lib/Math/Prime/Util.pm | 114 +++++-------- lib/Math/Prime/Util/ECAffinePoint.pm | 236 ++++++++++++++++++++++++++ lib/Math/Prime/Util/EllipticCurve.pm | 7 +- lib/Math/Prime/Util/PP.pm | 317 +++++++++++++++++++++++++++++------ 6 files changed, 554 insertions(+), 131 deletions(-) diff --git a/Changes b/Changes index f89a628..320f7e3 100644 --- a/Changes +++ b/Changes @@ -2,12 +2,20 @@ Revision history for Perl extension Math::Prime::Util. 0.26 xx April 2013 - - Pure Perl p-1 factoring, Fermat factoring, and speedup for pbrent. + - Pure Perl factoring: + - real p-1 -- much faster and more effective + - Fermat (no better than HOLF) + - speedup for pbrent + - simple affine single stage ECM + - redo factoring mix - New functions: prime_certificate produces a certificate of primality. verify_prime checks a primality certificate. + - Pure perl primality proof now uses BLS75 instead of Lucas, so some + numbers will be much faster. n-1 only needs factoring to (n/2)^1/3. + - Math::Prime::Util::EllipticCurve module. 0.25 19 March 2013 diff --git a/MANIFEST b/MANIFEST index bdd9012..416fbfa 100644 --- a/MANIFEST +++ b/MANIFEST @@ -5,6 +5,7 @@ lib/Math/Prime/Util/PrimeArray.pm lib/Math/Prime/Util/PP.pm lib/Math/Prime/Util/ZetaBigFloat.pm lib/Math/Prime/Util/EllipticCurve.pm +lib/Math/Prime/Util/ECAffinePoint.pm LICENSE Makefile.PL MANIFEST diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index 0243ade..ce49022 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -1600,6 +1600,7 @@ sub is_prob_prime { return ($n <= 18446744073709551615) ? 2 : 1; } + sub is_provable_prime { my($n, $ref_proof) = @_; return 0 if defined $n && $n < 2; @@ -1612,17 +1613,21 @@ sub is_provable_prime { # 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); - @$ref_proof = ($n) if defined $ref_proof && $isp; - return $isp; - } - if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime) { - return Math::Prime::Util::GMP::is_provable_prime($n) if !defined $ref_proof; - if (defined $ref_proof && $Math::Prime::Util::GMP::VERSION > 0.08) { - return Math::Prime::Util::GMP::is_provable_prime($n, $ref_proof); + 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; } - # proof needed but MPU::GMP too old to give it. + 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; } my $is_prob_prime = is_prob_prime($n); @@ -1640,63 +1645,24 @@ sub is_provable_prime { } } - # At this point we know it is almost certainly a prime, but we need to - # prove it. We should do ECPP or APR-CL now, or failing that, do the - # Brillhart-Lehmer-Selfridge test, or Pocklington-Lehmer. Until those - # are written here, we'll do a Lucas test, which is super simple but may - # be very slow. We have AKS code, but it's insanely slow. + # 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 + # AKS horribly slow # See http://primes.utm.edu/prove/merged.html or other sources. - # It shouldn't be possible to get here without BigInt already loaded. - 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); - # If not doing a proof, check all factors here. - if (!defined $ref_proof) { - if ( (scalar grep { is_provable_prime($_) != 2 } @factors) > 0) { - carp "could not prove primality of $n.\n"; - return 1; - } - } - { # remove duplicate factors - my %uf; - undef @uf{@factors}; - @factors = sort {$a <=> $b} keys %uf; - } + #my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_lucas($n); + my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_bls75($n); - 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; - # If doing a proof, we verify each factor here and add to proof. - if (defined $ref_proof) { - @$ref_proof = (); - my @fac_proofs; - foreach my $f (@factors) { - my @fproof; - if (is_provable_prime($f, \@fproof) != 2) { - carp "could not prove primality of $n.\n"; - return 1; - } - push @fac_proofs, (scalar @fproof == 1) ? @fproof : [@fproof]; - } - @$ref_proof = ("$n", "Pratt", [@fac_proofs], $a); - } - return 2; - } - carp "proved $n is not prime\n"; - return 0; + @$ref_proof = @$pref if defined $ref_proof; + carp "proved $n is not prime\n" if !$isp; + return $isp; } + sub prime_certificate { my($n) = @_; return () if defined $n && $n < 2; @@ -1937,21 +1903,23 @@ sub verify_prime { warn "verify_prime: AGKM block q is too small\n"; return 0; } - if (!defined $Math::Prime::Util::EllipticCurve::VERSION) { - eval { require Math::Prime::Util::EllipticCurve; 1; } - or do { croak "Cannot load Math::Prime::Util::EllipticCurve"; }; + # Final check, check that we've got a bound above and below (Hasse) + if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) { + eval { require Math::Prime::Util::ECAffinePoint; 1; } + or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; }; } - my $EC = Math::Prime::Util::EllipticCurve->new($a, $b, $ni); $m = Math::BigInt->new("$m") unless ref($m) eq 'Math::BigInt'; $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt'; - - # Assume P0 and P1 are affine. - my $Px = Math::BigInt->new($P->[0]); - my $Py = Math::BigInt->new($P->[1]); + my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $ni, $P->[0], $P->[1]); # Compute U = (m/q)P, check U != point at infinity - my($Ux,$Uy) = $EC->mul_a( int($m/$q), $Px, $Py ); # U = (m/q)P - my($Vx,$Vy) = $EC->mul_a( $q, $Ux, $Uy ); # V = qU - if ( (($Ux == 0) && ($Uy == 1)) || (($Vx != 0) || ($Vy != 1)) ) { + $ECP->mul( int($m/$q) ); + if ($ECP->is_infinity) { + warn "verify_prime: AGKM point does not multiply correctly.\n"; + return 0; + } + # Compute V = qU, check V = point at infinity + $ECP->mul( $q ); + if (! $ECP->is_infinity) { warn "verify_prime: AGKM point does not multiply correctly.\n"; return 0; } diff --git a/lib/Math/Prime/Util/ECAffinePoint.pm b/lib/Math/Prime/Util/ECAffinePoint.pm new file mode 100644 index 0000000..81c8d06 --- /dev/null +++ b/lib/Math/Prime/Util/ECAffinePoint.pm @@ -0,0 +1,236 @@ +package Math::Prime::Util::ECAffinePoint; +use strict; +use warnings; +use Carp qw/carp croak confess/; + +if (!defined $Math::BigInt::VERSION) { + eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; } + or do { croak "Cannot load Math::BigInt"; }; +} + +BEGIN { + $Math::Prime::Util::EllipticCurve::AUTHORITY = 'cpan:DANAJ'; + $Math::Prime::Util::EllipticCurve::VERSION = '0.26'; +} + +# Pure perl (with Math::BigInt) manipulation of Elliptic Curves +# in affine coordinates. Should be split into a point class. + +sub new { + my ($class, $a, $b, $n, $x, $y) = @_; + $a = Math::BigInt->new("$a") unless ref($a) eq 'Math::BigInt'; + $b = Math::BigInt->new("$b") unless ref($b) eq 'Math::BigInt'; + $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt'; + $x = Math::BigInt->new("$x") unless ref($x) eq 'Math::BigInt'; + $y = Math::BigInt->new("$y") unless ref($y) eq 'Math::BigInt'; + + croak "n must be >= 2" unless $n >= 2; + $a->bmod($n); + $b->bmod($n); + + my $self = { + a => $a, + b => $b, + n => $n, + x => $x, + y => $y, + f => $n->copy->bone, + }; + + bless $self, $class; + return $self; +} + +sub _add { + my ($self, $P1x, $P1y, $P2x, $P2y) = @_; + my $n = $self->{'n'}; + + if ($P1x == $P2x) { + my $t = ($P1y + $P2y) % $n; + return (Math::BigInt->bzero,Math::BigInt->bone) if $t == 0; + } + my $deltax = ($P2x - $P1x) % $n; + $deltax->bmodinv($n); + return (Math::BigInt->bzero,Math::BigInt->bone) if $deltax eq "NaN"; + + my $deltay = ($P2y - $P1y) % $n; + my $m = ($deltay * $deltax) % $n; # m = deltay / deltax + + my $x = ($m*$m - $P1x - $P2x) % $n; + my $y = ($m*($P1x - $x) - $P1y) % $n; + return ($x,$y); +} + +sub _double { + my ($self, $P1x, $P1y) = @_; + my $n = $self->{'n'}; + + my $m = 2*$P1y; + $m->bmodinv($n); + return (Math::BigInt->bzero,Math::BigInt->bone) if $m eq "NaN"; + + $m = ((3*$P1x*$P1x + $self->{'a'}) * $m) % $n; + + my $x = ($m*$m - 2*$P1x) % $n; + my $y = ($m*($P1x - $x) - $P1y) % $n; + return ($x,$y); +} + +sub mul { + my ($self, $k) = @_; + my $x = $self->{'x'}; + my $y = $self->{'y'}; + my $n = $self->{'n'}; + my $f = $self->{'f'}; + + my $Bx = $n->copy->bzero; + my $By = $n->copy->bone; + my $v = 1; + + while ($k > 0) { + if ( ($k % 2) != 0) { + $k--; + $f->bmul($Bx-$x)->bmod($n); + if ($x == 0 && $y == 1) { } + elsif ($Bx == 0 && $By == 1) { ($Bx,$By) = ($x,$y); } + else { ($Bx,$By) = $self->_add($x,$y,$Bx,$By); } + } else { + $k >>= 1; + $f->bmul(2*$y)->bmod($n); + ($x,$y) = $self->_double($x,$y); + } + } + $f = Math::BigInt::bgcd( $f, $n); + $self->{'x'} = $Bx; + $self->{'y'} = $By; + $self->{'f'} = $f; + return $self; +} + +sub add { + my ($self, $other) = @_; + croak "add takes a EC point" + unless ref($other) eq 'Math::Prime::Util::ECAffinePoint'; + croak "second point is not on the same curve" + unless $self->{'a'} == $other->{'a'} && + $self->{'b'} == $other->{'b'} && + $self->{'n'} == $other->{'n'}; + + ($self->{'x'}, $self->{'y'}) = $self->_add($self->{'x'}, $self->{'y'}, + $other->{'x'}, $other->{'y'}); + return $self; +} + + +sub a { return shift->{'a'}; } +sub b { return shift->{'b'}; } +sub n { return shift->{'n'}; } +sub x { return shift->{'x'}; } +sub y { return shift->{'y'}; } +sub f { return shift->{'f'}; } + +sub is_infinity { + my $self = shift; + return ($self->{'x'}->is_zero() && $self->{'y'}->is_one()); +} + +1; + +__END__ + + +# ABSTRACT: Elliptic curve operations for affine points + +=pod + +=encoding utf8 + + +=head1 NAME + +Math::Prime::Util::ECAffinePoint - Elliptic curve operations for affine points + + +=head1 VERSION + +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); + + # scalar multiplication by k. + $ECP->mul($k) + + # add two points on the same curve + $ECP->add($ECP2) + + say "P = O" if $ECP->is_infinity(); + +=head1 DESCRIPTION + +This really should just be in Math::EllipticCurve. + +To write. + + +=head1 FUNCTIONS + +=head2 new + + $point = Math::Prime::Util::EllipticCurve->new(a, b); + +Returns a new curve defined by a and b. + +=head2 a + +=head2 b + +=head2 n + +Returns the C<a>, C<b>, or C<n> values that describe the curve. + +=head2 x + +=head2 y + +Returns the C<x> or C<y> values that define the point on the curve. + +=head2 f + +Returns a possible factor found during EC multiplication. + +=head2 add + +Takes another point on the same curve as an argument and adds it this point. + +=head2 mul + +Takes an integer and performs scalar multiplication of the point. + +=head2 is_infinity + +Returns true if the point is (0,1), which is the point at infinity for +the affine coordinates. + + +=head1 SEE ALSO + +L<Math::EllipticCurve::Prime> + +This really should just be in a L<Math::EllipticCurve> module. + +=head1 AUTHORS + +Dana Jacobsen E<lt>d...@acm.orge<gt> + + +=head1 COPYRIGHT + +Copyright 2012-2013 by Dana Jacobsen E<lt>d...@acm.orge<gt> + +This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself. + +=cut diff --git a/lib/Math/Prime/Util/EllipticCurve.pm b/lib/Math/Prime/Util/EllipticCurve.pm index fff5fef..d04601b 100644 --- a/lib/Math/Prime/Util/EllipticCurve.pm +++ b/lib/Math/Prime/Util/EllipticCurve.pm @@ -136,23 +136,24 @@ sub mul_a { my $Bx = Math::BigInt->bzero; my $By = Math::BigInt->bone; my $v = 1; + my $d = 1; while ($v && $k > 0) { if ( ($k % 2) != 0) { $k--; - my $d = Math::BigInt::bgcd( ($Bx - $x) % $n, $n); + $d = Math::BigInt::bgcd( ($Bx - $x) % $n, $n); $v = ($d == 1 || $d == $n); if ($x == 0 && $y == 1) { } elsif ($Bx == 0 && $By == 1) { ($Bx,$By) = ($x,$y); } elsif ($v) { ($Bx,$By) = $self->add_a($x,$y,$Bx,$By); } } else { $k >>= 1; - my $d = Math::BigInt::bgcd( 2*$y % $n, $n); + $d = Math::BigInt::bgcd( 2*$y % $n, $n); $v = ($d == 1 || $d == $n); if ($v) { ($x,$y) = $self->double_a($x,$y); } } } - return ($Bx, $By); + return ($Bx, $By, $d); } 1; diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm index 8e86250..acf1096 100644 --- a/lib/Math/Prime/Util/PP.pm +++ b/lib/Math/Prime/Util/PP.pm @@ -1137,6 +1137,29 @@ sub trial_factor { @factors; } +my $_holf_r; +my @_fsublist = ( + sub { my $n = shift; return ($n) unless $n < $_half_word; + holf_factor ($n, 64*1024, $_holf_r); $_holf_r += 64*1024; }, + sub { prho_factor (shift, 8*1024, 3) }, + sub { pbrent_factor (shift, 32*1024, 1) }, + sub { pminus1_factor(shift, 10_000); }, + sub { pminus1_factor(shift, 600_000); }, + sub { pbrent_factor (shift, 512*1024, 7) }, + sub { ecm_factor (shift, 1_000, 1_000, 10) }, + sub { pminus1_factor(shift, 4_000_000); }, + sub { pbrent_factor (shift, 512*1024, 11) }, + sub { ecm_factor (shift, 10_000, 10_000, 10) }, + sub { holf_factor (shift, 256*1024, $_holf_r); $_holf_r += 256*1024; }, + sub { pminus1_factor(shift,20_000_000); }, + sub { ecm_factor (shift, 100_000, 100_000, 10) }, + sub { holf_factor (shift, 512*1024, $_holf_r); $_holf_r += 512*1024; }, + sub { pbrent_factor (shift, 2048*1024, 13) }, + sub { holf_factor (shift, 2048*1024, $_holf_r); $_holf_r += 2048*1024; }, + sub { ecm_factor (shift, 1_000_000, 1_000_000, 10) }, + sub { pminus1_factor(shift, 100_000_000, 500_000_000); }, +); + sub factor { my($n) = @_; _validate_positive_integer($n); @@ -1168,54 +1191,10 @@ sub factor { #print "Looking at $n with stack ", join(",",@nstack), "\n"; while ( ($n >= (31*31)) && !_is_prime7($n) ) { my @ftry; - my $holf_rounds = 0; - if ($n < $_half_word) { - $holf_rounds = 64*1024; - #warn "trying holf 64k on $n\n"; - @ftry = holf_factor($n, $holf_rounds); - } - if (scalar @ftry < 2) { - #warn "trying prho 8k {3} on $n\n"; - @ftry = prho_factor($n, 8*1024, 3); - } - if (scalar @ftry < 2) { - #warn "trying pbrent 32k {1} on $n\n"; - @ftry = pbrent_factor($n, 32*1024, 1); - } - if (scalar @ftry < 2) { - #warn "trying p-1 10k on $n\n"; - @ftry = pminus1_factor($n, 10_000); - } - if (scalar @ftry < 2) { - #warn "trying p-1 1M on $n\n"; - @ftry = pminus1_factor($n, 1_000_000); - } - if (scalar @ftry < 2) { - #warn "trying pbrent 512k {7} on $n\n"; - @ftry = pbrent_factor($n, 512*1024, 7); - } - if (scalar @ftry < 2) { - #warn "trying holf 128k on $n\n"; - @ftry = holf_factor($n, 128*1024, $holf_rounds); - $holf_rounds += 128*1024; - } - if (scalar @ftry < 2) { - #warn "trying pbrent 2M {13} on $n\n"; - @ftry = pbrent_factor($n, 2048*1024, 13); - } - if (scalar @ftry < 2) { - #warn "trying holf 256k on $n\n"; - @ftry = holf_factor($n, 256*1024, $holf_rounds); - $holf_rounds += 256*1024; - } - if (scalar @ftry < 2) { - #warn "trying p-1 100M on $n\n"; - @ftry = pminus1_factor($n, 100_000_000, 500_000_000); - } - if (scalar @ftry < 2) { - #warn "trying holf 512k on $n\n"; - @ftry = holf_factor($n, 512*1024, $holf_rounds); - $holf_rounds += 512*1024; + $_holf_r = 1; + foreach my $sub (@_fsublist) { + last if scalar @ftry >= 2; + @ftry = $sub->($n); } if (scalar @ftry > 1) { #print " split into ", join(",",@ftry), "\n"; @@ -1466,7 +1445,7 @@ sub pminus1_factor { push @factors, $n; return @factors; } - $B2 = 1*$B1 unless defined $B2; + $B2 = 1*$B1 unless defined $B2; my $one = $n->copy->bone; my ($j, $q, $saveq) = (1, 2, 2); @@ -1574,6 +1553,7 @@ sub holf_factor { _validate_positive_integer($n); $rounds = 64*1024*1024 unless defined $rounds; $startrounds = 1 unless defined $startrounds; + $startrounds = 1 if $startrounds < 1; my @factors = _basic_factor($n); return @factors if $n < 4; @@ -1582,8 +1562,17 @@ sub holf_factor { for my $i ($startrounds .. $rounds) { my $ni = $n->copy->bmul($i); my $s = $ni->copy->bsqrt->bfloor->as_int; - $s->binc if ($s * $s) != $ni; - my $m = $s->copy->bmul($s)->bmod($n); + if ($s * $s == $ni) { + # s^2 = n*i, so m = s^2 mod n = 0. Hence f = GCD(n, s) = GCD(n, n*i) + my $f = Math::BigInt::bgcd($ni, $n); + last if $f == 1 || $f == $n; # Should never happen + push @factors, $f; + push @factors, int($n/$f); + croak "internal error in HOLF" unless ($f * int($n/$f)) == $n; + return @factors; + } + $s->binc; + my $m = ($s * $s) - $ni; # Check for perfect square my $mc = int(($m & 31)->bstr); next unless $mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25; @@ -1591,10 +1580,9 @@ sub holf_factor { next unless ($f*$f) == $m; $f = Math::BigInt::bgcd( ($s > $f) ? $s-$f : $f-$s, $n); last if $f == 1 || $f == $n; # Should never happen - my $f2 = $n->copy->bdiv($f)->as_int; push @factors, $f; - push @factors, $f2; - croak "internal error in HOLF" unless ($f * $f2) == $n; + push @factors, int($n/$f); + croak "internal error in HOLF" unless ($f * int($n/$f)) == $n; # print "HOLF found factors in $i rounds\n"; return @factors; } @@ -1678,7 +1666,228 @@ sub fermat_factor { @factors; } +sub ecm_factor { + my($n, $B1, $B2, $ncurves) = @_; + _validate_positive_integer($n); + + my @factors = _basic_factor($n); + return @factors if $n < 4; + + $ncurves = 10 unless defined $ncurves; + + if (!defined $B1) { + for my $mul (1, 10, 100, 1000, 10_000, 100_000, 1_000_000) { + $B1 = 100 * $mul; + $B2 = 1*$B1; + #warn "Trying ecm with $B1 / $B2\n"; + my @nf = ecm_factor($n, $B1, $B2, $ncurves); + if (scalar @nf > 1) { + push @factors, @nf; + return @factors; + } + } + push @factors, $n; + return @factors; + } + + $B2 = 1*$B1 unless defined $B2; + my $sqrt_b1 = int(sqrt($B1)+1); + + if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) { + eval { require Math::Prime::Util::ECAffinePoint; 1; } + or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; }; + } + + # With multiple curves, it's better to get all the primes at once. The + # downside is this can kill memory with a very large B1. + my @bprimes = @{ primes(2, $B1) }; + my $irandf = Math::Prime::Util::_get_rand_func(); + + foreach my $curve (1 .. $ncurves) { + my $a = $irandf->($n-1); + my $b = 1; + my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $n, 0, 1); + + foreach my $q (@bprimes) { + my $k = $q; + if ($k < $sqrt_b1) { + my $kmin = int($B1 / $q); + while ($k <= $kmin) { $k *= $q; } + } + $ECP->mul($k); + my $f = $ECP->f; + if ($f != 1) { + last if $f == $n; + push @factors, $f; + push @factors, int($n/$f); + croak "internal error in ecm" unless ($f * int($n/$f)) == $n; + warn "ECM found factors with B1 = $B1 in curve $curve\n"; + return @factors; + } + last if $ECP->is_infinity; + } + } + push @factors, $n; + @factors; +} + + +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 ($n,[$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 @fproof; + if (Math::Prime::Util::is_provable_prime($f, \@fproof) != 2) { + carp "could not prove primality of $n.\n"; + return (1, []); + } + push @fac_proofs, (scalar @fproof == 1) ? @fproof : [@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 ($n,[$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"; }; + } + $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; + 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) { + if (($B % $f) == 0) { + push @factors, $f; + do { $B /= $f; $A *= $f; } while (($B % $f) == 0); + } + } + my @nstack; + # nstack should only hold composites + if (Math::Prime::Util::is_prob_prime($B)) { + push @factors, $B; + $A *= $B; $B /= $B; # completely factored already + } else { + push @nstack, $B; + } + 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; + + 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) { + last if scalar @ftry >= 2; + @ftry = $sub->($m); + } + # If we couldn't find a factor, skip it. + next unless scalar @ftry > 1; + # Process each factor + foreach my $f (@ftry) { + croak "Invalid factoring" 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; + undef @uf{@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 = ($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(); + 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 @fproof; + if (Math::Prime::Util::is_provable_prime($f, \@fproof) != 2) { + carp "could not prove primality of $n.\n"; + return (1, []); + } + push @fac_proofs, (scalar @fproof == 1) ? @fproof : [@fproof]; + } + my @proof = ("$n", "n-1", [@fac_proofs], [@as]); + return (2, [@proof]); +} my $_const_euler = 0.57721566490153286060651209008240243104215933593992; -- 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