1 /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
3 Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation,
6 This file is part of the GNU MP Library.
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
13 The GNU MP Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21 MA 02111-1307, USA. */
23 #include <stdio.h> /* for NULL */
30 mpz_powm_ui (mpz_ptr res, mpz_srcptr base, unsigned long int exp, mpz_srcptr mod)
32 mpz_powm_ui (res, base, exp, mod)
35 unsigned long int exp;
40 mp_size_t msize, bsize, rsize;
44 mp_limb_t *free_me = NULL;
48 msize = ABS (mod->_mp_size);
58 /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
59 depending on if MOD equals 1. */
60 res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1;
67 /* Normalize MOD (i.e. make its most significant bit set) as required by
68 mpn_divmod. This will make the intermediate values in the calculation
69 slightly larger, but the correct result is obtained after a final
70 reduction using the original MOD value. */
72 mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB);
73 count_leading_zeros (mod_shift_cnt, mod->_mp_d[msize - 1]);
74 if (mod_shift_cnt != 0)
75 mpn_lshift (mp, mod->_mp_d, msize, mod_shift_cnt);
77 MPN_COPY (mp, mod->_mp_d, msize);
79 bsize = ABS (base->_mp_size);
82 /* The base is larger than the module. Reduce it. */
84 /* Allocate (BSIZE + 1) with space for remainder and quotient.
85 (The quotient is (bsize - msize + 1) limbs.) */
86 bp = (mp_ptr) TMP_ALLOC ((bsize + 1) * BYTES_PER_MP_LIMB);
87 MPN_COPY (bp, base->_mp_d, bsize);
88 /* We don't care about the quotient, store it above the remainder,
90 mpn_divmod (bp + msize, bp, bsize, mp, msize);
92 /* Canonicalize the base, since we are going to multiply with it
94 MPN_NORMALIZE (bp, bsize);
106 if (res->_mp_alloc < size)
108 /* We have to allocate more space for RES. If any of the input
109 parameters are identical to RES, defer deallocation of the old
112 if (rp == mp || rp == bp)
115 free_me_size = res->_mp_alloc;
118 (*_mp_free_func) (rp, res->_mp_alloc * BYTES_PER_MP_LIMB);
120 rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB);
121 res->_mp_alloc = size;
126 /* Make BASE, EXP and MOD not overlap with RES. */
129 /* RES and BASE are identical. Allocate temp. space for BASE. */
130 bp = (mp_ptr) TMP_ALLOC (bsize * BYTES_PER_MP_LIMB);
131 MPN_COPY (bp, rp, bsize);
135 /* RES and MOD are identical. Allocate temporary space for MOD. */
136 mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB);
137 MPN_COPY (mp, rp, msize);
141 MPN_COPY (rp, bp, bsize);
145 mp_ptr xp = (mp_ptr) TMP_ALLOC (2 * (msize + 1) * BYTES_PER_MP_LIMB);
148 mp_limb_t carry_limb;
150 negative_result = (exp & 1) && base->_mp_size < 0;
153 count_leading_zeros (c, e);
154 e = (e << c) << 1; /* shift the exp bits to the left, lose msb */
155 c = BITS_PER_MP_LIMB - 1 - c;
159 Make the result be pointed to alternately by XP and RP. This
160 helps us avoid block copying, which would otherwise be necessary
161 with the overlap restrictions of mpn_divmod. With 50% probability
162 the result after this loop will be in the area originally pointed
163 by RP (==RES->_mp_d), and with 50% probability in the area originally
171 mpn_mul_n (xp, rp, rp, rsize);
173 xsize -= xp[xsize - 1] == 0;
176 mpn_divmod (xp + msize, xp, xsize, mp, msize);
180 tp = rp; rp = xp; xp = tp;
183 if ((mp_limb_signed_t) e < 0)
185 mpn_mul (xp, rp, rsize, bp, bsize);
186 xsize = rsize + bsize;
187 xsize -= xp[xsize - 1] == 0;
190 mpn_divmod (xp + msize, xp, xsize, mp, msize);
194 tp = rp; rp = xp; xp = tp;
201 /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
202 steps. Adjust the result by reducing it with the original MOD.
204 Also make sure the result is put in RES->_mp_d (where it already
205 might be, see above). */
207 if (mod_shift_cnt != 0)
209 carry_limb = mpn_lshift (res->_mp_d, rp, rsize, mod_shift_cnt);
213 rp[rsize] = carry_limb;
219 MPN_COPY (res->_mp_d, rp, rsize);
225 mpn_divmod (rp + msize, rp, rsize, mp, msize);
229 /* Remove any leading zero words from the result. */
230 if (mod_shift_cnt != 0)
231 mpn_rshift (rp, rp, rsize, mod_shift_cnt);
232 MPN_NORMALIZE (rp, rsize);
235 if (negative_result && rsize != 0)
237 if (mod_shift_cnt != 0)
238 mpn_rshift (mp, mp, msize, mod_shift_cnt);
239 mpn_sub (rp, mp, msize, rp, rsize);
241 MPN_NORMALIZE (rp, rsize);
243 res->_mp_size = rsize;
246 (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);