Commit Martin Boij's new mpn_perfect_power_p code.

1 files changed, 3 insertions, 257 deletions

diff --git a/mpz/perfpow.c b/mpz/perfpow.c
index e66734013..81c8884aa 100644
--- a/mpz/perfpow.c
+++ b/mpz/perfpow.c
@@ -1,7 +1,8 @@
 /* mpz_perfect_power_p(arg) -- Return non-zero if ARG is a perfect power,
    zero otherwise.
 
-Copyright 1998, 1999, 2000, 2001, 2005, 2008 Free Software Foundation, Inc.
+Copyright 1998, 1999, 2000, 2001, 2005, 2008, 2009 Free Software Foundation,
+Inc.
 
 This file is part of the GNU MP Library.
 
@@ -18,266 +19,11 @@ License for more details.
 You should have received a copy of the GNU Lesser General Public License
 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
 
-/*
-  We are to determine if c is a perfect power, c = a ^ b.
-  Assume c is divisible by 2^n and that codd = c/2^n is odd.
-  Assume a is divisible by 2^m and that aodd = a/2^m is odd.
-  It is always true that m divides n.
-
-  * If n is prime, either 1) a is 2*aodd and b = n
-		       or 2) a = c and b = 1.
-    So for n prime, we readily have a solution.
-  * If n is factorable into the non-trivial factors p1,p2,...
-    Since m divides n, m has a subset of n's factors and b = n / m.
-*/
-
-/* This is a naive approach to recognizing perfect powers.
-   Many things can be improved.  In particular, we should use p-adic
-   arithmetic for computing possible roots.  */
-
-#include <stdio.h> /* for NULL */
 #include "gmp.h"
 #include "gmp-impl.h"
