Integer root extracting
Hans Aberg
haberg-1 at telia.com
Fri Sep 3 12:46:02 CEST 2010
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.
More information about the gmp-discuss
mailing list