perfect power / nth root improvements

[Niels Möller]

arthur <arthur at ilvsun3.com> writes:

> Benefits of p-adic vs. floating point iterations:
> a) the exact number of iterations needed can be calculated precisely
> b) the number of p-adic digits of accuracy exactly doubles with each 
> iteration
> c) no divisions
> d) power of 2 mods can be implemented efficiently by bit manipulations
>
>
>
> Let me explain how I got here.

Do I understand the overall power testing algorithm correctly as
follows?

Input: Number n. Output: x, k, such that x^k = n and k >= 2, if such
numbers exist.

1. For each potential exponent k, 2 <= k <= log_2(n) (I guess it's
   sufficient to iterate over primes in that interval, since in the
   relation x^k = n with the smallest k >= 2, k is prime. If we really
   returned x and k, it would make some sense to find the largest k,
   but for power *testing*, that doesn't matter).

2. Do some quick test on the least significant word. If k can't
   possibly be an exponent, continue with the next one.

3. Compute a candidate k:th root r of n using 2-adic newton iteration.

4. Do some quick test based on the most significant words of r and n.

5. Actually compute r^k, and check if it equals n.

One trick you don't mention, but which might be useful, is to first
check for small prime factors of n. The benefits are that

* If there are no factors less than B, then x > B, which implies a
  better upper bound on k.

* If there are any small prime factors, that implies that x >= their
  product, also giving an improved upper bound.

* I'm not sure if this pays off, but if one also divides out all the
  small factors, then the remaining number corresponds to a base > B,
  giving an upper bound on k. Furthermore, k must divide the
  multiplicity of each of the small primes that were divided out. Let
  g = gcd(small prime multiplicities), then k must be a factor of g,
  and in particular if g = 1, then n is not a power.

The prime power detection in Cohen (computational algebraic number
theory, algorithm 1.7.5) might also give some inspiration, but I don't
see any easy way to extend it to non-prime powers (one of the key
ideas seem to be that for a prime power, any witness to compositeness
in the the miller-rabin test also produces a non-trivial factor, which
isn't true at all in general...)

> Suggestions, regardless of what is done with this p-adic routine:
>
> 1) provide efficient mpz_mod and mpz_powm for powers of 2, similar to 
> the bit shifting calls mpz_tdiv_r_2exp.

Better powm is in the works, I think. The idea is to write the modulo
n as 2^k m', with odd m', then do the exponentiations separately mod
2^k and m', and combine using the chinese remainder theorem.

> 2) provide signed mod routines as maple's "mods".  It obviously has a 
> use, as this routine needs it, and I had to cobble together a crude 
> solution.  Here's a definition of mods:

That might be useful, but I don't understand why you would need that
for the 2-adic newton iteration.

> 4) see about replacing bit-test with a more precise test, looking at the 
> trade-off in time of how long the test takes vs. the computation time 
> required for the candidates it excludes over and above the bit-test

I think you could du something reasonable by extracting the most
significant (normalized) word of r, exponentiate using ordinary
square-and-multiply (not sure if sliding window applies easily, maybe
it does) where for each multiplication, you keep the most significant
(normalized) word of the product. The output from the exponentation is
a single word and a shift count, to be compared to n. Analyzing the
truncation error, and/or doing something with interval arithmetic,
shouldn't be too hard, and my guess is that you would get sufficient
accuracy to rule out many candidate roots.

But using floating point (if you can manage the accuracy bounds) may
well be more efficient.

Regards,
/Niels