Integer root extracting
Hans Aberg
haberg-1 at telia.com
Fri Sep 3 16:09:19 CEST 2010
On 3 Sep 2010, at 15:57, Torbjorn Granlund wrote:
>> 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 thought you were interested in a usable algorithm, not merely a
> theoretically correct method. I meant thatt complete factorisation of
> large integers is in practice not possible, since it would take too
> much
> time.
Yes, that is the reason I was asking the question in the first place.
Even though small primes would be used, some large ones may come in
when converting from floats (which, in turn, is a bit corny).
>> 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.
>
> I am quite sure it is faster even for small numbers, and certainly so
> for large numbers. An added feature of my pseudo code is that its
> worst-case performance isn't terrible, while factoring takes wildly
> different time for different numbers.
That is another concern, but mostly, primes need not to be recomputed,
as addition is not used.
Anyway, I made a modification of your algorithm restricting to primes
and iterating mpz_root(), which you perhaps not have seen yet.
More information about the gmp-discuss
mailing list