2-adic roots (Re: bdiv vs redc)

[Niels Möller]

nisse at lysator.liu.se (Niels Möller) writes:

> 3. A reasonable square root algorithm is the standard Newton iteration
>    x_j which converges to a^{-1/2}, each step extending from \ell bits
>    of precision to 2 \ell - 1.

Correction: Each iteration gets from \ell to 2 \ell-2, and it needs 2
\ell-1 bits of precision for intermediate values.

  x'  <-- x - x * (a x^2 - 1)/ 2

I've implemented this, see attached code. One iteration needs one
squaring and two multiplies. The first operation could make good use of
wraparound multiplication since the low half is known in advance (but
the code doesn't do that), the second could also use wraparound
arithmetic, or possibly mulmid, while the third is a mullo but with many
trailing zeros on one of the operands (the code tries to take advantage
of that.

The code returns the square root which is 1 (mod 8), out of the four
possible.

Happpy hacking,
/Niels

-------------- next part --------------
A non-text attachment was scrubbed...
Name: bsqrt.c
Type: application/octet-stream
Size: 5961 bytes
Desc: not available
URL: <http://gmplib.org/list-archives/gmp-devel/attachments/20120723/99c0937c/attachment.obj>
-------------- next part --------------


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