[libmath-prime-util-perl] 43/54: Custom lgamma so improved AKS works on all platforms

[Partha P. Mukherjee]

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

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

commit 3783549f9d032d48354d74620c8ca4b2a8c2ebf8
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Feb 27 13:34:36 2014 -0800

    Custom lgamma so improved AKS works on all platforms
---
 Changes             |  4 ++--
 aks.c               | 34 +++++++++++++++++++++++++---------
 t/22-aks-prime.t    |  6 +-----
 xt/primality-aks.pl | 16 +++++++++++-----
 4 files changed, 39 insertions(+), 21 deletions(-)

diff --git a/Changes b/Changes
index 26238c2..78b1691 100644
--- a/Changes
+++ b/Changes
@@ -10,8 +10,8 @@ Revision history for Perl module Math::Prime::Util
 
     - Factoring powers is much faster.
 
-    - Add Bernstein+Voloch improvements to AKS.
-      Much faster, though of course still much slower than BPSW.
+    - Add Bernstein+Voloch improvements to AKS.  Much faster than the v6
+      implementation, though still terribly slow vs. BPSW or other proofs.
 
     [OTHER]
 
diff --git a/aks.c b/aks.c
index ed754b5..4154a8b 100644
--- a/aks.c
+++ b/aks.c
@@ -29,13 +29,8 @@
 
 #define SQRTN_SHORTCUT 1
 
-/* Use improvements from Bornemann's 2002 implementation if we have lgamma */
-#if !defined(_MSC_VER) && \
-    (defined(__USE_ISOC99) || defined(_ISOC99_SOURCE) || defined(_NETBSD_SOURCE) || defined(__cplusplus) || !defined(__STRICT_ANSI__) || __STDC_VERSION__ >= 199901L || _POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600)
+/* Use improvements from Bornemann's 2002 implementation */
 #define IMPL_BORNEMANN 1
-#else
-#define IMPL_BORNEMANN 0
-#endif
 
 #include "ptypes.h"
 #include "aks.h"
@@ -58,6 +53,28 @@ static int is_primitive_root(UV n, UV r)
   }
   return 1;
 }
+/* We could use lgamma, but it isn't in MSVC and not in pre-C99.  The only
+ * sure way to find if it is available is test compilation (ala autoconf).
+ * Instead, we'll just use our own implementation.
+ * See http://mrob.com/pub/ries/lanczos-gamma.html for alternates. */
+static double lanczos_coef[8+1] =
+{ 0.99999999999980993, 676.5203681218851, -1259.1392167224028,
+  771.32342877765313, -176.61502916214059, 12.507343278686905,
+  -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 };
+static double log_sqrt_two_pi =  0.91893853320467274178;
+static double log_gamma(double x)
+{
+  double base = x + 7 + 0.5;
+  double sum = 0;
+  int i;
+  for (i = 8; i >= 1; i--)
+    sum += lanczos_coef[i] / (x + (double)i);
+  sum += lanczos_coef[0];
+  sum = log_sqrt_two_pi + logl(sum/x) + ( (x+0.5)*logl(base) - base );
+  return sum;
+}
+#undef lgamma
+#define lgamma(x) log_gamma(x)
 #else
 /* Naive znorder.  Works well here because limit will be very small. */
 static UV order(UV r, UV n, UV limit) {
@@ -254,8 +271,6 @@ int _XS_is_aks_prime(UV n)
       return 1;
 
     s = (UV) floor(sqrt(r-1) * log2n);
-
-    if (verbose) { printf("# aks r = %lu  s = %lu\n", (unsigned long) r, (unsigned long) s); }
   }
 #else
   {
@@ -293,9 +308,10 @@ int _XS_is_aks_prime(UV n)
     if (slim >= HALF_WORD || (slim*slim) >= n)
       return 1;
   }
-  if (verbose) { printf("# aks r = %lu  s = %lu\n", (unsigned long) r, (unsigned long) s); }
 #endif
 
+  if (verbose) { printf("# aks r = %lu  s = %lu\n", (unsigned long) r, (unsigned long) s); }
+
   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 58089aa..7418599 100644
--- a/t/22-aks-prime.t
+++ b/t/22-aks-prime.t
@@ -37,11 +37,7 @@ SKIP: {
   # If we're pure Perl, then this is definitely too slow.
   # Arguably we should check to see if they have the GMP code.
   skip "Skipping PP AKS on PP -- just too slow.", 1 if $ispp;
-  # If we have 64-bit available in the compiler (e.g. uint64_t), this can
-  # still be quite fast.  However for pretty much everyone else, this is
-  # just far too slow for running in a test suite.
-  skip "Skipping PP AKS on 32-bit -- just too slow.", 1 if !$use64;
-  # The first number that makes it past the sqrt test to actually run.
+  # The least number that performs the full test with either implementation.
   is( is_aks_prime(69197), 1, "is_aks_prime(69197) is true" );
 }
 
diff --git a/xt/primality-aks.pl b/xt/primality-aks.pl
index c0cf7f5..4107acc 100755
--- a/xt/primality-aks.pl
+++ b/xt/primality-aks.pl
@@ -5,16 +5,22 @@ $| = 1;  # fast pipes
 
 my $limit = shift || 10_000_000;
 
-use Math::Prime::Util qw/is_aks_prime/;
-use Math::Prime::FastSieve;
-my $sieve = Math::Prime::FastSieve::Sieve->new($limit + 10_000);
+use Math::Prime::Util qw/is_aks_prime next_prime/;
+
+# 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);
+}
 
 if (1) {
   my $n = 2;
+  my $i = 1;
   while ($n <= $limit) {
-    print "$n\n" if $n > 69000; # unless $i++ % 262144;
+    print "$n\n" unless $i++ % 1000;
     die "$n should be prime" unless is_aks_prime($n);
-    my $next = $sieve->nearest_ge( $n+1 );
+    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) {

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