[libmath-prime-util-perl] 11/18: Projective ECM for pure Perl
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 1b374a4cf742a12c587028dc46b4b3e57f97e120
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Apr 15 18:50:21 2013 -0700
Projective ECM for pure Perl
---
Changes | 2 +-
MANIFEST | 2 +-
factor.c | 4 +-
lib/Math/Prime/Util/ECAffinePoint.pm | 8 +-
lib/Math/Prime/Util/ECProjectivePoint.pm | 284 +++++++++++++++++++++++++++++++
lib/Math/Prime/Util/EllipticCurve.pm | 221 ------------------------
lib/Math/Prime/Util/PP.pm | 269 ++++++++++++++++++-----------
7 files changed, 466 insertions(+), 324 deletions(-)
diff --git a/Changes b/Changes
index 320f7e3..1425355 100644
--- a/Changes
+++ b/Changes
@@ -6,7 +6,7 @@ Revision history for Perl extension Math::Prime::Util.
- real p-1 -- much faster and more effective
- Fermat (no better than HOLF)
- speedup for pbrent
- - simple affine single stage ECM
+ - simple ECM
- redo factoring mix
- New functions:
diff --git a/MANIFEST b/MANIFEST
index 416fbfa..10a9f2b 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -4,8 +4,8 @@ lib/Math/Prime/Util/MemFree.pm
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
+lib/Math/Prime/Util/ECProjectivePoint.pm
LICENSE
Makefile.PL
MANIFEST
diff --git a/factor.c b/factor.c
index db91e4c..6d0962a 100644
--- a/factor.c
+++ b/factor.c
@@ -8,6 +8,7 @@
#include "util.h"
#include "sieve.h"
#include "mulmod.h"
+#include "cache.h"
/*
* You need to remember to use UV for unsigned and IV for signed types that
@@ -623,7 +624,8 @@ int prho_factor(UV n, UV *factors, UV rounds)
/* Pollard's P-1 */
int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
{
- UV q, f;
+ UV f;
+ UV q = 2;
UV a = 2;
UV savea = 2;
UV saveq = 2;
diff --git a/lib/Math/Prime/Util/ECAffinePoint.pm b/lib/Math/Prime/Util/ECAffinePoint.pm
index 81c8d06..f3c9d29 100644
--- a/lib/Math/Prime/Util/ECAffinePoint.pm
+++ b/lib/Math/Prime/Util/ECAffinePoint.pm
@@ -9,8 +9,8 @@ if (!defined $Math::BigInt::VERSION) {
}
BEGIN {
- $Math::Prime::Util::EllipticCurve::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::EllipticCurve::VERSION = '0.26';
+ $Math::Prime::Util::ECAffinePoint::AUTHORITY = 'cpan:DANAJ';
+ $Math::Prime::Util::ECAffinePoint::VERSION = '0.26';
}
# Pure perl (with Math::BigInt) manipulation of Elliptic Curves
@@ -180,9 +180,9 @@ To write.
=head2 new
- $point = Math::Prime::Util::EllipticCurve->new(a, b);
+ $point = Math::Prime::Util::ECAffinePoint->new(a, b, n, x, y);
-Returns a new curve defined by a and b.
+Returns a new point at C<(x,y,1)> on the curve C<(a,b,n)>.
=head2 a
diff --git a/lib/Math/Prime/Util/ECProjectivePoint.pm b/lib/Math/Prime/Util/ECProjectivePoint.pm
new file mode 100644
index 0000000..5af473f
--- /dev/null
+++ b/lib/Math/Prime/Util/ECProjectivePoint.pm
@@ -0,0 +1,284 @@
+package Math::Prime::Util::ECProjectivePoint;
+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::ECProjectivePoint::AUTHORITY = 'cpan:DANAJ';
+ $Math::Prime::Util::ECProjectivePoint::VERSION = '0.26';
+}
+
+# Pure perl (with Math::BigInt) manipulation of Elliptic Curves
+# in projective coordinates.
+
+sub new {
+ my ($class, $a, $b, $n, $x, $z) = @_;
+ $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';
+ $z = Math::BigInt->new("$z") unless ref($z) 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,
+ z => $z,
+ f => $n-$n+1,
+ };
+
+ bless $self, $class;
+ return $self;
+}
+
+sub _add {
+ my ($self, $x2, $z2, $x1, $z1, $xin, $n) = @_;
+
+ my $u = ( ($x2 - $z2) * ($x1 + $z1) ) % $n;
+ my $v = ( ($x2 + $z2) * ($x1 - $z1) ) % $n;
+
+ my $puv = $u + $v;
+ my $muv = $u - $v;
+
+ return ( ($puv*$puv) % $n, ($muv*$muv * $xin) % $n );
+}
+
+sub _add3 {
+ my ($self, $x1, $z1, $x2, $z2, $xin, $zin, $n) = @_;
+
+ my $u = (($x2 - $z2) * ($x1 + $z1) ) % $n;
+ my $v = (($x2 + $z2) * ($x1 - $z1) ) % $n;
+
+ my $upv2 = (($u+$v) * ($u+$v)) % $n;
+ my $umv2 = (($u-$v) * ($u-$v)) % $n;
+
+ return ( ($upv2*$zin) % $n, ($umv2*$xin) % $n );
+}
+
+sub _double {
+ my ($self, $x, $z, $n) = @_;
+
+ my $u = $x + $z; $u = ($u * $u) % $n;
+ my $v = $x - $z; $v = ($v * $v) % $n;
+ my $w = $u - $v;
+
+ return ( ($u*$v)%$n , ($w*($v+$w*$self->{'b'}))%$n );
+}
+
+sub mul {
+ my ($self, $k) = @_;
+ my $x = $self->{'x'};
+ my $z = $self->{'z'};
+ my $n = $self->{'n'};
+
+ my ($x1, $x2, $z1, $z2);
+
+ my $r = --$k;
+ my $l = -1;
+ while ($r != 1) { $r >>= 1; $l++ }
+ if ($k & (1 << $l)) {
+ ($x2, $z2) = $self->_double($x, $z, $n);
+ ($x1, $z1) = $self->_add3($x2, $z2, $x, $z, $x, $z, $n);
+ ($x2, $z2) = $self->_double($x2, $z2, $n);
+ } else {
+ ($x1, $z1) = $self->_double($x, $z, $n);
+ ($x2, $z2) = $self->_add3($x, $z, $x1, $z1, $x, $z, $n);
+ }
+ $l--;
+ while ($l >= 1) {
+ if ($k & (1 << $l)) {
+ ($x1, $z1) = $self->_add3($x1, $z1, $x2, $z2, $x, $z, $n);
+ ($x2, $z2) = $self->_double($x2, $z2, $n);
+ } else {
+ ($x2, $z2) = $self->_add3($x2, $z2, $x1, $z1, $x, $z, $n);
+ ($x1, $z1) = $self->_double($x1, $z1, $n);
+ }
+ $l--;
+ }
+ if ($k & 1) {
+ ($x, $z) = $self->_double($x2, $z2, $n);
+ } else {
+ ($x, $z) = $self->_add3($x2, $z2, $x1, $z1, $x, $z, $n);
+ }
+
+ $self->{'x'} = $x;
+ $self->{'z'} = $z;
+ return $self;
+}
+
+sub add {
+ my ($self, $other) = @_;
+ croak "add takes a EC point"
+ unless ref($other) eq 'Math::Prime::Util::ECProjectivePoint';
+ 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->{'z'}) = $self->_add3($self->{'x'}, $self->{'z'},
+ $other->{'x'}, $other->{'z'},
+ $self->{'x'}, $self->{'z'},
+ $self->{'n'});
+ return $self;
+}
+
+sub double {
+ my ($self) = @_;
+ ($self->{'x'}, $self->{'z'}) = $self->_double($self->{'x'}, $self->{'z'}, $self->{'n'});
+ return $self;
+}
+
+sub _extended_gcd {
+ my ($a, $b) = @_;
+ my $zero = $a-$a;
+ my ($x, $lastx, $y, $lasty) = ($zero, $zero+1, $zero+1, $zero);
+ while ($b != 0) {
+ my $q = int($a/$b);
+ ($a, $b) = ($b, $a % $b);
+ ($x, $lastx) = ($lastx - $q*$x, $x);
+ ($y, $lasty) = ($lasty - $q*$y, $y);
+ }
+ return ($a, $lastx, $lasty);
+}
+
+sub normalize {
+ my ($self) = @_;
+ my $n = $self->{'n'};
+ my $z = $self->{'z'};
+ #my ($f, $u, undef) = _extended_gcd( $z, $n );
+ my $f = Math::BigInt::bgcd( $z, $n );
+ my $u = $z->copy->bmodinv($n);
+ $self->{'x'} = ( $self->{'x'} * $u ) % $n;
+ $self->{'z'} = $n-$n+1;
+ $self->{'f'} = ($f * $self->{'f'}) % $n;
+ return $self;
+}
+
+sub a { return shift->{'a'}; }
+sub b { return shift->{'b'}; }
+sub n { return shift->{'n'}; }
+sub x { return shift->{'x'}; }
+sub z { return shift->{'z'}; }
+sub f { return shift->{'f'}; }
+
+sub is_infinity {
+ my $self = shift;
+ return ($self->{'x'}->is_zero() && $self->{'z'}->is_one());
+}
+
+sub copy {
+ my $self = shift;
+ return Math::Prime::Util::ECProjectivePoint->new(
+ $self->{'a'}, $self->{'b'}, $self->{'n'}, $self->{'x'}, $self->{'z'});
+}
+
+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
deleted file mode 100644
index d04601b..0000000
--- a/lib/Math/Prime/Util/EllipticCurve.pm
+++ /dev/null
@@ -1,221 +0,0 @@
-package Math::Prime::Util::EllipticCurve;
-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 projective coordinates. Should be split into a point class.
-
-sub new {
- my ($class, $a, $b, $n) = @_;
- $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';
-
- my $self = {
- a => $a,
- b => $b,
- n => $n,
- };
-
- bless $self, $class;
- return $self;
-}
-
-sub double_p {
- my ($self, $x, $z) = @_;
- my $n = $self->{'n'};
-
- my $u = ( ($x+$z) * ($x+$z) ) % $n;
- my $v = ( ($x-$z) * ($x-$z) ) % $n;
- my $w = $u - $v;
-
- return ( ($u*$v)%$n , ($w*($v+$w*$self->{'b'}))%$n );
-}
-
-sub add3_p {
- my ($self, $x1, $z1, $x2, $z2, $xin, $zin) = @_;
- my $n = $self->{'n'};
-
- my $u = (($x2 - $z2) * ($x1 + $z1) ) % $n;
- my $v = (($x2 + $z2) * ($x1 - $z1) ) % $n;
-
- my $upv2 = (($u+$v) * ($u+$v)) % $n;
- my $umv2 = (($u-$v) * ($u-$v)) % $n;
-
- return ( ($upv2*$zin) % $n, ($umv2*$xin) % $n );
-}
-
-sub mul_p {
- my ($self, $k, $x, $z) = @_;
- $x = Math::BigInt->new("$x") unless ref($x) eq 'Math::BigInt';
- $z = Math::BigInt->new("$z") unless ref($z) eq 'Math::BigInt';
- my ($x1, $x2, $z1, $z2);
-
- my $r = --$k;
- my $l = -1;
- while ($r != 1) { $r >>= 1; $l++ }
- if ($k & (1 << $l)) {
- ($x2, $z2) = $self->double_p($x, $z);
- ($x1, $z1) = $self->add3_p($x2, $z2, $x, $z, $x, $z);
- ($x2, $z2) = $self->double_p($x2, $z2);
- } else {
- ($x1, $z1) = $self->double_p($x, $z);
- ($x2, $z2) = $self->add3_p($x, $z, $x1, $z1, $x, $z);
- }
- $l--;
- while ($l >= 1) {
- if ($k & (1 << $l)) {
- ($x1, $z1) = $self->add3_p($x1, $z1, $x2, $z2, $x, $z);
- ($x2, $z2) = $self->double_p($x2, $z2);
- } else {
- ($x2, $z2) = $self->add3_p($x2, $z2, $x1, $z1, $x, $z);
- ($x1, $z1) = $self->double_p($x1, $z1);
- }
- $l--;
- }
- if ($k & 1) {
- ($x, $z) = $self->double_p($x2, $z2);
- } else {
- ($x, $z) = $self->add3_p($x2, $z2, $x1, $z1, $x, $z);
- }
- return ($x, $z);
-}
-
-sub add_a {
- my ($self, $P1x, $P1y, $P2x, $P2y) = @_;
- my $n = $self->{'n'};
-
- if ($P1x == $P2x) {
- my $t = ($P1y + $P2y) % $n;
- return (0, 1) 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_a {
- 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_a {
- my ($self, $k, $x, $y) = @_;
- $x = Math::BigInt->new("$x") unless ref($x) eq 'Math::BigInt';
- $y = Math::BigInt->new("$y") unless ref($y) eq 'Math::BigInt';
- my $n = $self->{'n'};
-
- 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--;
- $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;
- $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, $d);
-}
-
-1;
-
-__END__
-
-
-# ABSTRACT: Elliptic curve operations
-
-=pod
-
-=encoding utf8
-
-
-=head1 NAME
-
-Math::Prime::Util::EllipticCurve - Elliptic curve operations
-
-
-=head1 VERSION
-
-Version 0.26
-
-
-=head1 SYNOPSIS
-
-Todo.
-
- Todo
- # TO DO
- To do
-
-=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.
-
-
-=head1 SEE ALSO
-
-L<Math::EllipticCurve::Prime>
-
-
-=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/PP.pm b/lib/Math/Prime/Util/PP.pm
index acf1096..b94e684 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1146,13 +1146,13 @@ my @_fsublist = (
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 { ecm_factor (shift, 1_000, 5_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 { ecm_factor (shift, 10_000, 50_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 { ecm_factor (shift, 100_000, 400_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; },
@@ -1213,6 +1213,18 @@ sub factor {
sort {$a<=>$b} @factors;
}
+sub _found_factor {
+ my($f, $n, $what, @factors) = @_;
+ if ($f == 1 || $f == $n) {
+ push @factors, $n;
+ } else {
+ push @factors, $f;
+ push @factors, int($n/$f);
+ croak "internal error in $what" unless ($f * int($n/$f)) == $n;
+ }
+ @factors;
+}
+
# TODO:
sub squfof_factor { trial_factor(@_) }
@@ -1242,11 +1254,7 @@ sub prho_factor {
if ($f == $n) {
last if $inloop++; # We've been here before
} elsif ($f != 1) {
- my $f2 = $n->copy->bdiv($f)->as_int;
- push @factors, $f;
- push @factors, $f2;
- croak "internal error in prho" unless ($f * $f2) == $n;
- return @factors;
+ return _found_factor($f, $n, "prho", @factors);
}
}
@@ -1260,10 +1268,7 @@ sub prho_factor {
if ($f == $n) {
last if $inloop++; # We've been here before
} elsif ($f != 1) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in prho" unless ($f * int($n/$f)) == $n;
- return @factors;
+ return _found_factor($f, $n, "prho", @factors);
}
}
@@ -1277,10 +1282,7 @@ sub prho_factor {
if ($f == $n) {
last if $inloop++; # We've been here before
} elsif ($f != 1) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in prho" unless ($f * int($n/$f)) == $n;
- return @factors;
+ return _found_factor($f, $n, "prho", @factors);
}
}
@@ -1342,11 +1344,7 @@ sub pbrent_factor {
} while ($f != 1 && $r-- != 0);
last if $f == 1 || $f == $n;
}
- my $f2 = $n->copy->bdiv($f)->as_int;
- push @factors, $f;
- push @factors, $f2;
- croak "internal error in pbrent" unless ($f * $f2) == $n;
- return @factors;
+ return _found_factor($f, $n, "pbrent", @factors);
}
} elsif ($n < $_half_word) {
@@ -1354,12 +1352,7 @@ sub pbrent_factor {
for my $i (1 .. $rounds) {
$Xi = ($Xi * $Xi + $a) % $n;
my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
- if ( ($f != 1) && ($f != $n) ) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
- return @factors;
- }
+ return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
$Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
}
@@ -1370,12 +1363,7 @@ sub pbrent_factor {
$Xi = _mulmod($Xi, $Xi, $n);
$Xi = (($n-$Xi) > $a) ? $Xi+$a : $Xi+$a-$n;
my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
- if ( ($f != 1) && ($f != $n) ) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
- return @factors;
- }
+ return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
$Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
}
@@ -1413,13 +1401,7 @@ sub pminus1_factor {
$a = _powmod($a, $k, $n);
if ($a == 0) { push @factors, $n; return @factors; }
my $f = _gcd_ui( $a-1, $n );
- if ($f == $n) { push @factors, $n; return @factors; }
- if ($f != 1) {
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in pminus1" unless ($f * int($n/$f)) == $n;
- return @factors;
- }
+ return _found_factor($f, $n, "pminus1", @factors) if $f != 1;
}
last if $pc_end >= $B1;
$pc_beg = $pc_end+1;
@@ -1540,12 +1522,7 @@ sub pminus1_factor {
}
$f = Math::BigInt::bgcd( $b, $n );
}
- if ($f != 1 && $f != $n) {
- push @factors, $f, $n/$f;
- return @factors;
- }
- push @factors, $n;
- @factors;
+ return _found_factor($f, $n, "pminus1", @factors);
}
sub holf_factor {
@@ -1565,11 +1542,7 @@ sub holf_factor {
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;
+ return _found_factor($f, $n, "HOLF", @factors);
}
$s->binc;
my $m = ($s * $s) - $ni;
@@ -1579,12 +1552,7 @@ sub holf_factor {
my $f = $m->copy->bsqrt->bfloor->as_int;
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
- push @factors, $f;
- 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;
+ return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
}
} else {
for my $i ($startrounds .. $rounds) {
@@ -1597,12 +1565,7 @@ sub holf_factor {
my $f = int(sqrt($m));
next unless $f*$f == $m;
$f = _gcd_ui( ($s > $f) ? $s - $f : $f - $s, $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;
- # print "HOLF found factors in $i rounds\n";
- return @factors;
+ return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
}
}
push @factors, $n;
@@ -1619,7 +1582,7 @@ sub fermat_factor {
if ( ref($n) eq 'Math::BigInt' ) {
my $a = $n->copy->bsqrt->bfloor->as_int;
- if ($a*$a == $n) { push @factors, $a, $a; return @factors; }
+ return _found_factor($a, $n, "Fermat", @factors) if $a*$a == $n;
$a++;
my $b2 = $a*$a - $n;
my $lasta = $a + $rounds;
@@ -1628,12 +1591,8 @@ sub fermat_factor {
if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
my $s = $b2->copy->bsqrt->bfloor->as_int;
if ($s*$s == $b2) {
- my $f = $a - $s;
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in Fermat" unless ($f * int($n/$f)) == $n;
- #warn "Fermat found factors in ",$a-($lasta-$rounds)+1," rounds\n";
- return @factors;
+ my $i = $a-($lasta-$rounds)+1;
+ return _found_factor($a - $s, $n, "Fermat ($i rounds)", @factors);
}
}
$a++;
@@ -1641,7 +1600,7 @@ sub fermat_factor {
}
} else {
my $a = int(sqrt($n));
- if ($a*$a == $n) { push @factors, $a, $a; return @factors; }
+ return _found_factor($a, $n, "Fermat", @factors) if $a*$a == $n;
$a++;
my $b2 = $a*$a - $n;
my $lasta = $a + $rounds;
@@ -1650,12 +1609,8 @@ sub fermat_factor {
if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
my $s = int(sqrt($b2));
if ($s*$s == $b2) {
- my $f = $a - $s;
- push @factors, $f;
- push @factors, int($n/$f);
- croak "internal error in Fermat" unless ($f * int($n/$f)) == $n;
- #warn "Fermat found factors in ",$a-($lasta-$rounds)+1," rounds\n";
- return @factors;
+ my $i = $a-($lasta-$rounds)+1;
+ return _found_factor($a - $s, $n, "Fermat ($i rounds)", @factors);
}
}
$a++;
@@ -1666,6 +1621,7 @@ sub fermat_factor {
@factors;
}
+
sub ecm_factor {
my($n, $B1, $B2, $ncurves) = @_;
_validate_positive_integer($n);
@@ -1678,7 +1634,7 @@ sub ecm_factor {
if (!defined $B1) {
for my $mul (1, 10, 100, 1000, 10_000, 100_000, 1_000_000) {
$B1 = 100 * $mul;
- $B2 = 1*$B1;
+ $B2 = 10*$B1;
#warn "Trying ecm with $B1 / $B2\n";
my @nf = ecm_factor($n, $B1, $B2, $ncurves);
if (scalar @nf > 1) {
@@ -1693,21 +1649,71 @@ sub ecm_factor {
$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"; };
+ # Affine code. About 3x slower than the projective, and no stage 2.
+ #
+ #if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) {
+ # eval { require Math::Prime::Util::ECAffinePoint; 1; }
+ # or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; };
+ #}
+ #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;
+ # warn "ECM found factors with B1 = $B1 in curve $curve\n";
+ # return _found_factor($f, $n, "ECM B1=$B1 curve $curve", @factors);
+ # }
+ # last if $ECP->is_infinity;
+ # }
+ #}
+
+ if (!defined $Math::Prime::Util::ECProjectivePoint::VERSION) {
+ eval { require Math::Prime::Util::ECProjectivePoint; 1; }
+ or do { croak "Cannot load Math::Prime::Util::ECProjectivePoint"; };
}
- # 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) };
+ # 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(3, $B1) };
+ my @b2primes = ($B2 > $B1) ? @{primes($B1+1, $B2)} : ();
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);
-
+ my $sigma;
+ do { $sigma = $irandf->($n-1) } while ($sigma <= 5);
+ my ($u, $v) = ( ($sigma*$sigma - 5) % $n, (4 * $sigma) % $n );
+ my ($x, $z) = ( ($u*$u*$u) % $n, ($v*$v*$v) % $n );
+ my $b = (4 * $x * $v) % $n;
+ my $a = ( (($v-$u)**3) * (3*$u + $v) ) % $n;
+ my $f = Math::BigInt::bgcd( $b, $n );
+ $f = Math::BigInt::bgcd( $z, $n ) if $f == 1;
+ next if $f == $n;
+ return _found_factor($f,$n, "ECM B1=$B1 curve $curve", @factors) if $f != 1;
+ $u = $b->copy->bmodinv($n);
+ $a = (($a*$u) - 2) % $n;
+ $b = $a+2;
+ $b = ( (($b % 2) == 0) ? $b : $b+$n ) >> 1;
+ $b = ( (($b % 2) == 0) ? $b : $b+$n ) >> 1;
+ #$u = $z->copy->bmodinv($n);
+ #$x = ($x * $u) % $n;
+ #$z = $n-$n+1;
+
+ my $ECP = Math::Prime::Util::ECProjectivePoint->new($a, $b, $n, $x, $z);
+ my $fm = $n-$n+1;
+ my $i = 15;
+
+ for (my $q = 2; $q < $B1; $q *= 2) { $ECP->double(); }
foreach my $q (@bprimes) {
my $k = $q;
if ($k < $sqrt_b1) {
@@ -1715,17 +1721,88 @@ sub ecm_factor {
while ($k <= $kmin) { $k *= $q; }
}
$ECP->mul($k);
- my $f = $ECP->f;
+ $fm = ($fm * $ECP->x() ) % $n;
+ if ($i++ % 32 == 0) {
+ $f = Math::BigInt::bgcd($fm, $n);
+ last if $f != 1;
+ }
+ }
+ $f = Math::BigInt::bgcd($fm, $n);
+ next if $f == $n;
+
+ if ($f == 1 && $B2 > $B1) { # BEGIN STAGE 2
+ my $D = int(sqrt($B2/2)); $D++ if $D % 2;
+ my $one = $n - $n + 1;
+ my $g = $one;
+
+ my $S2P = $ECP->copy->normalize;
+ $f = $S2P->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;
+ next if $f == $n;
+ #warn "ECM S2 normalize f=$f\n" if $f != 1;
+ return _found_factor($f, $n, "ECM S2 B1=$B1 curve $curve");
+ }
+ my $S2x = $S2P->x;
+ my @nqx = ($n-$n, $S2x);
+
+ foreach my $i (2 .. 2*$D) {
+ my($x2, $z2);
+ if ($i % 2) {
+ ($x2, $z2) = $S2P->_add($nqx[($i+1)/2], $one, $nqx[($i-1)/2], $one, $S2x, $n);
+ } else {
+ ($x2, $z2) = $S2P->_double($nqx[$i/2], $one, $n);
+ }
+ $nqx[$i] = $x2;
+ #($f, $u, undef) = _extended_gcd($z2, $n);
+ $f = Math::BigInt::bgcd( $z2, $n );
+ last if $f != 1;
+ $u = $z2->copy->bmodinv($n);
+ $nqx[$i] = ($x2 * $u) % $n;
+ }
+ if ($f != 1) {
+ next if $f == $n;
+ #warn "ECM S2 1: B1 $B1 B2 $B2 curve $curve f=$f\n";
+ return _found_factor($f, $n, "ECM S2 B1=$B1 curve $curve", @factors);
+ }
+
+ $x = $nqx[2*$D-1];
+ my $m = 1;
+ while ($m < ($B2+$D)) {
+ if ($m != 1) {
+ my $oldx = $S2x;
+ my ($x1, $z1) = $S2P->_add($S2x, $one, $nqx[2*$D], $one, $x, $n);
+ $f = Math::BigInt::bgcd( $z1, $n );
+ last if $f != 1;
+ $u = $z1->copy->bmodinv($n);
+ $S2x = ($x1 * $u) % $n;
+ $x = $oldx;
+ last if $f != 1;
+ }
+ if ($m+$D > $B1) {
+ my @p = grep { $_ >= $m-$D && $_ <= $m+$D } @b2primes;
+ foreach my $i (@p) {
+ last if $i >= $m;
+ $g = ($g * ($S2x - $nqx[$m+$D-$i])) % $n;
+ }
+ foreach my $i (@p) {
+ next unless $i > $m;
+ next if is_prime($m+$m-$i);
+ $g = ($g * ($S2x - $nqx[$i-$m])) % $n;
+ }
+ $f = Math::BigInt::bgcd($g, $n);
+ #warn "ECM S2 3: found $f in stage 2\n" if $f != 1;
+ last if $f != 1;
+ }
+ $m += 2*$D;
}
- last if $ECP->is_infinity;
+ } # END STAGE 2
+
+ next if $f == $n;
+ if ($f != 1) {
+ #warn "ECM found factors with B1 = $B1 in curve $curve\n";
+ return _found_factor($f, $n, "ECM B1=$B1 curve $curve", @factors);
}
+ # end of curve loop
}
push @factors, $n;
@factors;
--
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