1 /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
3 Copyright (C) 1991, 1992 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 General Public License as published by
9 the Free Software Foundation; either version 2, or (at your option)
12 The GNU MP Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
17 You should have received a copy of the GNU General Public License
18 along with the GNU MP Library; see the file COPYING. If not, write to
19 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
27 mpz_powm_ui (MP_INT *res, const MP_INT *base, unsigned long int exp,
30 mpz_powm_ui (res, base, exp, mod)
33 unsigned long int exp;
38 mp_size msize, bsize, rsize;
42 mp_limb *free_me = NULL;
45 msize = ABS (mod->size);
50 /* Normalize MOD (i.e. make its most significant bit set) as required by
51 mpn_div. This will make the intermediate values in the calculation
52 slightly larger, but the correct result is obtained after a final
53 reduction using the original MOD value. */
55 mp = (mp_ptr) alloca (msize * BYTES_PER_MP_LIMB);
56 count_leading_zeros (mod_shift_cnt, mod->d[msize - 1]);
57 if (mod_shift_cnt != 0)
58 (void) mpn_lshift (mp, mod->d, msize, mod_shift_cnt);
60 MPN_COPY (mp, mod->d, msize);
62 bsize = ABS (base->size);
65 /* The base is larger than the module. Reduce it. */
67 /* Allocate (BSIZE + 1) with space for remainder and quotient.
68 (The quotient is (bsize - msize + 1) limbs.) */
69 bp = (mp_ptr) alloca ((bsize + 1) * BYTES_PER_MP_LIMB);
70 MPN_COPY (bp, base->d, bsize);
71 /* We don't care about the quotient, store it above the remainder,
73 mpn_div (bp + msize, bp, bsize, mp, msize);
75 while (bsize > 0 && bp[bsize - 1] == 0)
81 bsize = ABS (base->size);
84 if (res->alloc < size)
86 /* We have to allocate more space for RES. If any of the input
87 parameters are identical to RES, defer deallocation of the old
90 if (rp == mp || rp == bp)
93 free_me_size = res->alloc;
96 (*_mp_free_func) (rp, res->alloc * BYTES_PER_MP_LIMB);
98 rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB);
104 /* Make BASE, EXP and MOD not overlap with RES. */
107 /* RES and BASE are identical. Allocate temp. space for BASE. */
108 bp = (mp_ptr) alloca (bsize * BYTES_PER_MP_LIMB);
109 MPN_COPY (bp, rp, bsize);
113 /* RES and MOD are identical. Allocate temporary space for MOD. */
114 mp = (mp_ptr) alloca (msize * BYTES_PER_MP_LIMB);
115 MPN_COPY (mp, rp, msize);
126 MPN_COPY (rp, bp, bsize);
131 mp_ptr dummyp = (mp_ptr) alloca ((msize + 1) * BYTES_PER_MP_LIMB);
132 mp_ptr xp = (mp_ptr) alloca (2 * (msize + 1) * BYTES_PER_MP_LIMB);
137 negative_result = (exp & 1) && base->size < 0;
140 count_leading_zeros (c, e);
141 e <<= (c + 1); /* shift the exp bits to the left, loose msb */
142 c = BITS_PER_MP_LIMB - 1 - c;
146 Make the result be pointed to alternately by XP and RP. This
147 helps us avoid block copying, which would otherwise be necessary
148 with the overlap restrictions of mpn_div. With 50% probability
149 the result after this loop will be in the area originally pointed
150 by RP (==RES->D), and with 50% probability in the area originally
158 xsize = mpn_mul (xp, rp, rsize, rp, rsize);
159 mpn_div (dummyp, xp, xsize, mp, msize);
161 /* Remove any leading zero words from the result. */
164 while (xsize > 0 && xp[xsize - 1] == 0)
167 tp = rp; rp = xp; xp = tp;
168 tsize = rsize; rsize = xsize; xsize = tsize;
170 if ((mp_limb_signed) e < 0)
173 xsize = mpn_mul (xp, rp, rsize, bp, bsize);
175 xsize = mpn_mul (xp, bp, bsize, rp, rsize);
176 mpn_div (dummyp, xp, xsize, mp, msize);
178 /* Remove any leading zero words from the result. */
181 while (xsize > 0 && xp[xsize - 1] == 0)
184 tp = rp; rp = xp; xp = tp;
185 tsize = rsize; rsize = xsize; xsize = tsize;
191 /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
192 steps. Adjust the result by reducing it with the original MOD.
194 Also make sure the result is put in RES->D (where it already
195 might be, see above). */
197 carry_limb = mpn_lshift (res->d, rp, rsize, mod_shift_cnt);
201 rp[rsize] = carry_limb;
204 mpn_div (dummyp, rp, rsize, mp, msize);
205 /* Remove any leading zero words from the result. */
208 while (rsize > 0 && rp[rsize - 1] == 0)
210 rsize = mpn_rshift (rp, rp, rsize, mod_shift_cnt);
213 res->size = negative_result >= 0 ? rsize : -rsize;
216 (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);