1 1.1 mrg /* mpn_bsqrtinv, compute r such that r^2 * y = 1 (mod 2^{b+1}). 2 1.1 mrg 3 1.1 mrg Contributed to the GNU project by Martin Boij (as part of perfpow.c). 4 1.1 mrg 5 1.1.1.2 mrg Copyright 2009, 2010, 2012, 2015 Free Software Foundation, Inc. 6 1.1 mrg 7 1.1 mrg This file is part of the GNU MP Library. 8 1.1 mrg 9 1.1 mrg The GNU MP Library is free software; you can redistribute it and/or modify 10 1.1.1.2 mrg it under the terms of either: 11 1.1.1.2 mrg 12 1.1.1.2 mrg * the GNU Lesser General Public License as published by the Free 13 1.1.1.2 mrg Software Foundation; either version 3 of the License, or (at your 14 1.1.1.2 mrg option) any later version. 15 1.1.1.2 mrg 16 1.1.1.2 mrg or 17 1.1.1.2 mrg 18 1.1.1.2 mrg * the GNU General Public License as published by the Free Software 19 1.1.1.2 mrg Foundation; either version 2 of the License, or (at your option) any 20 1.1.1.2 mrg later version. 21 1.1.1.2 mrg 22 1.1.1.2 mrg or both in parallel, as here. 23 1.1 mrg 24 1.1 mrg The GNU MP Library is distributed in the hope that it will be useful, but 25 1.1 mrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 26 1.1.1.2 mrg or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 27 1.1.1.2 mrg for more details. 28 1.1 mrg 29 1.1.1.2 mrg You should have received copies of the GNU General Public License and the 30 1.1.1.2 mrg GNU Lesser General Public License along with the GNU MP Library. If not, 31 1.1.1.2 mrg see https://www.gnu.org/licenses/. */ 32 1.1 mrg 33 1.1 mrg #include "gmp-impl.h" 34 1.1 mrg 35 1.1 mrg /* Compute r such that r^2 * y = 1 (mod 2^{b+1}). 36 1.1 mrg Return non-zero if such an integer r exists. 37 1.1 mrg 38 1.1 mrg Iterates 39 1.1 mrg r' <-- (3r - r^3 y) / 2 40 1.1 mrg using Hensel lifting. Since we divide by two, the Hensel lifting is 41 1.1 mrg somewhat degenerates. Therefore, we lift from 2^b to 2^{b+1}-1. 42 1.1 mrg 43 1.1 mrg FIXME: 44 1.1 mrg (1) Simplify to do precision book-keeping in limbs rather than bits. 45 1.1 mrg 46 1.1 mrg (2) Rewrite iteration as 47 1.1 mrg r' <-- r - r (r^2 y - 1) / 2 48 1.1 mrg and take advantage of zero low part of r^2 y - 1. 49 1.1 mrg 50 1.1 mrg (3) Use wrap-around trick. 51 1.1 mrg 52 1.1 mrg (4) Use a small table to get starting value. 53 1.1 mrg */ 54 1.1 mrg int 55 1.1 mrg mpn_bsqrtinv (mp_ptr rp, mp_srcptr yp, mp_bitcnt_t bnb, mp_ptr tp) 56 1.1 mrg { 57 1.1.1.2 mrg mp_ptr tp2; 58 1.1 mrg mp_size_t bn, order[GMP_LIMB_BITS + 1]; 59 1.1 mrg int i, d; 60 1.1 mrg 61 1.1 mrg ASSERT (bnb > 0); 62 1.1 mrg 63 1.1 mrg bn = 1 + bnb / GMP_LIMB_BITS; 64 1.1 mrg 65 1.1 mrg tp2 = tp + bn; 66 1.1 mrg 67 1.1 mrg rp[0] = 1; 68 1.1 mrg if (bnb == 1) 69 1.1 mrg { 70 1.1 mrg if ((yp[0] & 3) != 1) 71 1.1 mrg return 0; 72 1.1 mrg } 73 1.1 mrg else 74 1.1 mrg { 75 1.1 mrg if ((yp[0] & 7) != 1) 76 1.1 mrg return 0; 77 1.1 mrg 78 1.1 mrg d = 0; 79 1.1 mrg for (; bnb != 2; bnb = (bnb + 2) >> 1) 80 1.1 mrg order[d++] = bnb; 81 1.1 mrg 82 1.1 mrg for (i = d - 1; i >= 0; i--) 83 1.1 mrg { 84 1.1 mrg bnb = order[i]; 85 1.1 mrg bn = 1 + bnb / GMP_LIMB_BITS; 86 1.1 mrg 87 1.1.1.2 mrg mpn_sqrlo (tp, rp, bn); 88 1.1.1.2 mrg mpn_mullo_n (tp2, rp, tp, bn); /* tp2 <- rp ^ 3 */ 89 1.1.1.2 mrg 90 1.1.1.2 mrg mpn_mul_1 (tp, rp, bn, 3); 91 1.1 mrg 92 1.1 mrg mpn_mullo_n (rp, yp, tp2, bn); 93 1.1 mrg 94 1.1 mrg #if HAVE_NATIVE_mpn_rsh1sub_n 95 1.1 mrg mpn_rsh1sub_n (rp, tp, rp, bn); 96 1.1 mrg #else 97 1.1 mrg mpn_sub_n (tp2, tp, rp, bn); 98 1.1 mrg mpn_rshift (rp, tp2, bn, 1); 99 1.1 mrg #endif 100 1.1 mrg } 101 1.1 mrg } 102 1.1 mrg return 1; 103 1.1 mrg } 104
Indexes created Sat Apr 18 00:22:28 UTC 2026