[libmath-prime-util-perl] 07/18: Add primality certificates, elliptic curve start
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 592e6c3a106206e2579d3f6f3dfd6145fb6c3ddb
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Apr 12 20:14:21 2013 -0700
Add primality certificates, elliptic curve start
---
Changes | 6 +-
MANIFEST | 1 +
TODO | 10 ++
lib/Math/Prime/Util.pm | 316 +++++++++++++++++++++++++++++++++--
lib/Math/Prime/Util/EllipticCurve.pm | 220 ++++++++++++++++++++++++
t/81-bignum.t | 15 +-
6 files changed, 547 insertions(+), 21 deletions(-)
diff --git a/Changes b/Changes
index d2aaf4c..f89a628 100644
--- a/Changes
+++ b/Changes
@@ -4,7 +4,11 @@ Revision history for Perl extension Math::Prime::Util.
- Pure Perl p-1 factoring, Fermat factoring, and speedup for pbrent.
- - New verify_prime function that can check a primality proof.
+ - New functions:
+ prime_certificate produces a certificate of primality.
+ verify_prime checks a primality certificate.
+
+ - Math::Prime::Util::EllipticCurve module.
0.25 19 March 2013
diff --git a/MANIFEST b/MANIFEST
index 6ff97b3..bdd9012 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -4,6 +4,7 @@ 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
LICENSE
Makefile.PL
MANIFEST
diff --git a/TODO b/TODO
index c62ae60..d3db691 100644
--- a/TODO
+++ b/TODO
@@ -34,3 +34,13 @@
- Implement S2 calculation for LMO prime count.
- Add primality proof output to is_provable_prime.
+
+- Fixup EllipticCurve so we define points and a curve.
+
+- Change is_provable_prime from Lucas to BLS75.
+
+- Use EllipticCurve to make ecm_factor.
+
+- add recursive AGKM format.
+
+- Can we easily provide a certificate for Maurer primes?
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 9b25cea..0243ade 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -15,7 +15,8 @@ use base qw( Exporter );
our @EXPORT_OK =
qw( prime_get_config prime_set_config
prime_precalc prime_memfree
- is_prime is_prob_prime is_provable_prime verify_prime
+ is_prime is_prob_prime is_provable_prime
+ prime_certificate verify_prime
is_strong_pseudoprime is_strong_lucas_pseudoprime
is_aks_prime
miller_rabin
@@ -1600,17 +1601,44 @@ sub is_prob_prime {
}
sub is_provable_prime {
- my($n) = @_;
+ my($n, $ref_proof) = @_;
return 0 if defined $n && $n < 2;
_validate_positive_integer($n);
- # Shortcut some of the calls.
- return _XS_is_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
- return Math::Prime::Util::GMP::is_provable_prime($n) if $_HAVE_GMP
- && defined &Math::Prime::Util::GMP::is_provable_prime;
+ if (defined $ref_proof) {
+ croak "Second argument must be an array ref" if ref($ref_proof) ne 'ARRAY';
+ @$ref_proof = ();
+ }
- my $is_prob_prime = is_prob_prime($n);
- return $is_prob_prime unless $is_prob_prime == 1;
+ # 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);
+ }
+ # proof needed but MPU::GMP too old to give it.
+ }
+
+ my $is_prob_prime = is_prob_prime($n);
+ if ($is_prob_prime != 1) {
+ @$ref_proof = ($n) if defined $ref_proof && $is_prob_prime == 2;
+ return $is_prob_prime;
+ }
+ } else {
+ if ($n <= 10) {
+ if ($n==2||$n==3||$n==5||$n==7) {
+ @$ref_proof = ($n) if defined $ref_proof;
+ return 2;
+ }
+ return 0;
+ }
+ }
# 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
@@ -1626,31 +1654,67 @@ sub is_provable_prime {
}
my $nm1 = $n-1;
my @factors = factor($nm1);
- # Remember that we have to prove the primality of every factor.
- if ( (scalar grep { is_provable_prime($_) != 2 } @factors) > 0) {
- carp "could not prove primality of $n.\n";
- return 1;
+ # 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;
}
for (my $a = 2; $a < $nm1; $a++) {
my $ap = Math::BigInt->new("$a");
- # 1. a^(n-1) = 1 mod n.
- next if $ap->copy->bmodpow($nm1, $n) != 1;
- # 2. a^((n-1)/f) != 1 mod n for all f.
+ # 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;
}
+sub prime_certificate {
+ my($n) = @_;
+ return () if defined $n && $n < 2;
+ _validate_positive_integer($n);
+
+ my @cert;
+ my $is_prime = is_provable_prime($n, \@cert);
+ return () unless $is_prime == 2;
+ return @cert;
+}
+
sub verify_prime {
my @pdata = @_;
my $verbose = $_Config{'verbose'};
- croak "Parameter must be defined" if scalar @pdata == 0;
+ # 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]$/;
@@ -1662,7 +1726,7 @@ sub verify_prime {
eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
or do { croak "Cannot load Math::BigInt"; };
}
- $n = Math::BigInt->new("$n");
+ $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;
@@ -1689,8 +1753,65 @@ sub verify_prime {
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')) {
+ warn "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;
+ }
+ croak "Pratt error: n-1 not completely factored" unless $B == 1;
+
+ # 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
+ # BLS75 or generalized Pocklington
# http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
if (scalar @pdata != 2 || (ref($pdata[0]) ne 'ARRAY') || (ref($pdata[1]) ne 'ARRAY')) {
warn "verify_prime: incorrect n-1 format, must have factors and a values\n";
@@ -1771,6 +1892,77 @@ sub verify_prime {
return 1;
}
+ 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) {
+ warn "verify_prime: incorrect AGKM format\n";
+ return 0;
+ }
+ my ($ni, $a, $b, $m, $q, $P);
+ $q = $n;
+ foreach my $block (@pdata) {
+ if (ref($block) ne 'ARRAY' || scalar @$block != 6) {
+ warn "verify_prime: incorrect AGKM block format\n";
+ return 0;
+ }
+ if ($block->[0] != $q) {
+ warn "verify_prime: incorrect AGKM block format: block n != q\n";
+ return 0;
+ }
+ ($ni, $a, $b, $m, $q, $P) = @$block;
+ if (ref($P) ne 'ARRAY' || scalar @$P != 2) {
+ warn "verify_prime: incorrect AGKM block point format\n";
+ return 0;
+ }
+ $ni = $n->bzero->badd($ni) unless ref($ni) eq 'Math::BigInt';
+ if (Math::BigInt::bgcd($ni, 6) != 1) {
+ warn "verify_prime: AGKM block n '$ni' is divisible by 2 or 3\n";
+ return 0;
+ }
+ my $c = ($n-$n+4) * $a*$a*$a + ($n-$n+27)*$b*$b;
+ if (Math::BigInt::bgcd($c, $ni) != 1) {
+ warn "verify_prime: AGKM block gcd 4a^3+27b^2,n incorrect\n";
+ return 0;
+ }
+ if ($q <= $ni->copy->broot(4)->badd(1)->bpow(2)) {
+ 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"; };
+ }
+ 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]);
+ # 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)) ) {
+ warn "verify_prime: AGKM point does not multiply correctly.\n";
+ return 0;
+ }
+ }
+ # Check primality of last q using BPSW
+ return 0 unless verify_prime($q);
+
+ print "primality success: $n by A-K-G-M elliptic curve\n" if $verbose > 1;
+ return 1;
+ }
+
warn "verify_prime: Unknown method: '$method'.\n";
return 0;
}
@@ -2626,6 +2818,94 @@ be factored to the cube root of C<n>, which is more likely to succeed and will
usually take less time, but can still fail. Hence you should always test that
the result is C<2> to ensure the prime is proven.
+An optional second argument may be given, which must be an array reference,
+which will be filled in with a primality certificate. Normally this is used
+via the function L</prime_certificate> below.
+
+
+=head2 prime_certificate
+
+ my @cert = prime_certificate($n);
+ say verify_prime(@cert) ? "proven prime" : "not prime";
+
+Given a positive integer C<n> as input, returns either an empty array (we could
+not prove C<n> prime) or an array representing a certificate of primality.
+This may be examined or given to L</verify_prime> for verification. The latter
+function contains the description of the format.
+
+
+=head2 verify_prime
+
+ my @cert = prime_certificate($n);
+ say verify_prime(@cert) ? "proven prime" : "not prime";
+
+Given an array representing a certificate of primality, returns either 0 (not
+verified), or 1 (verified). The computations are all done using pure Perl
+Math::BigInt and should not be time consuming (the Pari or GMP backends will
+help with large inputs).
+
+A certificate is an array holding an C<n-cert>. An C<n-cert> is one of:
+
+ n
+ implies n,"BPSW"
+
+ n,"BPSW"
+ the number n is small enough to be proven with BPSW. This
+ currently means smaller than 2^64.
+
+ n,"Pratt",[n-cert, ...],a
+ A Pratt certificate. We are given n, the method "Pratt" or "Lucas",
+ a list of n-certs that indicate all the unique factors of n-1, and
+ an 'a' value to be used in the Lucas primality test.
+ The certificate passes if:
+ 1 all factor n-certs can be verified
+ 2 all n-certs are factors of n-1 and none are missing
+ 3 a is coprime to n
+ 4 a^(n-1) = 1 mod n
+ 5 a^((n-1)/f) != 1 mod n for each factor
+
+ n,"n-1",[n-cert, ...],[a,...]
+ An n-1 certificate suitable for the generalized Pocklington or the
+ BLS75 (Brillhart-Lehmer-Selfridge 1975, theorem 5) test. The proof
+ is performed using BLS75 theorem 5 which requires n-1 to be factored
+ up to (n/2)^1/3. If n-1 is factored to more than sqrt(n), then the
+ conditions are identical to the generalized Pocklington test.
+ The certificate passes if:
+ 1 all factor n-certs can be verified
+ 2 all factor n-certs are factors of n-1
+ 3 there must be a corresponding 'a' for each factor n-cert
+ 4 given A (the factored part of n-1), B = (n-1)/A (the unfactored
+ part), s = int(B/(2A)), r = B-s*2A:
+ - n < (A+1)(2*A*A+(r-a)A+a) [ n-1 factored to (n/2)^1/3 ]
+ - s = 0 or r*r-8s not a perfect square
+ - A and B are coprime
+ 5 for each pair (f,a) representing a factor n-cert and its 'a':
+ - a^(n-1) = 1 mod n
+ - gcd( a^((n-1)/f)-1, n ) = 1
+
+ n,"AGKM",[ec-block],[ec-block],...
+ An Elliptic Curve certificate. We are given n, the method "AGKM"
+ or "ECPP", and a one or more 6-element blocks representing a standard
+ ECPP or Atkin-Goldwasser-Kilian-Morain certificate. The format of
+ this n-cert is non-recursive so it can be easily used for similar
+ programs such as Sage and GMP-ECPP.
+ Every ec-block has 6 elements:
+ N the N value this block proves prime if q is prime
+ a value describing the elliptic curve to be used
+ b value describing the elliptic curve to be used
+ m order of the curve
+ q a probable prime > (N^1/4+1)^2
+ P a point [x,y] on the curve (affine coordinates)
+ The certificate passes if:
+ - the final q can be proved with BPSW.
+ - for each block:
+ - N is the same as the preceeding block's q
+ - N is not divisible by 2 or 3
+ - gcd( 4a^3 + 27b^2, N ) == 1;
+ - q > (N^1/4+1)^2
+ - U = (m/q)P is not the point at infinity
+ - V = qU is the point at infinity
+
=head2 is_aks_prime
diff --git a/lib/Math/Prime/Util/EllipticCurve.pm b/lib/Math/Prime/Util/EllipticCurve.pm
new file mode 100644
index 0000000..fff5fef
--- /dev/null
+++ b/lib/Math/Prime/Util/EllipticCurve.pm
@@ -0,0 +1,220 @@
+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;
+
+ while ($v && $k > 0) {
+ if ( ($k % 2) != 0) {
+ $k--;
+ my $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);
+ $v = ($d == 1 || $d == $n);
+ if ($v) { ($x,$y) = $self->double_a($x,$y); }
+ }
+ }
+ return ($Bx, $By);
+}
+
+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/t/81-bignum.t b/t/81-bignum.t
index 0d2b838..822291e 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -76,7 +76,7 @@ plan tests => 0
+ scalar(keys %allfactors)
+ 5 # moebius, euler_phi, jordan totient
+ 15 # random primes
- + 3 # verify_prime
+ + 5 # verify_prime
+ 0;
# Using GMP makes these tests run about 2x faster on some machines
@@ -294,6 +294,15 @@ sub check_pcbounds {
# Provable primes
{
my @proof;
+
+ @proof = (20907001, "Pratt", [ 2,
+ 3,
+ 5,
+ [23,"Pratt",[2,[11,"Pratt",[2,5],2]],5],
+ [101, "Pratt", [2, 5], 2],
+ ], 14 );
+ ok( verify_prime(@proof), "simple Lucas/Pratt proof verified" );
+
@proof = ('3364125245431456304736426076174232972735419017865223025179282077503701', 'n-1',
[2,5,127, ['28432789963853652887491983185920687231739655787', 'n-1',
[ 2,3,163,650933, [ '44662634059309451871488121651101494489', 'n-1',
@@ -304,7 +313,7 @@ sub check_pcbounds {
],
'9316417838190714313' ],
[ 2, 2, 2, 2, 2 ]);
- ok( verify_prime(@proof), "simple proof verified" );
+ ok( verify_prime(@proof), "simple n-1 proof verified" );
@proof = ('6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151', 'n-1',
[ 2,3,5,11,17,31,41,53,131,157,521,1613,61681,8191,42641,858001,51481, '7623851', '308761441' ],
@@ -320,4 +329,6 @@ sub check_pcbounds {
[ 3,5,3,2,3,3,3,3 ] );
ok( verify_prime(@proof), "2**607-1 primality proof verified" );
+ @proof = ('677826928624294778921',"AKGM", [677826928624294778921, 404277700094248015180, 599134911995823048257, 677826928656744857936, 104088901820753203, [2293544533, 356794037129589115041]], [104088901820753203, 0, 73704321689372825, 104088902465395836, 1112795797, [3380482019, 53320146243107032]], [1112795797, 0, 638297481, 1112860899, 39019, [166385704, 356512285]]);
+ ok( verify_prime(@proof), "ECPP primality proof of 677826928624294778921 verified" );
}
--
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