1 /* mpn_bsqrtinv, compute r such that r^2 * y = 1 (mod 2^{b+1}).
3 Contributed to the GNU project by Martin Boij (as part of perfpow.c).
5 Copyright 2009, 2010, 2012 Free Software Foundation, Inc.
7 This file is part of the GNU MP Library.
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MP Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
25 /* Compute r such that r^2 * y = 1 (mod 2^{b+1}).
26 Return non-zero if such an integer r exists.
29 r' <-- (3r - r^3 y) / 2
30 using Hensel lifting. Since we divide by two, the Hensel lifting is
31 somewhat degenerates. Therefore, we lift from 2^b to 2^{b+1}-1.
34 (1) Simplify to do precision book-keeping in limbs rather than bits.
36 (2) Rewrite iteration as
37 r' <-- r - r (r^2 y - 1) / 2
38 and take advantage of zero low part of r^2 y - 1.
40 (3) Use wrap-around trick.
42 (4) Use a small table to get starting value.
45 mpn_bsqrtinv (mp_ptr rp, mp_srcptr yp, mp_bitcnt_t bnb, mp_ptr tp)
49 mp_size_t bn, order[GMP_LIMB_BITS + 1];
54 bn = 1 + bnb / GMP_LIMB_BITS;
72 for (; bnb != 2; bnb = (bnb + 2) >> 1)
75 for (i = d - 1; i >= 0; i--)
78 bn = 1 + bnb / GMP_LIMB_BITS;
80 mpn_mul_1 (tp, rp, bn, k);
82 mpn_powlo (tp2, rp, &k, 1, bn, tp3);
83 mpn_mullo_n (rp, yp, tp2, bn);
85 #if HAVE_NATIVE_mpn_rsh1sub_n
86 mpn_rsh1sub_n (rp, tp, rp, bn);
88 mpn_sub_n (tp2, tp, rp, bn);
89 mpn_rshift (rp, tp2, bn, 1);