Anomaly in mpn_sqrtrem and mpn_rottrem

[Niels Möller]

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

> An alternative would be to keep the current
> algorithm, but use the fact that the previous divisor is a damn good
> staring approximation to the new one, and use Newton to maintain an
> inverse to it.

I vaguely remember that idea. That would be a loop that maintains
approximations of both x \approx sqrt(a) and v \approx 1/sqrt(a). It
makes my head spin a bit to think about the precision and order of the
needed updates.

When we have a n bit approximation x,
we need a 2n bit inverse v of x in order to compute the new x' with 2n
bits of precision. For the next iteration, we then need a 4n-bit inverse
v' of x'.

As inputs for the v' update we have the old v (2n bits), the new square
root approximation x' (2n bits), and we also have the old x (or the
difference x' - x) available.

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 also had a quick look at the math, and I realized (some of you surely
knew that already) that floor(sqrt(a)) is mostly independent of the
lowest half of the bits of a. So its essentially an operation with n
input bits and n output bits. Proof:


  sqrt(a+d) - sqrt(a) = d / (sqrt(a+d) + sqrt(a))

So if d <= sqrt(a), then

  sqrt(a+d) - sqrt(a) < 1/2.

The lowest half of the bits are needed only for computing the remainder,
and a final adjustment-by-one step. I think the iteration should totally
ignore the low half bits in each step.

Would it make sense with an mpn_sqrt_appr which takes a n-limb
normalized number A as input, and computes an n-limb square root
approximation of B^{n-1} A, with some well defined error bound, of at
most 2? (We might also need a function for the the (n-1)-limb sqrt of
B^{n-2} A).

And a corresponding mpn_root_appr would compute the kth root of
B^{(k-1) n - 1} A, I think.

Regards,
/Niels

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