Torbjorn Granlund <t...@gmplib.org> writes: > Used for computing x^3 mod B^n. It is really a very plain Newton > iteration being used here.
I see. I think the "right" way to compute the iteration x <-- x - x (a x^2 - 1) / 2 is as follows. Let the current x have \ell bits. 1. Compute the square x^2. Use wraparound, since the lowest \ell bits were computed in the same step in the previous iteration. 2. Compute a x^2. Use wraparound, since the low \ell bits are 0, ...,0, 1. Shift right, and ignore the low zero limbs. Denote the result as B^z e = (a x^2 - 1) / 2 3. Compute x e (mod 2^{2\ell} / B^z]). This is balanced (almost) plain mullo. 4. Subtract x -= B^z x e, extending the precision of x from \ell bits to 2\ell - 2. We could simplify the final subtraction (which is mostly a negation) if we negate a up front, really using a newton iteration converging to (-a)^{-1/2} (mod increasing powers of two). > I suppose one should for common k improve the starting value from 1 bit > to a few bits, and for any k iterate in mp_limb_t variables until > getting a full word of precision (using a fully unrolled loop). For square root I think the most practical is to tabulate square roots mod 2^10 (using a table with only 2^7 entries), and then iterate 10 -> 18 -> 34 -> 66. I don't think it is of much use to use a larger table unless we go up to 2^14 or 2^15 entries. I also think tabulating small kth roots for arbitrary (odd) k is practical, since only the low bits of k matter, but I haven't yet looked into the details. 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