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 aa43b9b989ab0d56daad9eeb0ba4e84d768f3861
Author: Dana Jacobsen <d...@acm.org>
Date:   Wed Apr 17 19:30:40 2013 -0700

    Primality proof updates
---
 lib/Math/Prime/Util.pm                   | 177 ++++++++++++++++++++-----------
 lib/Math/Prime/Util/ECProjectivePoint.pm |  25 +++--
 lib/Math/Prime/Util/PP.pm                |  58 +++++++---
 3 files changed, 179 insertions(+), 81 deletions(-)

diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index e023cb2..e7e4d47 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -15,7 +15,7 @@ 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
+      is_prime is_prob_prime is_provable_prime is_provable_prime_with_cert
       prime_certificate verify_prime
       is_strong_pseudoprime is_strong_lucas_pseudoprime
       is_aks_prime
@@ -26,7 +26,7 @@ our @EXPORT_OK =
       prime_count_lower prime_count_upper prime_count_approx
       nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
       random_prime random_ndigit_prime random_nbit_prime
-      random_strong_prime random_maurer_prime
+      random_strong_prime random_maurer_prime random_maurer_prime_with_cert
       primorial pn_primorial consecutive_integer_lcm
       factor all_factors
       moebius mertens euler_phi jordan_totient exp_mangoldt
@@ -834,10 +834,21 @@ sub primes {
   }
 
   sub random_maurer_prime {
+    my ($n, $cert) = random_maurer_prime_with_cert(@_);
+    croak "maurer prime $n failed certificate verification!"
+          unless verify_prime(@$cert);
+    return $n;
+  }
+
+  sub random_maurer_prime_with_cert {
     my($k) = @_;
     _validate_positive_integer($k, 2);
+    my @cert;
     if ($] < 5.008 && $_Config{'maxbits'} > 32) {
-      return random_nbit_prime($k) if $k <= 49;
+      if ($k <= 49) {
+        my $n = random_nbit_prime($k);
+        return ($n, [$n]);
+      }
       croak "Random Maurer not supported on old Perl";
     }
 
@@ -846,7 +857,12 @@ sub primes {
     # returning 2.  This should be reasonably fast to ~128 bits with MPU::GMP.
     my $p0 = $_Config{'maxbits'};
 
-    return random_nbit_prime($k) if $k <= $p0;
+    if ($k <= $p0) {
+      my $n = random_nbit_prime($k);
+      my ($isp, $cert) = is_provable_prime_with_cert($n);
+      croak "small nbit prime could not be proven" if $isp != 2;
+      return ($n, $cert);
+    }
 
     if (!defined $Math::BigInt::VERSION) {
       eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
@@ -877,7 +893,7 @@ sub primes {
     }
 
     # I've seen +0, +1, and +2 here.  Maurer uses +0.  Menezes uses +1.
-    my $q = random_maurer_prime( ($r * $k)->bfloor + 1 );
+    my ($q, $certref) = random_maurer_prime_with_cert( ($r * $k)->bfloor + 1 );
     $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt';
     my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int();
     print "B = $B  r = $r  k = $k  q = $q  I = $I\n" if $verbose && $verbose 
!= 3;
@@ -955,7 +971,15 @@ sub primes {
         #  2) make history by finding the first known BPSW pseudo-prime
         croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n);
 
-        return $n;
+        # Build up cert, knowing n-1 = 2*q*R, q > sqrt(n).
+        # We'll need to find the right a value for the factor 2.
+        foreach my $f2a (2 .. 200) {
+          $a = Math::BigInt->new($f2a);
+          next unless $a->copy->bmodpow($n-1, $n) == 1;
+          next unless Math::BigInt::bgcd($a->copy->bmodpow(($n-1)/2, 
$n)->bsub(1), $n) == 1;
+          @cert = ("$n", "n-1", [2, [@$certref]], [$f2a, $try_a]);
+          return ($n, \@cert);
+        }
       }
       # Didn't pass the selected a values.  Try another R.
     }
@@ -1600,46 +1624,51 @@ sub is_prob_prime {
   return ($n <= 18446744073709551615)  ?  2  :  1;
 }
 