-#include "longlong.h"
-
-static unsigned long int gcd __GMP_PROTO ((unsigned long int, unsigned long int));
-static int isprime __GMP_PROTO ((unsigned long int));
-
-static const unsigned short primes[] =
-{  2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
-  59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,
- 137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,
- 227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,
- 313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,
- 419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,
- 509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,
- 617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,
- 727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,
- 829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
- 947,953,967,971,977,983,991,997,0
-};
-#define SMALLEST_OMITTED_PRIME 1009
-
-
-#define POW2P(a) (((a) & ((a) - 1)) == 0)
 
 int
 mpz_perfect_power_p (mpz_srcptr u)
 {
-  unsigned long int prime;
-  unsigned long int n, n2;
-  int i;
-  unsigned long int rem;
-  mpz_t u2, q;
-  int exact;
-  mp_size_t uns;
-  mp_size_t usize = SIZ (u);
-  TMP_DECL;
-
-  if (mpz_cmpabs_ui (u, 1) <= 0)
-    return 1;			/* -1, 0, and +1 are perfect powers */
-
-  n2 = mpz_scan1 (u, 0);
-  if (n2 == 1)
-    return 0;			/* 2 divides exactly once.  */
-
-  TMP_MARK;
-
-  uns = ABS (usize) - n2 / BITS_PER_MP_LIMB;
-  MPZ_TMP_INIT (q, uns);
-  MPZ_TMP_INIT (u2, uns);
-
-  mpz_tdiv_q_2exp (u2, u, n2);
-  mpz_abs (u2, u2);
-
-  if (mpz_cmp_ui (u2, 1) == 0)
-    {
-      TMP_FREE;
-      /* factoring completed; consistent power */
-      return ! (usize < 0 && POW2P(n2));
-    }
-
-  if (isprime (n2))
-    goto n2prime;
-
-  for (i = 1; primes[i] != 0; i++)
-    {
-      prime = primes[i];
-
-      if (mpz_cmp_ui (u2, prime) < 0)
-	break;
-
-      if (mpz_divisible_ui_p (u2, prime))	/* divisible by this prime? */
-	{
-	  rem = mpz_tdiv_q_ui (q, u2, prime * prime);
-	  if (rem != 0)
-	    {
-	      TMP_FREE;
-	      return 0;		/* prime divides exactly once, reject */
-	    }
-	  mpz_swap (q, u2);
-	  for (n = 2;;)
-	    {
-	      rem = mpz_tdiv_q_ui (q, u2, prime);
-	      if (rem != 0)
-		break;
-	      mpz_swap (q, u2);
-	      n++;
-	    }
-
-	  n2 = gcd (n2, n);
-	  if (n2 == 1)
-	    {
-	      TMP_FREE;
-	      return 0;		/* we have multiplicity 1 of some factor */
-	    }
-
-	  if (mpz_cmp_ui (u2, 1) == 0)
-	    {
-	      TMP_FREE;
-	      /* factoring completed; consistent power */
-	      return ! (usize < 0 && POW2P(n2));
-	    }
-
-	  /* As soon as n2 becomes a prime number, stop factoring.
-	     Either we have u=x^n2 or u is not a perfect power.  */
-	  if (isprime (n2))
-	    goto n2prime;
-	}
-    }
-
-  if (n2 == 0)
-    {
-      /* We found no factors above; have to check all values of n.  */
-      unsigned long int nth;
-      for (nth = usize < 0 ? 3 : 2;; nth++)
-	{
-	  if (! isprime (nth))
-	    continue;
-#if 0
-	  exact = mpz_padic_root (q, u2, nth, PTH);
-	  if (exact)
-#endif
-	    exact = mpz_root (q, u2, nth);
-	  if (exact)
-	    {
-	      TMP_FREE;
-	      return 1;
-	    }
-	  if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0)
-	    {
-	      TMP_FREE;
-	      return 0;
-	    }
-	}
-    }
-  else
-    {
-      unsigned long int nth;
-
-      if (usize < 0 && POW2P(n2))
-	{
-	  TMP_FREE;
-	  return 0;
-	}
-
-      /* We found some factors above.  We just need to consider values of n
-	 that divides n2.  */
-      for (nth = 2; nth <= n2; nth++)
-	{
-	  if (! isprime (nth))
-	    continue;
-	  if (n2 % nth != 0)
-	    continue;
-#if 0
-	  exact = mpz_padic_root (q, u2, nth, PTH);
-	  if (exact)
-#endif
-	    exact = mpz_root (q, u2, nth);
-	  if (exact)
-	    {
-	      if (! (usize < 0 && POW2P(nth)))
-		{
-		  TMP_FREE;
-		  return 1;
-		}
-	    }
-	  if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0)
-	    {
-	      TMP_FREE;
-	      return 0;
-	    }
-	}
-
-      TMP_FREE;
-      return 0;
-    }
-
-n2prime:
-  if (usize < 0 && POW2P(n2))
-    {
-      TMP_FREE;
-      return 0;
-    }
-
-  exact = mpz_root (NULL, u2, n2);
-  TMP_FREE;
-  return exact;
-}
-
-static unsigned long int
-gcd (unsigned long int a, unsigned long int b)
-{
-  int an2, bn2, n2;
-
-  if (a == 0)
-    return b;
-  if (b == 0)
-    return a;
-
-  count_trailing_zeros (an2, a);
-  a >>= an2;
-
-  count_trailing_zeros (bn2, b);
-  b >>= bn2;
-
-  n2 = MIN (an2, bn2);
-
-  while (a != b)
-    {
-      if (a > b)
-	{
-	  a -= b;
-	  do
-	    a >>= 1;
-	  while ((a & 1) == 0);
-	}
-      else /*  b > a.  */
-	{
-	  b -= a;
-	  do
-	    b >>= 1;
-	  while ((b & 1) == 0);
-	}
-    }
-
-  return a << n2;
-}
-
-static int
-isprime (unsigned long int t)
-{
-  unsigned long int q, r, d;
-
-  if (t < 3 || (t & 1) == 0)
-    return t == 2;
-
-  for (d = 3, r = 1; r != 0; d += 2)
-    {
-      q = t / d;
-      r = t - q * d;
-      if (q < d)
-	return 1;
-    }
-  return 0;
+  return mpn_perfect_power_p (PTR (u), SIZ (u));
 }