Add n_is_perfect_power.

fredrik-johansson · Jan 8, 2017 · 1e85c09 · 1e85c09
1 parent 9dfe24d
commit 1e85c09
Show file tree

Hide file tree

Showing 4 changed files with 386 additions and 0 deletions.
diff --git a/ulong_extras.h b/ulong_extras.h
@@ -332,6 +332,8 @@ FLINT_DLL ulong n_cbrtrem(ulong* remainder, ulong n);
 
 FLINT_DLL int n_is_perfect_power235(ulong n);
 
+FLINT_DLL int n_is_perfect_power(ulong * root, ulong n);
+
 FLINT_DLL int n_is_oddprime_small(ulong n);
 
 FLINT_DLL int n_is_oddprime_binary(ulong n);

diff --git a/ulong_extras/doc/ulong_extras.txt b/ulong_extras/doc/ulong_extras.txt
@@ -954,6 +954,12 @@ int n_is_perfect_power235(ulong n)
     root can be taken, if indicated, to determine whether the power
     of that root is exactly equal to $n$.
 
+int n_is_perfect_power(ulong * root, ulong n)
+
+    If $n = r^k$, return $k$ and set \code{root} to $r$. Note that $0$ and
+    $1$ are considered squares. No guarantees are made about $r$ or $k$
+    being the minimum possible value.
+
 ulong n_rootrem(ulong* remainder, ulong n, ulong root)
 
     This function uses the Newton iteration method to calculate the nth root of