-
+# Return just the non-cert portion.
 sub is_provable_prime {
-  my($n, $ref_proof) = @_;
+  my($n) = @_;
   return 0 if defined $n && $n < 2;
   _validate_positive_integer($n);
 
-  if (defined $ref_proof) {
-    croak "Second argument must be an array ref" if ref($ref_proof) ne 'ARRAY';
-    @$ref_proof = ();
-  }
+  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_prime, $cert) = is_provable_prime_with_cert($n);
+  return $is_prime;
+}
+
+# Return just the cert portion.
+sub prime_certificate {
+  my($n) = @_;
+  my ($is_prime, $cert) = is_provable_prime_with_cert($n);
+  return @$cert;
+}
+
+
+sub is_provable_prime_with_cert {
+  my($n) = @_;
+  return 0 if defined $n && $n < 2;
+  _validate_positive_integer($n);
 
   # Set to 0 if you want the proof to go down to 11.
   if (1) {
-    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;
-      }
-      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;
+    if (ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL) {
+      my $isp = _XS_is_prime($n);
+      return ($isp == 2) ? ($isp, [$n]) : ($isp, []);
+    }
+    if ($_HAVE_GMP && defined 
&Math::Prime::Util::GMP::is_provable_prime_with_cert) {
+      return Math::Prime::Util::GMP::is_provable_prime_with_cert($n);
     }
 
-    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;
+    my $isp = is_prob_prime($n);
+    if ($isp != 1) {
+      return ($isp == 2) ? ($isp, [$n]) : ($isp, []);
     }
   } else {
     if ($n <= 10) {
       if ($n==2||$n==3||$n==5||$n==7) {
-        @$ref_proof = ($n) if defined $ref_proof;
-        return 2;
+        return (2, [$n]);
       }
       return 0;
     }
@@ -1648,32 +1677,23 @@ sub is_provable_prime {
   # 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
+  #   BLS75 combo  Corollary 11 of BLS75.  Trial factor n-1 and n+1 to B, find
+  #                factors F1 of n-1 and F2 of n+1.  Quit when:
+  #                B > (N/(F1*F1*(F2/2)))^1/3 or B > (N/((F1/2)*F2*F2))^1/3
+  #   BLS75 n+1    Requires factoring n+1 to (n/2)^1/3 (theorem 19)
+  #   BLS75 n-1    Requires factoring n-1 to (n/2)^1/3 (theorem 5 or 7)
+  #   Pocklington  Requires factoring n-1 to n^1/2 (BLS75 theorem 4)
+  #   Lucas        Easy, requires factoring of n-1 (BLS75 theorem 1)
   #   AKS          horribly slow
   # See http://primes.utm.edu/prove/merged.html or other sources.
 
   #my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_lucas($n);
   my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_bls75($n);
-
-  @$ref_proof = @$pref if defined $ref_proof;
   carp "proved $n is not prime\n" if !$isp;
-  return $isp;
+  return ($isp, $pref);
 }
 
 
-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'};
@@ -1681,6 +1701,9 @@ sub verify_prime {
   # Empty input = no certificate = not prime
   return 0 if scalar @pdata == 0;
 
+  # Handle case of being handed a reference to the certificate.
+  @pdata = @{$pdata[0]} if scalar @pdata == 1 && ref($pdata[0]) eq 'ARRAY';
+
   my $n = shift @pdata;
   if (length($n) == 1) {
     return 1 if $n =~ /^[2357]$/;
@@ -1804,7 +1827,7 @@ sub verify_prime {
         $f = Math::BigInt->new("$farray->[0]");
         return 0 unless verify_prime(@$farray);
       } else {
-        $f = $farray;
+        $f = Math::BigInt->new("$farray");
         return 0 unless verify_prime($f);
       }
       next if defined $factors_seen{"$f"};   # repeated factors
@@ -1821,7 +1844,13 @@ sub verify_prime {
       $factors_seen{"$f"} = 1;
     }
     croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
-    croak "BLS75 error: $A not even" unless $A->is_even();
+
+    # The theorems state that A is the even portion, so we are requiring 2 be
+    # listed as a factor.
+    if ($A->is_odd) {
+      print "primality fail: 2 must be included as a factor" if $verbose;
+      return 0;
+    }
 
     # TODO: consider: if B=1 and a single a is given, then Lucas test.
 
@@ -2777,18 +2806,16 @@ exist, there is a weak conjecture (Martin) that none 
exist under 10000 digits.
 
 Takes a positive number as input and returns back either 0 (composite),
 2 (definitely prime), or 1 (probably prime).  This gives it the same return
-values as C<is_prime> and C<is_prob_prime>.
+values as L</is_prime> and L</is_prob_prime>.
 
-The current implementation uses a Lucas test requiring a complete factorization
-of C<n-1>, which may not be possible in a reasonable amount of time.  The GMP
-version uses the BLS (Brillhart-Lehmer-Selfridge) method, requiring C<n-1> to
-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.
+The current implementation of both the Perl and GMP proofs is using theorem 5
+of BLS75 (Brillhart-Lehmer-Selfridge), requiring C<n-1> to be factored to
+C<(n/2)^(1/3))>.  This takes less time than factoring to C<n^0.5> as required
+by the generalized Pocklington test or C<n-1> for the Lucas test.  However it
+is possible a factor cannot be found in a reasonable amount of time, so you
+should always test that the result in C<2> to ensure it was 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.
+A later implementation will use an ECPP test for larger inputs.
 
 
 =head2 prime_certificate
@@ -2802,6 +2829,14 @@ This may be examined or given to L</verify_prime> for 
verification.  The latter
 function contains the description of the format.
 
 
+=head2 is_provable_prime_with_cert
+
+Given a positive integer as input, returns a two element array containing the
+result of L</is_provable_prime> and an array reference containing the primality
+certificate like L</prime_certificate>.  The certificate will be an empty
+array reference if the result is not 2 (definitely prime).
+
+
 =head2 verify_prime
 
   my @cert = prime_certificate($n);
@@ -3234,6 +3269,7 @@ Similar to L</random_nbit_prime>, the result will be a 
BigInt if the
 number of bits is greater than the native bit size.  For better performance
 with large bit sizes, install L<Math::Prime::Util::GMP>.
 
+
 =head2 random_maurer_prime
 
   my $bigprime = random_maurer_prime(512);
@@ -3247,6 +3283,23 @@ with large bit sizes, install L<Math::Prime::Util::GMP>.
 The differences between this function and that in L<Crypt::Primes> are
 described in the L</"SEE ALSO"> section.
 
+Internally this additionally runs the BPSW probable prime test on every
+partial result, and constructs a primality certificate for the final
+result, which is verified.  These add additional checks that the resulting
+value has been properly constructed.
+
+
+=head2 random_maurer_prime_with_cert
+
+  my($n, $cert_ref) = random_maurer_prime_with_cert(512)
+
+As with L</random_maurer_prime>, but returns a two-element array containing
+the n-bit provable prime along with a primality certificate.  The certificate
+is the same as produced by L</prime_certificate> or
+L</is_provable_prime_with_cert>, and can be parsed by L</verify_prime> or
+any other software that can parse the certificate (the "n-1" form is described
+in detail in L</verify_prime>).
+
 
 
 =head1 UTILITY FUNCTIONS
diff --git a/lib/Math/Prime/Util/ECProjectivePoint.pm 
b/lib/Math/Prime/Util/ECProjectivePoint.pm
index 5af473f..1d062d8 100644
--- a/lib/Math/Prime/Util/ECProjectivePoint.pm
+++ b/lib/Math/Prime/Util/ECProjectivePoint.pm
@@ -187,7 +187,7 @@ sub copy {
 __END__
 
 
-# ABSTRACT: Elliptic curve operations for affine points
+# ABSTRACT: Elliptic curve operations for projective points
 
 =pod
 
@@ -196,7 +196,7 @@ __END__
 
 =head1 NAME
 
-Math::Prime::Util::ECAffinePoint - Elliptic curve operations for affine points
+Math::Prime::Util::ECProjectivePoint - Elliptic curve operations for 
projective points
 
 
 =head1 VERSION
@@ -207,7 +207,7 @@ 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);
+  my $ECP = Math::Prime::Util::ECProjectivePoint->new($a, $b, $n, 0, 1);
 
   # scalar multiplication by k.
   $ECP->mul($k)
@@ -228,7 +228,7 @@ To write.
 
 =head2 new
 
-  $point = Math::Prime::Util::EllipticCurve->new(a, b);
+  $point = Math::Prime::Util::ECProjectivePoint->new(a, b);
 
 Returns a new curve defined by a and b.
 
@@ -242,9 +242,9 @@ Returns the C<a>, C<b>, or C<n> values that describe the 
curve.
 
 =head2 x
 
-=head2 y
+=head2 z
 
-Returns the C<x> or C<y> values that define the point on the curve.
+Returns the C<x> or C<z> values that define the point on the curve.
 
 =head2 f
 
@@ -254,6 +254,10 @@ Returns a possible factor found during EC multiplication.
 
 Takes another point on the same curve as an argument and adds it this point.
 
+=head2 double
+
+Double the current point on the curve.
+
 =head2 mul
 
 Takes an integer and performs scalar multiplication of the point.
@@ -263,6 +267,15 @@ Takes an integer and performs scalar multiplication of the 
point.
 Returns true if the point is (0,1), which is the point at infinity for
 the affine coordinates.
 
+=head2 copy
+
+Returns a copy of the point.
+
+=head2 normalize
+
+Performs an extended gcd operation to make C<z=1>.  If a factor of C<n> is
+found it is put in C<f>.
+
 
 =head1 SEE ALSO
 
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index b94e684..d8a7bcd 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1815,7 +1815,7 @@ sub primality_proof_lucas {
 
   # 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 (2, [$n]) if $n < 4;
   return @composite if miller_rabin($n,2,15,325) == 0;
 
   if (!defined $Math::BigInt::VERSION) {
@@ -1843,12 +1843,12 @@ sub primality_proof_lucas {
     # 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) {
+      my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
+      if ($isp != 2) {
         carp "could not prove primality of $n.\n";
         return (1, []);
       }
-      push @fac_proofs, (scalar @fproof == 1) ? @fproof : [@fproof];
+      push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
     }
     my @proof = ("$n", "Pratt", [@fac_proofs], $a);
     return (2, [@proof]);
@@ -1856,13 +1856,15 @@ sub primality_proof_lucas {
   return @composite;
 }
 
+use Data::Dump qw/dump/;
 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 (2, [$n]) if $n < 4;
+  return @composite if ($n & 1) == 0;
   return @composite if miller_rabin($n,2,15,325) == 0;
 
   if (!defined $Math::BigInt::VERSION) {
@@ -1874,9 +1876,11 @@ sub primality_proof_bls75 {
   my $nm1 = $n-1;
   my $A = $nm1->copy->bone;   # factored part
   my $B = $nm1->copy;         # unfactored part
-  my @factors;
+  my @factors = (2);
+  croak "BLS75 error: n-1 not even" unless $nm1->is_even();
   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) {
+    last if $f*$f > $B;
     if (($B % $f) == 0) {
       push @factors, $f;
       do { $B /= $f;  $A *= $f; } while (($B % $f) == 0);
@@ -1884,7 +1888,9 @@ sub primality_proof_bls75 {
   }
   my @nstack;
   # nstack should only hold composites
-  if (Math::Prime::Util::is_prob_prime($B)) {
+  if ($B == 1) {
+    # Completely factored.  Nothing.
+  } elsif (Math::Prime::Util::is_prob_prime($B)) {
     push @factors, $B;
     $A *= $B;  $B /= $B;   # completely factored already
   } else {
@@ -1925,8 +1931,7 @@ sub primality_proof_bls75 {
     }
   }
   { # remove duplicate factors and make a sorted array of bigints
-    my %uf;
-    undef @uf{@factors};
+    my %uf = map { $_ => 1 } @factors;
     @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
   }
   # Did we factor enough?
@@ -1955,14 +1960,14 @@ sub primality_proof_bls75 {
       last;
     }
     return @composite unless $success;
-    my @fproof;
-    if (Math::Prime::Util::is_provable_prime($f, \@fproof) != 2) {
+    my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
+    if ($isp != 2) {
       carp "could not prove primality of $n.\n";
       return (1, []);
     }
-    push @fac_proofs, (scalar @fproof == 1) ? @fproof : [@fproof];
+    push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
   }
-  my @proof = ("$n", "n-1", [@fac_proofs], [@as]);
+  my @proof = ($n, "n-1", [@fac_proofs], [@as]);
   return (2, [@proof]);
 }
 
@@ -2625,6 +2630,20 @@ Takes a positive number as input, and returns 1 if the 
input can be proven
 prime using the AKS primality test.  This code is included for completeness
 and as an example, but is impractically slow.
 
+=head2 primality_proof_lucas
+
+Given a positive number C<n> as input, performs a full factorization of C<n-1>,
+then attempts a Lucas test on the result.  A Pratt-style certificate is
+returned.  Note that if the input is composite, this will take a B<very> long
+time to return.
+
+=head2 primality_proof_bls75
+
+Given a positive number C<n> as input, performs a partial factorization of
+C<n-1>, then attempts a proof using theorem 5 of Brillhart, Lehmer, and
+Selfridge's 1975 paper.  This can take a long time to return if given a
+composite, though it should not be anywhere near as long as the Lucas test.
+
 
 =head1 UTILITY FUNCTIONS
 
@@ -2730,6 +2749,19 @@ given is to run with increasing B1 factors.
 This method can rapidly find a factor C<p> of C<n> where C<p-1> is smooth
 (i.e. C<p-1> has no large factors).
 
+=head2 ecm_factor
+
+  my @factors = ecm_factor($n);
+  my @factors = ecm_factor($n, 1000);           # B1 = B2 = 1000, curves = 10
+  my @factors = ecm_factor($n, 1000, 5000, 10); # Set B1, B2, and ncurves
+
+Given a positive number input, tries to discover a factor using ECM.  The
+resulting array will contain either two factors (it succeeded) or the original
+number (no factor was found).  In either case, multiplying @factors yields the
+original input.  The B1 and B2 smoothness factors for stage 1 and stage 2
+respectively may be supplied, as can be number of curves to try.
+
+
 
 
 =head1 MATHEMATICAL FUNCTIONS

-- 
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