mpn_sqrtrem1

[paul zimmermann]

       Torbjörn,

> I am not proud of the mpn_sqrtrem1 code.  It was made with some analysis
> and lots of testing.  The many undocumented magic constants are ugly.
> Perhaps I should never have checked in this code.

ah ok. I have an alternate implementation which is almost as efficient,
but although it passes make check with both GMP and MPFR, I have not yet
proven it is correct :-(

> Furthermore, for a real performance win one would need to rewrite
> mpn_sqrtrem2.

I know how to improve mpn_sqrtrem2. For a 64-bit computer, start from a 32-bit
approximation, say x, of 1/sqrt(a), which is already done by mpn_sqrtrem1. Then
use Karp-Markstein's trick for the square root:

y = approx(a*x) # 32-bit approximation of sqrt(a)
t = a - y^2     
return y + x/2*t # 64-bit approximation of sqrt(a)

I will try to implement that.

Paul