Integer root extracting

[Hans Aberg]

On 3 Sep 2010, at 12:27, Torbjorn Granlund wrote:

>  GMP has a function mpz_perfect_power_p() that for given integer a
>  tells if a = x^y is solvable in integers x, y with y > 1. But I need
>  to find the largest y for which this is true (and also finding x).  
> One
>  way to do this is to prime factor a and take the gcd of the  
> exponents,
>  but is there a faster method?
>
> Completely factoring a is of course not always possible.

I'm not sure what you mean here: a is an integer, in fact positive in  
my application. So it always has a prime factorization  
p_1^k_1*...*p_n^k_n.

I merely want to write a as x^y where x is does not have a root, in  
order to get a unique representation. The data type used it in is  
rational numbers > 0 augmented with rational exponents. The  
application is for musical intervals - this is the reason it should be  
fast.

> A more tractable algorithm is the following:
>
> for y = floor(log(a)/log(2)) to 2 by -1
>  If a^(1/y) is an integer, return y
>
> The starting value is the largest power y for which a^(1/y) >= 2.
>
> Testing if a^(1/y) can be done efficienty by computing mod some small
> numbers; thus most non-solutions can be quickly rejected.

I am using the prime number decomposition in Haskell, which is very  
fast. But I get some very large primes. For example,
   pi = 2^-48*5*7*11159*2264103847    from pi =  
884279719003555/281474976710656
   e = 2^-51*3*29*283*1483*167639911  from e =  
6121026514868073/2251799813685248
It does not take that long time to compute pi.

But looking at your code, recurring in log time, it makes me wonder if  
it is faster.