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 <d...@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>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/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

_______________________________________________
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

Reply via email to