1 /* mpz_powm(res,base,exp,mod) -- Set R to (U^E) mod M. 2 3 Contributed to the GNU project by Torbjorn Granlund. 4 5 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2005, 2008, 6 2009, 2011, 2012 Free Software Foundation, Inc. 7 8 This file is part of the GNU MP Library. 9 10 The GNU MP Library is free software; you can redistribute it and/or modify 11 it under the terms of the GNU Lesser General Public License as published by 12 the Free Software Foundation; either version 3 of the License, or (at your 13 option) any later version. 14 15 The GNU MP Library is distributed in the hope that it will be useful, but 16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 17 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 18 License for more details. 19 20 You should have received a copy of the GNU Lesser General Public License 21 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 22 23 24 #include "gmp.h" 25 #include "gmp-impl.h" 26 #include "longlong.h" 27 28 29 /* TODO 30 31 * Improve handling of buffers. It is pretty ugly now. 32 33 * For even moduli, we compute a binvert of its odd part both here and in 34 mpn_powm. How can we avoid this recomputation? 35 */ 36 37 /* 38 b ^ e mod m res 39 0 0 0 ? 40 0 e 0 ? 41 0 0 m ? 42 0 e m 0 43 b 0 0 ? 44 b e 0 ? 45 b 0 m 1 mod m 46 b e m b^e mod m 47 */ 48 49 #define HANDLE_NEGATIVE_EXPONENT 1 50 51 void 52 mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m) 53 { 54 mp_size_t n, nodd, ncnt; 55 int cnt; 56 mp_ptr rp, tp; 57 mp_srcptr bp, ep, mp; 58 mp_size_t rn, bn, es, en, itch; 59 mpz_t new_b; /* note: value lives long via 'b' */ 60 TMP_DECL; 61 62 n = ABSIZ(m); 63 if (UNLIKELY (n == 0)) 64 DIVIDE_BY_ZERO; 65 66 mp = PTR(m); 67 68 TMP_MARK; 69 70 es = SIZ(e); 71 if (UNLIKELY (es <= 0)) 72 { 73 if (es == 0) 74 { 75 /* b^0 mod m, b is anything and m is non-zero. 76 Result is 1 mod m, i.e., 1 or 0 depending on if m = 1. */ 77 SIZ(r) = n != 1 || mp[0] != 1; 78 PTR(r)[0] = 1; 79 TMP_FREE; /* we haven't really allocated anything here */ 80 return; 81 } 82 #if HANDLE_NEGATIVE_EXPONENT 83 MPZ_TMP_INIT (new_b, n + 1); 84 85 if (UNLIKELY (! mpz_invert (new_b, b, m))) 86 DIVIDE_BY_ZERO; 87 b = new_b; 88 es = -es; 89 #else 90 DIVIDE_BY_ZERO; 91 #endif 92 } 93 en = es; 94 95 bn = ABSIZ(b); 96 97 if (UNLIKELY (bn == 0)) 98 { 99 SIZ(r) = 0; 100 TMP_FREE; 101 return; 102 } 103 104 ep = PTR(e); 105 106 /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case. */ 107 if (UNLIKELY (en == 1 && ep[0] == 1)) 108 { 109 rp = TMP_ALLOC_LIMBS (n); 110 bp = PTR(b); 111 if (bn >= n) 112 { 113 mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1); 114 mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n); 115 rn = n; 116 MPN_NORMALIZE (rp, rn); 117 118 if (SIZ(b) < 0 && rn != 0) 119 { 120 mpn_sub (rp, mp, n, rp, rn); 121 rn = n; 122 MPN_NORMALIZE (rp, rn); 123 } 124 } 125 else 126 { 127 if (SIZ(b) < 0) 128 { 129 mpn_sub (rp, mp, n, bp, bn); 130 rn = n; 131 rn -= (rp[rn - 1] == 0); 132 } 133 else 134 { 135 MPN_COPY (rp, bp, bn); 136 rn = bn; 137 } 138 } 139 goto ret; 140 } 141 142 /* Remove low zero limbs from M. This loop will terminate for correctly 143 represented mpz numbers. */ 144 ncnt = 0; 145 while (UNLIKELY (mp[0] == 0)) 146 { 147 mp++; 148 ncnt++; 149 } 150 nodd = n - ncnt; 151 cnt = 0; 152 if (mp[0] % 2 == 0) 153 { 154 mp_ptr newmp = TMP_ALLOC_LIMBS (nodd); 155 count_trailing_zeros (cnt, mp[0]); 156 mpn_rshift (newmp, mp, nodd, cnt); 157 nodd -= newmp[nodd - 1] == 0; 158 mp = newmp; 159 ncnt++; 160 } 161 162 if (ncnt != 0) 163 { 164 /* We will call both mpn_powm and mpn_powlo. */ 165 /* rp needs n, mpn_powlo needs 4n, the 2 mpn_binvert might need more */ 166 mp_size_t n_largest_binvert = MAX (ncnt, nodd); 167 mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert); 168 itch = 3 * n + MAX (itch_binvert, 2 * n); 169 } 170 else 171 { 172 /* We will call just mpn_powm. */ 173 mp_size_t itch_binvert = mpn_binvert_itch (nodd); 174 itch = n + MAX (itch_binvert, 2 * n); 175 } 176 tp = TMP_ALLOC_LIMBS (itch); 177 178 rp = tp; tp += n; 179 180 bp = PTR(b); 181 mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp); 182 183 rn = n; 184 185 if (ncnt != 0) 186 { 187 mp_ptr r2, xp, yp, odd_inv_2exp; 188 unsigned long t; 189 int bcnt; 190 191 if (bn < ncnt) 192 { 193 mp_ptr newbp = TMP_ALLOC_LIMBS (ncnt); 194 MPN_COPY (newbp, bp, bn); 195 MPN_ZERO (newbp + bn, ncnt - bn); 196 bp = newbp; 197 } 198 199 r2 = tp; 200 201 if (bp[0] % 2 == 0) 202 { 203 if (en > 1) 204 { 205 MPN_ZERO (r2, ncnt); 206 goto zero; 207 } 208 209 ASSERT (en == 1); 210 t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt; 211 212 /* Count number of low zero bits in B, up to 3. */ 213 bcnt = (0x1213 >> ((bp[0] & 7) << 1)) & 0x3; 214 /* Note that ep[0] * bcnt might overflow, but that just results 215 in a missed optimization. */ 216 if (ep[0] * bcnt >= t) 217 { 218 MPN_ZERO (r2, ncnt); 219 goto zero; 220 } 221 } 222 223 mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt); 224 225 zero: 226 if (nodd < ncnt) 227 { 228 mp_ptr newmp = TMP_ALLOC_LIMBS (ncnt); 229 MPN_COPY (newmp, mp, nodd); 230 MPN_ZERO (newmp + nodd, ncnt - nodd); 231 mp = newmp; 232 } 233 234 odd_inv_2exp = tp + n; 235 mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n); 236 237 mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd); 238 239 xp = tp + 2 * n; 240 mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt); 241 242 if (cnt != 0) 243 xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1; 244 245 yp = tp; 246 if (ncnt > nodd) 247 mpn_mul (yp, xp, ncnt, mp, nodd); 248 else 249 mpn_mul (yp, mp, nodd, xp, ncnt); 250 251 mpn_add (rp, yp, n, rp, nodd); 252 253 ASSERT (nodd + ncnt >= n); 254 ASSERT (nodd + ncnt <= n + 1); 255 } 256 257 MPN_NORMALIZE (rp, rn); 258 259 if ((ep[0] & 1) && SIZ(b) < 0 && rn != 0) 260 { 261 mpn_sub (rp, PTR(m), n, rp, rn); 262 rn = n; 263 MPN_NORMALIZE (rp, rn); 264 } 265 266 ret: 267 MPZ_REALLOC (r, rn); 268 SIZ(r) = rn; 269 MPN_COPY (PTR(r), rp, rn); 270 271 TMP_FREE; 272 } 273