Anomaly in mpn_sqrtrem and mpn_rottrem

[Niels Möller]

tg at gmplib.org (Torbjörn Granlund) writes:

> nisse at lysator.liu.se (Niels Möller) writes:
>
>   As far as I see, a plain Newton iteration on won't produce an x' inverse
>   with 4n bits of accuracy. One would need aither two iterations, or some
>   other trick, maybe mixing in some interpolation of x' - x.
>   
> I don't think one should struggle in order to get n -> 2n bits in a
> Newtonian iteration.  It is enough to achieve n -> 2n-k for some
> constant k.

I don't think one can do that for any small k in this case. I'll try to
describe the problem.

In the Newton iteration for the square root, which extends the precision
of x from n bits to (close to) 2n bits, I'd expect that we 2n bits of
precision for all inputs, including the inverse of the input x. Do you
agree? So as input, we need a 2n-bit approximative inverse of the n-bit
input x.

As output, we want the new 2n-bit x and a corresponding 4n-bit
approximative inverse.

A single newton iteration can produce a 4n-bit approximation of the
*n*-bit x, but then we also have new x bits from the square root
approximation to take into account. One way to get the needed inverse
approximation is to

1. Adjust the 2n-bit inverse of n-bit x to a 2n-bit inverse of the new
   2n-bit x. Could be done with a Newton step (which seems a bit
   overkill), or perhaps one can use some cheaper linear interpolation.

2. Do a Newton step to extend that inverse of the 2n-bit x to 4n
   bits of precision.

Hmm. Or maybe this is stupid. I could stop insisting on using a full
size inverse (so that A / x or E / x can be computed as a *single*
mulhi), and instead work with a half-size inverse, so that the quotient
is computed in two steps. Then inverse computation using a plain
Newton-step per iteration should work fine.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.