[libmath-prime-util-perl] 05/55: Better AKS testing, change to znorder from naive order

[Partha P. Mukherjee]

This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.41
in repository libmath-prime-util-perl.

commit 0286fdb38539e63b34b15f17f8ea4c6bed83653e
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Apr 24 15:08:28 2014 -0700

    Better AKS testing, change to znorder from naive order
---
 Changes             |  2 ++
 aks.c               | 14 ++++++++------
 t/22-aks-prime.t    |  8 ++++++--
 xt/primality-aks.pl | 45 +++++++++++++++++++++------------------------
 4 files changed, 37 insertions(+), 32 deletions(-)

diff --git a/Changes b/Changes
index 5f8c641..47db7fa 100644
--- a/Changes
+++ b/Changes
@@ -4,6 +4,8 @@ Revision history for Perl module Math::Prime::Util
 
     - Small improvement to twin_prime_count_approx and nth_twin_prime_approx.
 
+    - Better AKS testing in xt/primality-aks.pl.
+
 
 0.40  2014-04-21
 
diff --git a/aks.c b/aks.c
index 81d9d06..e478503 100644
--- a/aks.c
+++ b/aks.c
@@ -82,8 +82,10 @@ static double log_gamma(double x)
 }
 #undef lgamma
 #define lgamma(x) log_gamma(x)
-#else
-/* Naive znorder.  Works well here because limit will be very small. */
+#endif
+
+#if 0
+/* Naive znorder.  Works well if limit is small.  Note arguments.  */
 static UV order(UV r, UV n, UV limit) {
   UV j;
   UV t = 1;
@@ -94,9 +96,6 @@ static UV order(UV r, UV n, UV limit) {
   }
   return j;
 }
-#endif
-
-#if 0
 static void poly_print(UV* poly, UV r)
 {
   int i;
@@ -270,7 +269,7 @@ int _XS_is_aks_prime(UV n)
       if (r > sqrtn)
         return 1;
 #endif
-      if (order(r, n, limit) > limit)
+      if (znorder(n, r) > limit)
         break;
     }
 
@@ -319,6 +318,9 @@ int _XS_is_aks_prime(UV n)
 
   if (verbose) { printf("# aks r = %lu  s = %lu\n", (unsigned long) r, (unsigned long) s); }
 
+  /* Almost every composite will get recognized by the first test.
+   * However, we need to run 's' tests to have the result proven for all n
+   * based on the theorems we have available at this time. */
   for (a = 1; a <= s; a++) {
     if (! test_anr(a, n, r) )
       return 0;
diff --git a/t/22-aks-prime.t b/t/22-aks-prime.t
index 8aa19bd..fdba79f 100644
--- a/t/22-aks-prime.t
+++ b/t/22-aks-prime.t
@@ -15,6 +15,12 @@ plan tests =>   6   # range
               + 1*$extra
               + 0;
 
+# Note: AKS testing is *extremely* sparse due to its lack of speed.
+#       This does almost nothing to test whether AKS is working properly.
+# 
+# If you are concerned about AKS correctness, you really need to use
+# the xt/primality-aks.pl test.
+
 ok(!eval { is_aks_prime(undef); }, "is_prime(undef)");
 ok( is_aks_prime(2),  '2 is prime');
 ok(!is_aks_prime(1),  '1 is not prime');
@@ -29,8 +35,6 @@ is( is_aks_prime(877), 1, "is_aks_prime(877) is true" );
 # UltraSPARC).  With the B+V improvements this should be fast enough for
 # the little example that we are ok.
 
-#diag "Unfortunately these tests are very slow.";
-
 SKIP: {
   # If we're pure Perl, then this is definitely too slow.
   # Arguably we should check to see if they have the GMP code.
diff --git a/xt/primality-aks.pl b/xt/primality-aks.pl
index 4107acc..371583a 100755
--- a/xt/primality-aks.pl
+++ b/xt/primality-aks.pl
@@ -1,34 +1,31 @@
 #!/usr/bin/env perl
 use strict;
 use warnings;
+use Math::Prime::Util qw/is_aks_prime is_prime primes/;
 $| = 1;  # fast pipes
 
-my $limit = shift || 10_000_000;
+my $limit = shift || 2_000_000_000;
+my $nrand = 8000;
 
-use Math::Prime::Util qw/is_aks_prime next_prime/;
+my %isprime = map { $_ => 1 } @{primes(160_000)};
 
-# It doesn't really matter for this use, but try to use independent code for
-# next_prime.  Math::Prime::FastSieve works well.
-my $sieve;
-if (eval { require Math::Prime::FastSieve; Math::Prime::FastSieve->import(); require Math::Prime::FastSieve::Sieve; 1; }) {
-  $sieve = Math::Prime::FastSieve::Sieve->new($limit + 10_000);
+print "Testing AKS for all numbers from 1 to 160,000:\n";
+foreach my $n (1 .. 160_000) {
+  print "." unless $n % 2000;
+  if ($isprime{$n}) {
+    die "\n$n is prime\n" unless is_aks_prime($n);
+  } else {
+    die "\n$n is composite\n" if is_aks_prime($n);
+  }
 }
-
-if (1) {
-  my $n = 2;
-  my $i = 1;
-  while ($n <= $limit) {
-    print "$n\n" unless $i++ % 1000;
-    die "$n should be prime" unless is_aks_prime($n);
-    my $next = (defined $sieve) ? $sieve->nearest_ge( $n+1 ) : next_prime($n);
-    my $diff = ($next - $n) >> 1;
-    if ($diff > 1) {
-      foreach my $d (1 .. $diff-1) {
-        my $cn = $n + 2*$d;
-        die "$cn should be composite" if is_aks_prime($cn);
-      }
-    }
-    $n = $next;
+print "\n";
+print "Testing $nrand random numbers from 1 to $limit:\n";
+for (1 .. $nrand) {
+  print "." unless $_ % 100;
+  my $n = 1 + int(rand($limit));
+  if (is_prime($n)) {
+    die "\n$n is prime\n" unless is_aks_prime($n);
+  } else {
+    die "\n$n is composite\n" if is_aks_prime($n);
   }
-  print "Success to $limit!\n";
 }

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git