1/* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. 2 3Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2005 Free Software 4Foundation, Inc. 5 6This file is part of the GNU MP Library. 7 8The GNU MP Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MP Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 20 21 22#include "gmp.h" 23#include "gmp-impl.h" 24#include "longlong.h" 25 26/* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and 27 t is defined by (tp,mn). */ 28static void 29reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn) 30{ 31 mp_ptr qp; 32 TMP_DECL; 33 34 TMP_MARK; 35 qp = TMP_ALLOC_LIMBS (an - mn + 1); 36 37 mpn_tdiv_qr (qp, tp, 0L, ap, an, mp, mn); 38 39 TMP_FREE; 40} 41 42void 43mpz_powm_ui (mpz_ptr r, mpz_srcptr b, unsigned long int el, mpz_srcptr m) 44{ 45 mp_ptr xp, tp, qp, mp, bp; 46 mp_size_t xn, tn, mn, bn; 47 int m_zero_cnt; 48 int c; 49 mp_limb_t e; 50 TMP_DECL; 51 52 mp = PTR(m); 53 mn = ABSIZ(m); 54 if (mn == 0) 55 DIVIDE_BY_ZERO; 56 57 if (el == 0) 58 { 59 /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 60 depending on if MOD equals 1. */ 61 SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1; 62 PTR(r)[0] = 1; 63 return; 64 } 65 66 TMP_MARK; 67 68 /* Normalize m (i.e. make its most significant bit set) as required by 69 division functions below. */ 70 count_leading_zeros (m_zero_cnt, mp[mn - 1]); 71 m_zero_cnt -= GMP_NAIL_BITS; 72 if (m_zero_cnt != 0) 73 { 74 mp_ptr new_mp = TMP_ALLOC_LIMBS (mn); 75 mpn_lshift (new_mp, mp, mn, m_zero_cnt); 76 mp = new_mp; 77 } 78 79 bn = ABSIZ(b); 80 bp = PTR(b); 81 if (bn > mn) 82 { 83 /* Reduce possibly huge base. Use a function call to reduce, since we 84 don't want the quotient allocation to live until function return. */ 85 mp_ptr new_bp = TMP_ALLOC_LIMBS (mn); 86 reduce (new_bp, bp, bn, mp, mn); 87 bp = new_bp; 88 bn = mn; 89 /* Canonicalize the base, since we are potentially going to multiply with 90 it quite a few times. */ 91 MPN_NORMALIZE (bp, bn); 92 } 93 94 if (bn == 0) 95 { 96 SIZ(r) = 0; 97 TMP_FREE; 98 return; 99 } 100 101 tp = TMP_ALLOC_LIMBS (2 * mn + 1); 102 xp = TMP_ALLOC_LIMBS (mn); 103 104 qp = TMP_ALLOC_LIMBS (mn + 1); 105 106 MPN_COPY (xp, bp, bn); 107 xn = bn; 108 109 e = el; 110 count_leading_zeros (c, e); 111 e = (e << c) << 1; /* shift the exp bits to the left, lose msb */ 112 c = GMP_LIMB_BITS - 1 - c; 113 114 /* Main loop. */ 115 116 /* If m is already normalized (high bit of high limb set), and b is the 117 same size, but a bigger value, and e==1, then there's no modular 118 reductions done and we can end up with a result out of range at the 119 end. */ 120 if (c == 0) 121 { 122 if (xn == mn && mpn_cmp (xp, mp, mn) >= 0) 123 mpn_sub_n (xp, xp, mp, mn); 124 goto finishup; 125 } 126 127 while (c != 0) 128 { 129 mpn_sqr (tp, xp, xn); 130 tn = 2 * xn; tn -= tp[tn - 1] == 0; 131 if (tn < mn) 132 { 133 MPN_COPY (xp, tp, tn); 134 xn = tn; 135 } 136 else 137 { 138 mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn); 139 xn = mn; 140 } 141 142 if ((mp_limb_signed_t) e < 0) 143 { 144 mpn_mul (tp, xp, xn, bp, bn); 145 tn = xn + bn; tn -= tp[tn - 1] == 0; 146 if (tn < mn) 147 { 148 MPN_COPY (xp, tp, tn); 149 xn = tn; 150 } 151 else 152 { 153 mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn); 154 xn = mn; 155 } 156 } 157 e <<= 1; 158 c--; 159 } 160 161 finishup: 162 /* We shifted m left m_zero_cnt steps. Adjust the result by reducing 163 it with the original MOD. */ 164 if (m_zero_cnt != 0) 165 { 166 mp_limb_t cy; 167 cy = mpn_lshift (tp, xp, xn, m_zero_cnt); 168 tp[xn] = cy; xn += cy != 0; 169 170 if (xn < mn) 171 { 172 MPN_COPY (xp, tp, xn); 173 } 174 else 175 { 176 mpn_tdiv_qr (qp, xp, 0L, tp, xn, mp, mn); 177 xn = mn; 178 } 179 mpn_rshift (xp, xp, xn, m_zero_cnt); 180 } 181 MPN_NORMALIZE (xp, xn); 182 183 if ((el & 1) != 0 && SIZ(b) < 0 && xn != 0) 184 { 185 mp = PTR(m); /* want original, unnormalized m */ 186 mpn_sub (xp, mp, mn, xp, xn); 187 xn = mn; 188 MPN_NORMALIZE (xp, xn); 189 } 190 MPZ_REALLOC (r, xn); 191 SIZ (r) = xn; 192 MPN_COPY (PTR(r), xp, xn); 193 194 TMP_FREE; 195} 196