2-adic roots (Re: bdiv vs redc)

[Niels Möller]

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

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

Yet another correction, after more work on the implementation.

For numbers x = 1 (mod 8), there are four square roots mod 2^k. If x is
one of them, the other three are

  x + 2^{k-1}
  -x
  -x + 2^{k-1}

So the top bit doesn't matter. This implies that no extra bit is needed
for temporary values: If we start with \ell bits, we can do the
iteration using 2\ell-2 bits for all values. And it doesn't matter
which value we get in the most significant bit after the final shift
right, we get a square root (mod 2^{2\ell - 2}) in either case.

And we can canonicalize the returned (mod 2^k) square root by saying
that it must be = 1 (mod 4), and < 2^{k-1} (i.e., most significant bit
always zero).

I'm attaching my current code for both square root and nth root. Appears
to work fine. None of them yet use any wraparound tricks, but they do
take some advantage from cancellation.

In the nth-root code, I use the iteration 

     x' <-- x - x * (a^{n-1} x^n - 1) / n

which converges to a^{1/n-1}, so the nth root is recovered as a * x. The
number of correct bits is *exactly* doubled in each iteration, so I
didn't see much point in using bit counts, instead I use a limb count
for the desired precision. Perhaps this iteration could be useful in the
euclidean case too, I haven't investigated that.

The factor a^{n-1} is a loop invariant, so it can be computed (at the
highest needed precision) before the loop.

For the power x^n, I currently use mpn_powlo. But the least significant
half is known from the previous iteration, so wraparound would be
desirable. To me it would make some sense with a pow function for (mod
(2^k - 1)), i.e., using mpn_mulmod_bnm1 and mpn_sqrmod_bnm1. Or is there
any easier way to take advantage of wraparound for that computation?

Maybe it would be a good idea to write a general or "abstract" pow
function taking multiplication and squaring functions as arguments?

We currently have modular exponentation, powlo and "regular" powering
with no reduction of any kind. I'm suggesting a pow_modbnm1. For
euclidean square root, and for mpfr, it might also be useful with a
pow_high, keeping only the n most significant limbs of each product, and
returning the number of discarded low limbs. Another potential use is
powm with moduli of special form, where the reduction can be done
cheaper than with montgomery redc. Maybe the function could even be used
for more complicated groups, e.g., if implementing elliptic curve
operations on top of gmp. Or maybe it's easy enough to duplicate the
code for each useful case, perhaps sharing some bit extraction macros in
gmp-impl.h.

Regards,
/Niels

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: broot.c
URL: <http://gmplib.org/list-archives/gmp-devel/attachments/20120730/f94c8e89/attachment-0002.c>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: bsqrt.c
URL: <http://gmplib.org/list-archives/gmp-devel/attachments/20120730/f94c8e89/attachment-0003.c>
-------------- next part --------------

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