broot.c revision 1.1.1.1
1/* mpn_broot -- Compute hensel sqrt 2 3 Contributed to the GNU project by Niels M�ller 4 5 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY 6 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 7 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. 8 9Copyright 2012 Free Software Foundation, Inc. 10 11This file is part of the GNU MP Library. 12 13The GNU MP Library is free software; you can redistribute it and/or modify 14it under the terms of the GNU Lesser General Public License as published by 15the Free Software Foundation; either version 3 of the License, or (at your 16option) any later version. 17 18The GNU MP Library is distributed in the hope that it will be useful, but 19WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 20or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 21License for more details. 22 23You should have received a copy of the GNU Lesser General Public License 24along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 25 26#include "gmp.h" 27#include "gmp-impl.h" 28 29/* Computes a^e (mod B). Uses right-to-left binary algorithm, since 30 typical use will have e small. */ 31static mp_limb_t 32powlimb (mp_limb_t a, mp_limb_t e) 33{ 34 mp_limb_t r = 1; 35 mp_limb_t s = a; 36 37 for (r = 1, s = a; e > 0; e >>= 1, s *= s) 38 if (e & 1) 39 r *= s; 40 41 return r; 42} 43 44/* Computes a^{1/k - 1} (mod B^n). Both a and k must be odd. 45 46 Iterates 47 48 r' <-- r - r * (a^{k-1} r^k - 1) / n 49 50 If 51 52 a^{k-1} r^k = 1 (mod 2^m), 53 54 then 55 56 a^{k-1} r'^k = 1 (mod 2^{2m}), 57 58 Compute the update term as 59 60 r' = r - (a^{k-1} r^{k+1} - r) / k 61 62 where we still have cancelation of low limbs. 63 64 */ 65void 66mpn_broot_invm1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t k) 67{ 68 mp_size_t sizes[GMP_LIMB_BITS * 2]; 69 mp_ptr akm1, tp, rnp, ep, scratch; 70 mp_limb_t a0, r0, km1, kp1h, kinv; 71 mp_size_t rn; 72 unsigned i; 73 74 TMP_DECL; 75 76 ASSERT (n > 0); 77 ASSERT (ap[0] & 1); 78 ASSERT (k & 1); 79 ASSERT (k >= 3); 80 81 TMP_MARK; 82 83 akm1 = TMP_ALLOC_LIMBS (4*n); 84 tp = akm1 + n; 85 86 km1 = k-1; 87 /* FIXME: Could arrange the iteration so we don't need to compute 88 this up front, computing a^{k-1} * r^k as (a r)^{k-1} * r. Note 89 that we can use wraparound also for a*r, since the low half is 90 unchanged from the previous iteration. Or possibly mulmid. Also, 91 a r = a^{1/k}, so we get that value too, for free? */ 92 mpn_powlo (akm1, ap, &km1, 1, n, tp); /* 3 n scratch space */ 93 94 a0 = ap[0]; 95 binvert_limb (kinv, k); 96 97 /* 4 bits: a^{1/k - 1} (mod 16): 98 99 a % 8 100 1 3 5 7 101 k%4 +------- 102 1 |1 1 1 1 103 3 |1 9 9 1 104 */ 105 r0 = 1 + (((k << 2) & ((a0 << 1) ^ (a0 << 2))) & 8); 106 r0 = kinv * r0 * (k+1 - akm1[0] * powlimb (r0, k & 0x7f)); /* 8 bits */ 107 r0 = kinv * r0 * (k+1 - akm1[0] * powlimb (r0, k & 0x7fff)); /* 16 bits */ 108 r0 = kinv * r0 * (k+1 - akm1[0] * powlimb (r0, k)); /* 32 bits */ 109#if GMP_NUMB_BITS > 32 110 { 111 unsigned prec = 32; 112 do 113 { 114 r0 = kinv * r0 * (k+1 - akm1[0] * powlimb (r0, k)); 115 prec *= 2; 116 } 117 while (prec < GMP_NUMB_BITS); 118 } 119#endif 120 121 rp[0] = r0; 122 if (n == 1) 123 { 124 TMP_FREE; 125 return; 126 } 127 128 /* For odd k, (k+1)/2 = k/2+1, and the latter avoids overflow. */ 129 kp1h = k/2 + 1; 130 131 /* FIXME: Special case for two limb iteration. */ 132 rnp = TMP_ALLOC_LIMBS (2*n + 1); 133 ep = rnp + n; 134 135 /* FIXME: Possible to this on the fly with some bit fiddling. */ 136 for (i = 0; n > 1; n = (n + 1)/2) 137 sizes[i++] = n; 138 139 rn = 1; 140 141 while (i-- > 0) 142 { 143 /* Compute x^{k+1}. */ 144 mpn_sqr (ep, rp, rn); /* For odd n, writes n+1 limbs in the 145 final iteration.*/ 146 mpn_powlo (rnp, ep, &kp1h, 1, sizes[i], tp); 147 148 /* Multiply by a^{k-1}. Can use wraparound; low part equals 149 r. */ 150 151 mpn_mullo_n (ep, rnp, akm1, sizes[i]); 152 ASSERT (mpn_cmp (ep, rp, rn) == 0); 153 154 ASSERT (sizes[i] <= 2*rn); 155 mpn_pi1_bdiv_q_1 (rp + rn, ep + rn, sizes[i] - rn, k, kinv, 0); 156 mpn_neg (rp + rn, rp + rn, sizes[i] - rn); 157 rn = sizes[i]; 158 } 159 TMP_FREE; 160} 161 162/* Computes a^{1/k} (mod B^n). Both a and k must be odd. */ 163void 164mpn_broot (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t k) 165{ 166 mp_ptr tp; 167 TMP_DECL; 168 169 ASSERT (n > 0); 170 ASSERT (ap[0] & 1); 171 ASSERT (k & 1); 172 173 if (k == 1) 174 { 175 MPN_COPY (rp, ap, n); 176 return; 177 } 178 179 TMP_MARK; 180 tp = TMP_ALLOC_LIMBS (n); 181 182 mpn_broot_invm1 (tp, ap, n, k); 183 mpn_mullo_n (rp, tp, ap, n); 184 185 TMP_FREE; 186} 187