* Niels Möller <ni...@lysator.liu.se> [Jul 20. 2012 08:20]: > ni...@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@gmplib.org > http://gmplib.org/mailman/listinfo/gmp-devel _______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org http://gmplib.org/mailman/listinfo/gmp-devel