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, 2009 6 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 #ifdef BERKELEY_MP 28 #include "mp.h" 29 #endif 30 31 32 /* TODO 33 34 * Improve handling of buffers. It is pretty ugly now. 35 36 * For even moduli, we compute a binvert of its odd part both here and in 37 mpn_powm. How can we avoid this recomputation? 38 */ 39 40 /* 41 b ^ e mod m res 42 0 0 0 ? 43 0 e 0 ? 44 0 0 m ? 45 0 e m 0 46 b 0 0 ? 47 b e 0 ? 48 b 0 m 1 mod m 49 b e m b^e mod m 50 */ 51 52 #define HANDLE_NEGATIVE_EXPONENT 1 53 54 void 55 #ifndef BERKELEY_MP 56 mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m) 57 #else /* BERKELEY_MP */ 58 pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) 59 #endif /* BERKELEY_MP */ 60 { 61 mp_size_t n, nodd, ncnt; 62 int cnt; 63 mp_ptr rp, tp; 64 mp_srcptr bp, ep, mp; 65 mp_size_t rn, bn, es, en, itch; 66 TMP_DECL; 67 68 n = ABSIZ(m); 69 if (n == 0) 70 DIVIDE_BY_ZERO; 71 72 mp = PTR(m); 73 74 TMP_MARK; 75 76 es = SIZ(e); 77 if (UNLIKELY (es <= 0)) 78 { 79 mpz_t new_b; 80 if (es == 0) 81 { 82 /* b^0 mod m, b is anything and m is non-zero. 83 Result is 1 mod m, i.e., 1 or 0 depending on if m = 1. */ 84 SIZ(r) = n != 1 || mp[0] != 1; 85 PTR(r)[0] = 1; 86 TMP_FREE; /* we haven't really allocated anything here */ 87 return; 88 } 89 #if HANDLE_NEGATIVE_EXPONENT 90 MPZ_TMP_INIT (new_b, n + 1); 91 92 if (! mpz_invert (new_b, b, m)) 93 DIVIDE_BY_ZERO; 94 b = new_b; 95 es = -es; 96 #else 97 DIVIDE_BY_ZERO; 98 #endif 99 } 100 en = es; 101 102 bn = ABSIZ(b); 103 104 if (UNLIKELY (bn == 0)) 105 { 106 SIZ(r) = 0; 107 TMP_FREE; 108 return; 109 } 110 111 ep = PTR(e); 112 113 /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case. */ 114 if (UNLIKELY (en == 1 && ep[0] == 1)) 115 { 116 rp = TMP_ALLOC_LIMBS (n); 117 bp = PTR(b); 118 if (bn >= n) 119 { 120 mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1); 121 mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n); 122 rn = n; 123 MPN_NORMALIZE (rp, rn); 124 125 if (SIZ(b) < 0 && rn != 0) 126 { 127 mpn_sub (rp, mp, n, rp, rn); 128 rn = n; 129 MPN_NORMALIZE (rp, rn); 130 } 131 } 132 else 133 { 134 if (SIZ(b) < 0) 135 { 136 mpn_sub (rp, mp, n, bp, bn); 137 rn = n; 138 rn -= (rp[rn - 1] == 0); 139 } 140 else 141 { 142 MPN_COPY (rp, bp, bn); 143 rn = bn; 144 } 145 } 146 goto ret; 147 } 148 149 /* Remove low zero limbs from M. This loop will terminate for correctly 150 represented mpz numbers. */ 151 ncnt = 0; 152 while (UNLIKELY (mp[0] == 0)) 153 { 154 mp++; 155 ncnt++; 156 } 157 nodd = n - ncnt; 158 cnt = 0; 159 if (mp[0] % 2 == 0) 160 { 161 mp_ptr new = TMP_ALLOC_LIMBS (nodd); 162 count_trailing_zeros (cnt, mp[0]); 163 mpn_rshift (new, mp, nodd, cnt); 164 nodd -= new[nodd - 1] == 0; 165 mp = new; 166 ncnt++; 167 } 168 169 if (ncnt != 0) 170 { 171 /* We will call both mpn_powm and mpn_powlo. */ 172 /* rp needs n, mpn_powlo needs 4n, the 2 mpn_binvert might need more */ 173 mp_size_t n_largest_binvert = MAX (ncnt, nodd); 174 mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert); 175 itch = 3 * n + MAX (itch_binvert, 2 * n); 176 } 177 else 178 { 179 /* We will call just mpn_powm. */ 180 mp_size_t itch_binvert = mpn_binvert_itch (nodd); 181 itch = n + MAX (itch_binvert, 2 * n); 182 } 183 tp = TMP_ALLOC_LIMBS (itch); 184 185 rp = tp; tp += n; 186 187 bp = PTR(b); 188 mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp); 189 190 rn = n; 191 192 if (ncnt != 0) 193 { 194 mp_ptr r2, xp, yp, odd_inv_2exp; 195 unsigned long t; 196 int bcnt; 197 198 if (bn < ncnt) 199 { 200 mp_ptr new = TMP_ALLOC_LIMBS (ncnt); 201 MPN_COPY (new, bp, bn); 202 MPN_ZERO (new + bn, ncnt - bn); 203 bp = new; 204 } 205 206 r2 = tp; 207 208 if (bp[0] % 2 == 0) 209 { 210 if (en > 1) 211 { 212 MPN_ZERO (r2, ncnt); 213 goto zero; 214 } 215 216 ASSERT (en == 1); 217 t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt; 218 219 /* Count number of low zero bits in B, up to 3. */ 220 bcnt = (0x1213 >> ((bp[0] & 7) << 1)) & 0x3; 221 /* Note that ep[0] * bcnt might overflow, but that just results 222 in a missed optimization. */ 223 if (ep[0] * bcnt >= t) 224 { 225 MPN_ZERO (r2, ncnt); 226 goto zero; 227 } 228 } 229 230 mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt); 231 232 zero: 233 if (nodd < ncnt) 234 { 235 mp_ptr new = TMP_ALLOC_LIMBS (ncnt); 236 MPN_COPY (new, mp, nodd); 237 MPN_ZERO (new + nodd, ncnt - nodd); 238 mp = new; 239 } 240 241 odd_inv_2exp = tp + n; 242 mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n); 243 244 mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd); 245 246 xp = tp + 2 * n; 247 mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt); 248 249 if (cnt != 0) 250 xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1; 251 252 yp = tp; 253 if (ncnt > nodd) 254 mpn_mul (yp, xp, ncnt, mp, nodd); 255 else 256 mpn_mul (yp, mp, nodd, xp, ncnt); 257 258 mpn_add (rp, yp, n, rp, nodd); 259 260 ASSERT (nodd + ncnt >= n); 261 ASSERT (nodd + ncnt <= n + 1); 262 } 263 264 MPN_NORMALIZE (rp, rn); 265 266 if ((ep[0] & 1) && SIZ(b) < 0 && rn != 0) 267 { 268 mpn_sub (rp, PTR(m), n, rp, rn); 269 rn = n; 270 MPN_NORMALIZE (rp, rn); 271 } 272 273 ret: 274 MPZ_REALLOC (r, rn); 275 SIZ(r) = rn; 276 MPN_COPY (PTR(r), rp, rn); 277 278 TMP_FREE; 279 } 280