perfect power / nth root improvements

[Torbjorn Granlund]

arthur <arthur at ilvsun3.com> writes:

  Torbjorn Granlund wrote:
  > The current code is O(M(n)(log n)^2),

  > it is possible to do O(M(n)) in
  > the worst case, or O(M(n/k)) on average, for computing the kth root of
  > an n-bit number.

  Is this the method outlined by Bernstein in: Detecting Perfect Powers
  In Essentially Linear Time?

I am talking about *computing* the truncated nth root.  In LaTeX the
function that mpn_root should compute can be described as
$\lfloor\sqrt[n]{A}\rfloor$.

For detecting perfect powers, Bernstein's algorithm is probably better
than what we have in GMP today, even with the commented-out padic
calls.  But improving the GMP nth root code could improve GMP's
perfect power recognition by a large factor, and of course help the
direct root calls.

My claim is that the kth root of a n-bit number A essentially depends
on just the n/k most significant bits of A.  An algorithm that
computes a result based on just these most significant bits, and
detects if the result is too small (by at most 1) quickly, should be
possible to design.

Detecting if the result is wrong might not be hard.  If we call the
kth root of A for B, we could compute (B+1)^k mod R, where R is the
architecture imposed word base (typically 2^32 or 2^64).  Then compare
the result to A (the least significant bits).

  In my sadness, I got busy and implemented a padic root routine for
  perfect power testing that runs 20% faster than the current code up to
  inputs around 10^100, then it really improves.

Cool!

Could you define, in mathematical terms, exactly what you've
implemented?

  Seconds to check numbers for being a perfect power.

  starting   check   how many  seconds  seconds
  at         this    perfect   w/ orig  w/ new
             many    powers    code     code
  10^0    10,000,000   3395     8.94     7.28
  10^10   10,000,000     51    11.11     8.17
  10^100  10,000,000      1    43.98    37.23
  10^1000    100,000      1    16.04     3.81
  10^10000    10,000      1   385.76     7.79

Excellent numbers!

Does the code avoid calls to mpz_root/mpn_root?

  I'm not familiar with the process of making code submissions.  I put
  the code and instructions here: www.ilvsun3.com/files/artroot.tar

The process is:

(1) finish the code
(2) assign the code to the FSF and get your employer to sign a disclaimer
    http://www.gnu.org/licenses/why-assign.html
    http://www.gnu.org/prep/maintain/html_node/Copyright-Papers.html

  It could easily be improved by implementing the root algorithm with
  GMP internal calls instead of calls through the high level interface
  (and better signed "mod" for mod powers-of-2).

There are several calls at the mpz level: mpz_cdiv_r_2exp,
mpz_fdiv_r_2exp, mpz_tdiv_r_2exp.

The mpn level is being expanded with many new functions for GMP 5, if
you want to implement at the mon level your code you might want to use
these calls.  (GMP 5 is still a few years away, but I publish things
at gmplib.org/devel as they get ready.)

  I tested it against the current algorithm for the range of integers
  -2^30..2^30 and various large numbers with identical results.  I don't
  know if this is the best nth root algorithm to use as a long-term
  choice for GMP, but it seems to be working quickly and correctly.

We don't need to have optimal code, a great improvement is always
welcome.

[I am busy this week and the next week with exams, please accept slow
response from me.]

-- 
Torbjörn