Permalink
Browse files

Merge remote-tracking branch 'danaj/master'

  • Loading branch information...
2 parents 3bba000 + dd9fb91 commit 31cc8fdf1367a905bea34a2d4da423b49ba3830d @leto committed Feb 19, 2013
Showing with 195 additions and 165 deletions.
  1. +4 −0 Changes
  2. +1 −1 lib/Math/Primality.pm
  3. +58 −54 lib/Math/Primality/AKS.pm
  4. +91 −107 lib/Math/Primality/BigPolynomial.pm
  5. +33 −0 t/big_polynomial.t
  6. +8 −3 t/is_aks_prime.t
View
@@ -7,6 +7,10 @@
The AKS implementation was only doing trial division. Fix to do the
actual AKS test. Sadly this makes it incredibly slow.
+ Change AKS to the latest standard implementation (v6 of the paper).
+
+ Rewrite big parts of BigPolynomial.
+
=head2 0.08
Use a small GCD in the prime_count loop, making it much faster.
View
@@ -510,7 +510,7 @@ sub is_prime($) {
return 0 if Rmpz_cmp_ui($n, 2) == -1;
# is it a small prime? (<=257)
- return _is_small_prime($n) if Rmpz_cmp_ui($n, 257) == -1;
+ return _is_small_prime($n) if Rmpz_cmp_ui($n, 257) <= 0;
# Divisible by 2?
return 0 if Rmpz_even_p($n);
View
@@ -56,92 +56,96 @@ sub is_aks_prime($) {
# http://www.cs.cmu.edu/afs/cs/user/mjs/ftp/thesis-program/2005/rotella.pdf
# http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf
- # This code follows the Rotella 2005 implementation. Unfortunately that
- # paper was based on the original AKS algorithm, which is *very* slow. His
- # graphs show hundreds of seconds to test the primality of 3-digit numbers,
- # with a C+GMP implementation. We therefore expect this Perl+GMPz version
- # to also be completely unusable due to performance, and should be surprised
- # if it proves primality of a 5-digit prime in under an hour.
+ # The previous version of this code followed the Rotella 2005 implementation,
+ # which was based on the original AKS algorithm, hence *very* slow. His
+ # graphs show hundreds of seconds to test the primality of 3-digit numbers
+ # with a C+GMP implementation.
#
- # The newer algorithm described in the primality_v6.pdf paper is much faster,
- # though still totally impractical for any actual work without a *lot* of
- # optimization work (see Crandall and Papadopoulos for example). The
- # variant algorithms of Bernstein become practical in terms of speed, though
- # still do not match the speed of APR-CL or ECPP.
+ # The current code is done using the v6 paper so is much faster, though
+ # still impractically slow compared to other methods. See the paper of
+ # Crandall and Papadopoulos as well as Bernstein's papers for more ideas.
+ # In particular, the current bottleneck is polynomial multiplication.
my $n = GMP->new($_[0]);
+ # Step 0 - check that n is a positive integer >= 2;
+ if (Rmpz_cmp_ui($n, 2) < 0) {
+ debug "fails step 0 of aks - $n is in range\n";
+ return 0; # negative, 0, or 1.
+ }
# Step 1 - check if $n = m^d for some m, d
# if $n is a power then return 0
if (Rmpz_perfect_power_p($n)) {
- debug "fails step 1 of aks - $n is a perfect power";
+ debug "fails step 1 of aks - $n is a perfect power\n";
return 0; # composite
}
- my $r = Math::GMPz->new(2);
+
+ # This will overestimate r a little, but it's the only good way to do it
+ # in GMP without MPFR. What we want:
+ # floor( log2(n) * log2(n) )
+ # What this does:
+ # ceil(log2(n)) * ceil(log2(n))
+ # Ex: For 85991 r = 293 instead of 283, with polylimit 290 instead of 275.
my $logn = Rmpz_sizeinbase($n, 2);
my $limit = Math::GMPz->new($logn * $logn);
- Rmpz_mul_ui($limit, $limit, 4);
- # Witness search
+ # Search for first r where order(r, n) > limit
+ my $r = Math::GMPz->new(2);
while (Rmpz_cmp($r, $n) == -1) {
if(Rmpz_divisible_p($n, $r)) {
debug "$n is divisible by $r\n";
return 0;
}
- # We're following Rotella exactly. If we weren't, then doing the
- # order test followed by $r = next_prime($r) would make more sense.
- if(Rmpz_probab_prime_p($r, 5)) {
- my $i = Math::GMPz->new(1);
- my $res = Math::GMPz->new(0);
-
- my $failed = 0;
- for ( ; Rmpz_cmp($i, $limit) <= 0; Rmpz_add_ui($i, $i, 1)) {
- Rmpz_powm($res, $n, $i, $r);
- if (Rmpz_cmp_ui($res, 1) == 0) {
- $failed = 1;
- last;
- }
+ my $i = Math::GMPz->new(1);
+ my $res = Math::GMPz->new(0);
+ my $failed = 0;
+ for ( ; Rmpz_cmp($i, $limit) <= 0; Rmpz_add_ui($i, $i, 1)) {
+ Rmpz_powm($res, $n, $i, $r);
+ if (Rmpz_cmp_ui($res, 1) == 0) {
+ $failed = 1;
+ last;
}
- last if !$failed;
}
- Rmpz_add_ui($r, $r, 1);
+ last if !$failed;
+ Rmpz_nextprime($r, $r);
}
- if (Rmpz_cmp($r, $n) == 0) {
+ # We've performed trial division for every prime <= r.
+ # Hence if n > r*r then n is prime.
+ if (Rmpz_cmp($r*$r, $n) >= 0) {
debug "Found $n is prime while checking for r\n";
return 1;
}
- # Polynomial check
- my $a;
- my $sqrtr = Math::GMPz->new($r);
-
- Rmpz_sqrt($sqrtr, $r);
- my $polylimit = Math::GMPz->new(0);
- Rmpz_add_ui($polylimit, $sqrtr, 1);
- Rmpz_mul_ui($polylimit, $polylimit, $logn);
- Rmpz_mul_ui($polylimit, $polylimit, 2);
+ # Polynomial check. Limit is floor(log2n * sqrt(phi(r))).
+ # r is prime hence phi(r) = r-1.
+ my $sqrt_phi_r = Math::GMPz->new($r - 1);
+ my $plus_one = Rmpz_perfect_square_p($sqrt_phi_r) ? 0 : 1;
+ Rmpz_sqrt($sqrt_phi_r, $sqrt_phi_r);
+ my $polylimit = $logn * $sqrt_phi_r + $plus_one;
my $intr = Rmpz_get_ui($r);
+ debug "Running poly check for $n with r = $intr polylimit = $polylimit\n";
- for($a = 1; Rmpz_cmp_ui($polylimit, $a) >= 0; $a++) {
- debug "Checking at $a\n";
- my $final_size = Math::GMPz->new(0);
- Rmpz_mod($final_size, $n, $r);
- my $compare = Math::Primality::BigPolynomial->new(Rmpz_get_ui($final_size));
- $compare->setCoef(Math::GMPz->new(1), $final_size);
- $compare->setCoef(Math::GMPz->new($a), 0);
- my $res = Math::Primality::BigPolynomial->new($intr);
- my $base = Math::Primality::BigPolynomial->new(1);
- $base->setCoef(Math::GMPz->new(0), $a);
- $base->setCoef(Math::GMPz->new(1), 1);
+ my $final_size = Math::GMPz->new(0);
+ Rmpz_mod($final_size, $n, $r);
+ my $compare = Math::Primality::BigPolynomial->new();
+ $compare->setCoef(Math::GMPz->new(1), $final_size);
- Math::Primality::BigPolynomial::mpz_poly_mod_power($res, $base, $n, $n, $intr);
+ for(my $a = 1; Rmpz_cmp_ui($polylimit, $a) >= 0; $a++) {
+ debug " checking at $a / $polylimit\n";
+ my $poly = Math::Primality::BigPolynomial->new( [$a, 1 ]);
+ #print "Start: ", $poly->string, "\n";
+ $poly->powmod($n, $n, $intr);
+ #print "Finish: ", $poly->string, "\n";
+ $compare->setCoef($a, 0);
+ #print "Compare: ", $compare->string, "\n";
- if($res->isEqual($compare)) {
- debug "Found not prime at $a\n";
+ if ( ! $poly->isEqual($compare) ) {
+ debug "Found $n not prime at $a\n";
return 0;
}
}
+ debug "Found $n prime after checking $polylimit polynomials\n";
return 1;
}
Oops, something went wrong.

0 comments on commit 31cc8fd

Please sign in to comment.