GMP mpz_root() improvement (fwd)

[Ivo Moravec]

Forwarding message to  entire list.

So, when I send the code, where do I send it?  (or do you want a URL?)

Ivo

---------- Forwarded message ----------
Date: Fri, 23 Jan 2004 10:24:26 -0500 (EST)
From: Ivo Moravec <imoravec at scg.math.uwaterloo.ca>
To: Kevin Ryde <user42 at zip.com.au>
Subject: Re: GMP mpz_root() improvement



On Thu, 22 Jan 2004, Kevin Ryde wrote:

> Ivo Moravec <imoravec at scg.math.uwaterloo.ca> writes:
> >
> > Halley Iterations (a third order) for computing integer nth roots.
>
> I think I'd seen that (and higher order forms) at one time, but wasn't
> smart enough to know if it would be an overall saving, since of course
> each iteration does more work.

Well, to tell the truth, the difference between Halley an Newton
iterations are not THAT great - at least in my tests.  I modified the GMP
implementation to work at partial precision, and Halley did come out
ahead, but not by a huge amount.  The difference would really start
to really appear with numbers many tens of thousands of bits and up.

> > Basically, I've written a replacement for
> > your mpz_root() and mpn_rootrem() functions,
>
> Good, fast, clean, well-tested code is always of interest.

I'll clean it up and send it - give me a couple days to strip out the
debug code and all that.

> > The new code runs about 10 to 30 or more times faster than the
> > current implementation in gmp-4.1
>
> The main problem with the current code, as noted in the projects file,
> is that it uses full precision at every step.
>
> A comparison with sqrt might be more indicative.

A comparison against sqrt is not really fair either.  For one thing, there
is no exponentiation in a sqrt root, and second you use look up tables for
your sqrt function which are not possible for nth roots.  As a result, I
just call the mpn sqrt function from the nth root code if a sqrt is
requested.  What I can tell you is that the code is faster than the
current Newton iteration modified to work at partial precision.

I need to know your policy on two things however:

1)
mpn_rootrem() calculates the residual of the root (a-x^n), which
mpz_root() is very happy to throw away.  I have made 2 routines.  One
which computes the residual, and one that doesn't.  The one that doesn't
is about twice as quick as the one that does.  How willing would you be to
adding a "mpn_root" function which would enable the mpz_root() function to
run much faster?

2)
I played around using floating point numbers to get an inital estimate for
use in the refinement.  It worked quite well since my machine is an Intel
based processor which has hardware log and exp instructions.  It's my
understanding that not all hardware does have them however.  What's your
policy on this?  I have a regular integer based guesstimator too, so I can
use that if the floating point is not suitable.