[libmath-prime-util-perl] 03/04: New perfect power implementation

[Partha P. Mukherjee]

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

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

commit e2794ed6220b89a425ba1a0e141baa93913104f1
Author: Dana Jacobsen <dana at acm.org>
Date:   Sat Feb 2 23:54:33 2013 -0800

    New perfect power implementation
---
 Changes |  2 ++
 TODO    |  3 ---
 aks.c   | 52 +++++++++++++++++++++++++++++++++++++++++++---------
 3 files changed, 45 insertions(+), 12 deletions(-)

diff --git a/Changes b/Changes
index 11d6ef6..839f080 100644
--- a/Changes
+++ b/Changes
@@ -7,6 +7,8 @@ Revision history for Perl extension Math::Prime::Util.
     - Replaced fast sqrt detection in PP.pm with a slightly slower version.
       The bloom filter doesn't work right in 32-bit Perl.
 
+    - Fix is_perfect_power in XS AKS.
+
 0.19  1 February 2012
 
     - Update MR bases with newest from http://miller-rabin.appspot.com/.
diff --git a/TODO b/TODO
index 5477e1b..0957ee4 100644
--- a/TODO
+++ b/TODO
@@ -45,7 +45,4 @@
   each), then use it to get 32-bit irands.  This would only be used if they
   didn't give us a RNG (so they don't care about strict crypto).
 
-- Perfect power from Dietzfelbinger algorithm 2.3.5 works a bit better.
-  Newton's method would be even better.
-
 - Add PE 35 Circular primes and PE 291 Panatoipal primes to bin/primes.pl
diff --git a/aks.c b/aks.c
index b9be560..3097ecb 100644
--- a/aks.c
+++ b/aks.c
@@ -2,6 +2,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
+#include <float.h>
 
 /*
  * The AKS v6 algorithm, for native integers.  Based on the AKS v6 paper.
@@ -37,16 +38,49 @@ static UV log2floor(UV n) {
   return log2n;
 }
 
-/* See Bach and Sorenson (1993) for much better algorithm */
-/* This isn't generally going to work for numbers > 2**53 */
-static int is_perfect_power(UV x) {
+/* Bach and Sorenson (1993) would be better */
+static int is_perfect_power(UV n) {
   UV b, last;
-  if ((x & (x-1)) == 0)  return 1;          /* powers of 2    */
-  b = sqrt(x)+0.5; if (b*b == x)  return 1; /* perfect square */
-  last = log2floor(x) + 1;
-  for (b = 3; b < last; b = _XS_next_prime(b)) {
-    UV root = pow(x, 1.0 / (double)b) + 0.5;
-    if ( ((UV)(pow(root, b)+0.5)) == x)  return 1;
+  if ((n <= 3) || (n == UV_MAX)) return 0;
+  if ((n & (n-1)) == 0)  return 1;          /* powers of 2    */
+  last = log2floor(n-1) + 1;
+#if (BITS_PER_WORD == 32) || (DBL_DIG > 19)
+  if (1) {
+#elif DBL_DIG == 10
+  if (n < UVCONST(10000000000)) {
+#elif DBL_DIG == 15
+  if (n < UVCONST(1000000000000000)) {
+#else
+  if ( n < (UV) pow(10, DBL_DIG) ) {
+#endif
+    /* Simple floating point method.  Fast, but need enough mantissa. */
+    b = sqrt(n)+0.5; if (b*b == n)  return 1; /* perfect square */
+    for (b = 3; b < last; b = _XS_next_prime(b)) {
+      UV root = pow(n, 1.0 / (double)b) + 0.5;
+      if ( ((UV)(pow(root, b)+0.5)) == n)  return 1;
+    }
+  } else {
+    /* Dietzfelbinger, algorithm 2.3.5 (without optimized exponential) */
+    for (b = 2; b <= last; b++) {
+      UV a = 1;
+      UV c = n;
+      while (c >= HALF_WORD) c = (1+c)>>1;
+      while ((c-a) >= 2) {
+        UV m, maxm, p, i;
+        m = (a+c) >> 1;
+        maxm = UV_MAX / m;
+        p = m;
+        for (i = 2; i <= b; i++) {
+          if (p > maxm) { p = n+1; break; }
+          p *= m;
+        }
+        if (p == n)  return 1;
+        if (p < n)
+          a = m;
+        else
+          c = m;
+      }
+    }
   }
   return 0;
 }

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