[libmath-prime-util-perl] 09/18: PP: add simple ECM factoring and BLS75 primality proof
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 b83566207a51be1f6a34f575c92ce179e6f80279
Author: Dana Jacobsen <dana at 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>dana at acm.orgE<gt>
+
+
+=head1 COPYRIGHT
+
+Copyright 2012-2013 by Dana Jacobsen E<lt>dana at 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(",", at 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(",", at 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
More information about the Pkg-perl-cvs-commits
mailing list