1 /* mpz_powm(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
3 Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Library General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or (at your
10 option) any later version.
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15 License for more details.
17 You should have received a copy of the GNU Library General Public License
18 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20 MA 02111-1307, USA. */
29 mpz_powm (mpz_ptr res, mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod)
31 mpz_powm (res, base, exp, mod)
37 #else /* BERKELEY_MP */
40 pow (mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod, mpz_ptr res)
42 pow (base, exp, mod, res)
48 #endif /* BERKELEY_MP */
50 mp_ptr rp, ep, mp, bp;
51 mp_size_t esize, msize, bsize, rsize;
55 mp_limb_t *free_me = NULL;
59 esize = ABS (exp->_mp_size);
60 msize = ABS (mod->_mp_size);
67 msize = 1 / msize; /* provoke a signal */
71 /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
72 depending on if MOD equals 1. */
74 res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1;
80 /* Normalize MOD (i.e. make its most significant bit set) as required by
81 mpn_divmod. This will make the intermediate values in the calculation
82 slightly larger, but the correct result is obtained after a final
83 reduction using the original MOD value. */
85 mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB);
86 count_leading_zeros (mod_shift_cnt, mod->_mp_d[msize - 1]);
87 if (mod_shift_cnt != 0)
88 mpn_lshift (mp, mod->_mp_d, msize, mod_shift_cnt);
90 MPN_COPY (mp, mod->_mp_d, msize);
92 bsize = ABS (base->_mp_size);
95 /* The base is larger than the module. Reduce it. */
97 /* Allocate (BSIZE + 1) with space for remainder and quotient.
98 (The quotient is (bsize - msize + 1) limbs.) */
99 bp = (mp_ptr) TMP_ALLOC ((bsize + 1) * BYTES_PER_MP_LIMB);
100 MPN_COPY (bp, base->_mp_d, bsize);
101 /* We don't care about the quotient, store it above the remainder,
103 mpn_divmod (bp + msize, bp, bsize, mp, msize);
105 /* Canonicalize the base, since we are going to multiply with it
106 quite a few times. */
107 MPN_NORMALIZE (bp, bsize);
119 if (res->_mp_alloc < size)
121 /* We have to allocate more space for RES. If any of the input
122 parameters are identical to RES, defer deallocation of the old
125 if (rp == ep || rp == mp || rp == bp)
128 free_me_size = res->_mp_alloc;
131 (*_mp_free_func) (rp, res->_mp_alloc * BYTES_PER_MP_LIMB);
133 rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB);
134 res->_mp_alloc = size;
139 /* Make BASE, EXP and MOD not overlap with RES. */
142 /* RES and BASE are identical. Allocate temp. space for BASE. */
143 bp = (mp_ptr) TMP_ALLOC (bsize * BYTES_PER_MP_LIMB);
144 MPN_COPY (bp, rp, bsize);
148 /* RES and EXP are identical. Allocate temp. space for EXP. */
149 ep = (mp_ptr) TMP_ALLOC (esize * BYTES_PER_MP_LIMB);
150 MPN_COPY (ep, rp, esize);
154 /* RES and MOD are identical. Allocate temporary space for MOD. */
155 mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB);
156 MPN_COPY (mp, rp, msize);
160 MPN_COPY (rp, bp, bsize);
165 mp_ptr xp = (mp_ptr) TMP_ALLOC (2 * (msize + 1) * BYTES_PER_MP_LIMB);
168 mp_limb_t carry_limb;
170 negative_result = (ep[0] & 1) && base->_mp_size < 0;
174 count_leading_zeros (c, e);
175 e = (e << c) << 1; /* shift the exp bits to the left, lose msb */
176 c = BITS_PER_MP_LIMB - 1 - c;
180 Make the result be pointed to alternately by XP and RP. This
181 helps us avoid block copying, which would otherwise be necessary
182 with the overlap restrictions of mpn_divmod. With 50% probability
183 the result after this loop will be in the area originally pointed
184 by RP (==RES->_mp_d), and with 50% probability in the area originally
194 mpn_mul_n (xp, rp, rp, rsize);
198 mpn_divmod (xp + msize, xp, xsize, mp, msize);
202 tp = rp; rp = xp; xp = tp;
205 if ((mp_limb_signed_t) e < 0)
207 mpn_mul (xp, rp, rsize, bp, bsize);
208 xsize = rsize + bsize;
211 mpn_divmod (xp + msize, xp, xsize, mp, msize);
215 tp = rp; rp = xp; xp = tp;
226 c = BITS_PER_MP_LIMB;
229 /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
230 steps. Adjust the result by reducing it with the original MOD.
232 Also make sure the result is put in RES->_mp_d (where it already
233 might be, see above). */
235 if (mod_shift_cnt != 0)
237 carry_limb = mpn_lshift (res->_mp_d, rp, rsize, mod_shift_cnt);
241 rp[rsize] = carry_limb;
247 MPN_COPY (res->_mp_d, rp, rsize);
253 mpn_divmod (rp + msize, rp, rsize, mp, msize);
257 /* Remove any leading zero words from the result. */
258 if (mod_shift_cnt != 0)
259 mpn_rshift (rp, rp, rsize, mod_shift_cnt);
260 MPN_NORMALIZE (rp, rsize);
263 if (negative_result && rsize != 0)
265 if (mod_shift_cnt != 0)
266 mpn_rshift (mp, mp, msize, mod_shift_cnt);
267 mpn_sub (rp, mp, msize, rp, rsize);
269 MPN_NORMALIZE (rp, rsize);
271 res->_mp_size = rsize;
274 (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);