yann@1250: Original patch from: perfpow.c.diff
yann@1250: 
yann@1250: -= BEGIN original header =-
yann@1250: Copyright 1998, 1999, 2000, 2001, 2005, 2008 Free Software Foundation, Inc.
yann@1250: 
yann@1250: This file is part of the GNU MP Library.
yann@1250: 
yann@1250: The GNU MP Library is free software; you can redistribute it and/or modify
yann@1250: it under the terms of the GNU Lesser General Public License as published by
yann@1250: the Free Software Foundation; either version 3 of the License, or (at your
yann@1250: option) any later version.
yann@1250: 
yann@1250: The GNU MP Library is distributed in the hope that it will be useful, but
yann@1250: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
yann@1250: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
yann@1250: License for more details.
yann@1250: 
yann@1250: You should have received a copy of the GNU Lesser General Public License
yann@1250: along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
yann@1250: 
yann@1250: -= END original header =-
yann@1250: 
yann@1250: diff -durN gmp-4.2.4.orig/mpz/perfpow.c gmp-4.2.4/mpz/perfpow.c
yann@1250: --- gmp-4.2.4.orig/mpz/perfpow.c	2007-08-30 20:31:41.000000000 +0200
yann@1250: +++ gmp-4.2.4/mpz/perfpow.c	2009-03-08 18:36:16.000000000 +0100
yann@1250: @@ -1,7 +1,7 @@
yann@1250:  /* mpz_perfect_power_p(arg) -- Return non-zero if ARG is a perfect power,
yann@1250:     zero otherwise.
yann@1250:  
yann@1250: -Copyright 1998, 1999, 2000, 2001, 2005 Free Software Foundation, Inc.
yann@1250: +Copyright 1998, 1999, 2000, 2001, 2005, 2008 Free Software Foundation, Inc.
yann@1250:  
yann@1250:  This file is part of the GNU MP Library.
yann@1250:  
yann@1250: @@ -59,6 +59,8 @@
yann@1250:  #define SMALLEST_OMITTED_PRIME 1009
yann@1250:  
yann@1250:  
yann@1250: +#define POW2P(a) (((a) & ((a) - 1)) == 0)
yann@1250: +
yann@1250:  int
yann@1250:  mpz_perfect_power_p (mpz_srcptr u)
yann@1250:  {
yann@1250: @@ -72,16 +74,13 @@
yann@1250:    mp_size_t usize = SIZ (u);
yann@1250:    TMP_DECL;
yann@1250:  
yann@1250: -  if (usize == 0)
yann@1250: -    return 1;			/* consider 0 a perfect power */
yann@1250: +  if (mpz_cmpabs_ui (u, 1) <= 0)
yann@1250: +    return 1;			/* -1, 0, and +1 are perfect powers */
yann@1250:  
yann@1250:    n2 = mpz_scan1 (u, 0);
yann@1250:    if (n2 == 1)
yann@1250:      return 0;			/* 2 divides exactly once.  */
yann@1250:  
yann@1250: -  if (n2 != 0 && (n2 & 1) == 0 && usize < 0)
yann@1250: -    return 0;			/* 2 has even multiplicity with negative U */
yann@1250: -
yann@1250:    TMP_MARK;
yann@1250:  
yann@1250:    uns = ABS (usize) - n2 / BITS_PER_MP_LIMB;
yann@1250: @@ -89,6 +88,14 @@
yann@1250:    MPZ_TMP_INIT (u2, uns);
yann@1250:  
yann@1250:    mpz_tdiv_q_2exp (u2, u, n2);
yann@1250: +  mpz_abs (u2, u2);
yann@1250: +
yann@1250: +  if (mpz_cmp_ui (u2, 1) == 0)
yann@1250: +    {
yann@1250: +      TMP_FREE;
yann@1250: +      /* factoring completed; consistent power */
yann@1250: +      return ! (usize < 0 && POW2P(n2));
yann@1250: +    }
yann@1250:  
yann@1250:    if (isprime (n2))
yann@1250:      goto n2prime;
yann@1250: @@ -97,6 +104,9 @@
yann@1250:      {
yann@1250:        prime = primes[i];
yann@1250:  
yann@1250: +      if (mpz_cmp_ui (u2, prime) < 0)
yann@1250: +	break;
yann@1250: +
yann@1250:        if (mpz_divisible_ui_p (u2, prime))	/* divisible by this prime? */
yann@1250:  	{
yann@1250:  	  rem = mpz_tdiv_q_ui (q, u2, prime * prime);
yann@1250: @@ -115,12 +125,6 @@
yann@1250:  	      n++;
yann@1250:  	    }
yann@1250:  
yann@1250: -	  if ((n & 1) == 0 && usize < 0)
yann@1250: -	    {
yann@1250: -	      TMP_FREE;
yann@1250: -	      return 0;		/* even multiplicity with negative U, reject */
yann@1250: -	    }
yann@1250: -
yann@1250:  	  n2 = gcd (n2, n);
yann@1250:  	  if (n2 == 1)
yann@1250:  	    {
yann@1250: @@ -128,10 +132,11 @@
yann@1250:  	      return 0;		/* we have multiplicity 1 of some factor */
yann@1250:  	    }
yann@1250:  
yann@1250: -	  if (mpz_cmpabs_ui (u2, 1) == 0)
yann@1250: +	  if (mpz_cmp_ui (u2, 1) == 0)
yann@1250:  	    {
yann@1250:  	      TMP_FREE;
yann@1250: -	      return 1;		/* factoring completed; consistent power */
yann@1250: +	      /* factoring completed; consistent power */
yann@1250: +	      return ! (usize < 0 && POW2P(n2));
yann@1250:  	    }
yann@1250:  
yann@1250:  	  /* As soon as n2 becomes a prime number, stop factoring.
yann@1250: @@ -169,6 +174,10 @@
yann@1250:    else
yann@1250:      {
yann@1250:        unsigned long int nth;
yann@1250: +
yann@1250: +      if (usize < 0 && POW2P(n2))
yann@1250: +	return 0;
yann@1250: +
yann@1250:        /* We found some factors above.  We just need to consider values of n
yann@1250:  	 that divides n2.  */
yann@1250:        for (nth = 2; nth <= n2; nth++)
yann@1250: @@ -184,8 +193,11 @@
yann@1250:  	    exact = mpz_root (q, u2, nth);
yann@1250:  	  if (exact)
yann@1250:  	    {
yann@1250: -	      TMP_FREE;
yann@1250: -	      return 1;
yann@1250: +	      if (! (usize < 0 && POW2P(nth)))
yann@1250: +		{
yann@1250: +		  TMP_FREE;
yann@1250: +		  return 1;
yann@1250: +		}
yann@1250:  	    }
yann@1250:  	  if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0)
yann@1250:  	    {
yann@1250: @@ -199,6 +211,9 @@
yann@1250:      }
yann@1250:  
yann@1250:  n2prime:
yann@1250: +  if (usize < 0 && POW2P(n2))
yann@1250: +    return 0;
yann@1250: +
yann@1250:    exact = mpz_root (NULL, u2, n2);
yann@1250:    TMP_FREE;
yann@1250:    return exact;