2-adic roots (Re: bdiv vs redc)

[Joerg Arndt]

* Niels Möller <nisse at lysator.liu.se> [Jul 20. 2012 08:20]:
> nisse at lysator.liu.se (Niels Möller) writes:
> 
> > I haven't really looked into the perfpow code or theory, but let me
> > think aloud...
> 
> I've looked a bit more into this. This is my current understanding:
> 
> 1. Powers of two have to be handled specially, so consider only odd
>    numbers, we'll be working with the multiplicative group Z_{2^k}^*.
> 
> 2. An odd number has a square root (or in fact two) if and only if it's
>    = 1 (mod 4).

IIRC that should be = 1 (mod 8)

Here are my routines for 2-adic sqrt() and inverse sqrt()
(for type unsigned long):

static inline ulong invsqrt2adic(ulong d)
// Return inverse square root modulo 2**BITS_PER_LONG
// Must have:  d==1 mod 8
// The number of correct bits is doubled with each step
// ==> loop is executed prop. log_2(BITS_PER_LONG) times
// precision is 4, 8, 16, 32, 64, ... bits (or better)
{
    if ( 1 != (d&7) )  return 0;  // no inverse sqrt
    // start value: if d == ****10001 ==> x := ****1001
    ulong x = (d >> 1) | 1;
    ulong p, y;
    do
    {
        y = x;
        p = (3 - d * y * y);
        x = (y * p) >> 1;
    }
    while ( x!=y );
    return  x;
}
// -------------------------


static inline ulong sqrt2adic(ulong d)
// Return square root modulo 2**BITS_PER_LONG
// Must have: d==1 mod 8  or  d==4 mod 32,  d==16 mod 128
//   ... d==4**k mod 4**(k+3)
// Result undefined if condition does not hold
{
    if ( 0==d )  return 0;
    ulong s = 0;
    while ( 0==(d&1) )  { d >>= 1; ++s; }
    d *= invsqrt2adic(d);
    d <<= (s>>1);
    return   d;
}
// -------------------------

For GMP an iteration along the lines of Schoenhage's coupled method
may be the best way to proceed.


> [...]

> 
> Does this seem right?
> 
> Regards,
> /Niels
> 
> -- 
> Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
> Internet email is subject to wholesale government surveillance.
> _______________________________________________
> gmp-devel mailing list
> gmp-devel at gmplib.org
> https://gmplib.org/mailman/listinfo/gmp-devel