diff --git a/ulong_extras/is_perfect_power.c b/ulong_extras/is_perfect_power.c
@@ -0,0 +1,243 @@
+/*
+    Copyright (C) 2009 Thomas Boothby
+    Copyright (C) 2009, 2016 William Hart
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include <gmp.h>
+#define ulong ulongxx /* interferes with system includes */
+#include <math.h>
+#undef ulong
+#include "flint.h"
+#include "ulong_extras.h"
+
+int n_is_perfect_power(ulong * root, ulong n)
+{
+    static unsigned char mod63[63] = {7,7,4,0,5,4,0,5,6,5,4,4,0,4,4,0,5,4,5,4,
+                              4,0,5,4,0,5,4,6,7,4,0,4,4,0,4,6,7,5,4,0,4,4,0,5,
+                              4,4,5,4,0,5,4,0,4,4,4,6,4,0,5,4,0,4,6};
+    static unsigned char mod61[61] = {7,7,0,3,1,1,0,0,2,3,0,6,1,5,5,1,1,0,0,1,
+                              3,4,1,2,2,1,0,3,2,4,0,0,4,2,3,0,1,2,2,1,4,3,1,0,
+                              0,1,1,5,5,1,6,0,3,2,0,0,1,1,3,0,7};
+    static unsigned char mod44[44] = {7,7,0,2,3,3,0,2,2,3,0,6,7,2,0,2,3,2,0,2,
+                              3,6,0,6,2,3,0,2,2,2,0,2,6,7,0,2,3,3,0,2,2,2,0,6};
+    static unsigned char mod31[31] = {7,7,3,0,3,5,4,1,3,1,1,0,0,0,1,2,3,0,1,1,
+                              1,0,0,2,0,5,4,2,1,2,6};
+    static unsigned char mod72[72] = {7,7,0,0,0,7,0,7,7,7,0,7,0,7,0,0,7,7,0,7,
+                              0,0,0,7,0,7,0,7,0,7,0,7,7,0,0,7,0,7,0,0,7,7,0,7,
+                              0,7,0,7,0,7,0,0,0,7,0,7,7,0,0,7,0,7,0,7,7,7,0,7,
+                              0,0,0,7};
+    static unsigned char mod49[49] = {1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,
+                              0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
+                              0,0,0,0,1};
+    static unsigned char mod67[67] = {2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+                              0,0,0,0,0,0,0,0,0,2,2,0,0,0,0,0,0,2,2,0,0,0,0,0,
+                              0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2};
+    static unsigned char mod79[79] = {4,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+                              0,0,0,4,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+                              0,0,0,0,0,0,0,0,0,0,0,4,4,0,0,0,0,0,0,0,0,0,0,0,
+                              0,0,0,0,0,0,0,0,0,0,4};
+
+    unsigned char t;
+    ulong count, exp, r;
+
+    /* check for powers 2, 3, 5 */
+    t = mod31[n%31];
+    t &= mod44[n%44];
+    t &= mod61[n%61];
+    t &= mod63[n%63];
+
+    if (t & 1) 
+    {
+        ulong y = n_sqrtrem(&r, n);
+        if (r == 0)
+        {
+           *root = y;
+           return 2;
+        }
+    }
+
+    if (t & 2) 
+    {
+        ulong y = n_cbrtrem(&r, n);
+        if (r == 0)
+        if (n == n_pow(y, 3))
+        {
+           *root = y;
+           return 3;
+        }
+    }
+
+    if (t & 4) 
+    {
+        ulong y = n_rootrem(&r, n, 5);
+        if (r == 0)
+        {
+           *root = y;
+           return 5;
+        }
+    }
+
+    /* check for power 7, 11, 13 */
+    t = mod49[n%49];
+    t |= mod67[n%67];
+    t |= mod79[n%79];
+    t &= mod72[n%72];
+
+    if (t & 1) 
+    {
+        ulong y = n_rootrem(&r, n, 7);
+        if (r == 0)
+        {
+           *root = y;
+           return 7;
+        }
+    }
+
+    if (t & 2) 
+    {
+        ulong y = n_rootrem(&r, n, 11);
+        if (r == 0)
+        {
+           *root = y;
+           return 11;
+        }
+    }
+
+    if (t & 13) 
+    {
+        ulong y = n_rootrem(&r, n, 13);
+        if (r == 0)
+        {
+           *root = y;
+           return 13;
+        }
+    }
+
+    /* highest power of 2 */
+    count_trailing_zeros(count, n);
+    n >>= count;
+
+    if (n == 1)
+    {
+       if (count == 1)
+          return 0;
+       *root = 2;
+       return count;
+    }
+
+    /* check other powers (exp >= 17, root <= 13 and odd) */
+    exp = 0;
+    while ((n % 3) == 0)
+    {
+       n /= 3;
+       exp += 1;
+    }
+    if (exp > 0)
+    {
+       if (n == 1 && exp > 1)
+       {
+          if (count == 0)
+          {
+             *root = 3;
+             return exp;
+          } else if (count == exp)
+          {
+             *root = 6;
+             return exp;
+          } else if (count == 2*exp)
+          {
+             *root = 12;
+             return exp;
+          }
+       }
+       return 0;
+    }
+
+#if FLINT64
+
+    exp = 0;
+    while ((n % 5) == 0)
+    {
+       n /= 5;
+       exp += 1;
+    }
+    if (exp > 0)
+    {
+       if (n == 1 && exp > 1)
+       {
+          if (count == 0)
+          {
+             *root = 5;
+             return exp;
+          } else if (count == exp)
+          {
+             *root = 10;
+             return exp;
+          }
+       }
+       return 0;
+    }
+
+    if (count > 0)
+       return 0;
+
+    exp = 0;
+    while ((n % 7) == 0)
+    {
+       n /= 7;
+       exp += 1;
+    }
+    if (exp > 0)
+    {
+       if (n == 1 && exp > 1)
+       {
+          *root = 7;
+          return exp;
+       } 
+       return 0;
+    }
+
+    exp = 0;
+    while ((n % 11) == 0)
+    {
+       n /= 11;
+       exp += 1;
+    }
+    if (exp > 0)
+    {
+       if (n == 1 && exp > 1)
+       {
+          *root = 11;
+          return exp;
+       } 
+       return 0;
+    }
+
+    exp = 0;
+    while ((n % 13) == 0)
+    {
+       n /= 13;
+       exp += 1;
+    }
+    if (exp > 0)
+    {
+       if (n == 1 && exp > 1)
+       {
+          *root = 13;
+          return exp;
+       } 
+       return 0;
+    }
+
+#endif
+
+    return 0;
+}
+