=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpz/Attic/powm_ui.c,v retrieving revision 1.1.1.2 retrieving revision 1.1.1.3 diff -u -p -r1.1.1.2 -r1.1.1.3 --- OpenXM_contrib/gmp/mpz/Attic/powm_ui.c 2000/09/09 14:12:56 1.1.1.2 +++ OpenXM_contrib/gmp/mpz/Attic/powm_ui.c 2003/08/25 16:06:33 1.1.1.3 @@ -1,7 +1,7 @@ /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. -Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation, -Inc. +Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002 Free Software +Foundation, Inc. This file is part of the GNU MP Library. @@ -20,229 +20,178 @@ along with the GNU MP Library; see the file COPYING.LI the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include /* for NULL */ + #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" -void -#if __STDC__ -mpz_powm_ui (mpz_ptr res, mpz_srcptr base, unsigned long int exp, mpz_srcptr mod) -#else -mpz_powm_ui (res, base, exp, mod) - mpz_ptr res; - mpz_srcptr base; - unsigned long int exp; - mpz_srcptr mod; -#endif +/* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and + t is defined by (tp,mn). */ +static void +reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn) { - mp_ptr rp, mp, bp; - mp_size_t msize, bsize, rsize; - mp_size_t size; - int mod_shift_cnt; - int negative_result; - mp_limb_t *free_me = NULL; - size_t free_me_size; + mp_ptr qp; TMP_DECL (marker); - msize = ABS (mod->_mp_size); - size = 2 * msize; + TMP_MARK (marker); + qp = TMP_ALLOC_LIMBS (an - mn + 1); - rp = res->_mp_d; + mpn_tdiv_qr (qp, tp, 0L, ap, an, mp, mn); - if (msize == 0) + TMP_FREE (marker); +} + +void +mpz_powm_ui (mpz_ptr r, mpz_srcptr b, unsigned long int el, mpz_srcptr m) +{ + mp_ptr xp, tp, qp, mp, bp; + mp_size_t xn, tn, mn, bn; + int m_zero_cnt; + int c; + mp_limb_t e; + TMP_DECL (marker); + + mp = PTR(m); + mn = ABSIZ(m); + if (mn == 0) DIVIDE_BY_ZERO; - if (exp == 0) + if (el == 0) { /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 depending on if MOD equals 1. */ - res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1; - rp[0] = 1; + SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1; + PTR(r)[0] = 1; return; } TMP_MARK (marker); - /* Normalize MOD (i.e. make its most significant bit set) as required by - mpn_divmod. This will make the intermediate values in the calculation - slightly larger, but the correct result is obtained after a final - reduction using the original MOD value. */ - - mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB); - count_leading_zeros (mod_shift_cnt, mod->_mp_d[msize - 1]); - if (mod_shift_cnt != 0) - mpn_lshift (mp, mod->_mp_d, msize, mod_shift_cnt); - else - MPN_COPY (mp, mod->_mp_d, msize); - - bsize = ABS (base->_mp_size); - if (bsize > msize) + /* Normalize m (i.e. make its most significant bit set) as required by + division functions below. */ + count_leading_zeros (m_zero_cnt, mp[mn - 1]); + m_zero_cnt -= GMP_NAIL_BITS; + if (m_zero_cnt != 0) { - /* The base is larger than the module. Reduce it. */ + mp_ptr new_mp = TMP_ALLOC_LIMBS (mn); + mpn_lshift (new_mp, mp, mn, m_zero_cnt); + mp = new_mp; + } - /* Allocate (BSIZE + 1) with space for remainder and quotient. - (The quotient is (bsize - msize + 1) limbs.) */ - bp = (mp_ptr) TMP_ALLOC ((bsize + 1) * BYTES_PER_MP_LIMB); - MPN_COPY (bp, base->_mp_d, bsize); - /* We don't care about the quotient, store it above the remainder, - at BP + MSIZE. */ - mpn_divmod (bp + msize, bp, bsize, mp, msize); - bsize = msize; - /* Canonicalize the base, since we are going to multiply with it - quite a few times. */ - MPN_NORMALIZE (bp, bsize); + bn = ABSIZ(b); + bp = PTR(b); + if (bn > mn) + { + /* Reduce possibly huge base. Use a function call to reduce, since we + don't want the quotient allocation to live until function return. */ + mp_ptr new_bp = TMP_ALLOC_LIMBS (mn); + reduce (new_bp, bp, bn, mp, mn); + bp = new_bp; + bn = mn; + /* Canonicalize the base, since we are potentially going to multiply with + it quite a few times. */ + MPN_NORMALIZE (bp, bn); } - else - bp = base->_mp_d; - if (bsize == 0) + if (bn == 0) { - res->_mp_size = 0; + SIZ(r) = 0; TMP_FREE (marker); return; } - if (res->_mp_alloc < size) + tp = TMP_ALLOC_LIMBS (2 * mn + 1); + xp = TMP_ALLOC_LIMBS (mn); + + qp = TMP_ALLOC_LIMBS (mn + 1); + + MPN_COPY (xp, bp, bn); + xn = bn; + + e = el; + count_leading_zeros (c, e); + e = (e << c) << 1; /* shift the exp bits to the left, lose msb */ + c = BITS_PER_MP_LIMB - 1 - c; + + /* Main loop. */ + + /* If m is already normalized (high bit of high limb set), and b is the + same size, but a bigger value, and e==1, then there's no modular + reductions done and we can end up with a result out of range at the + end. */ + if (c == 0) { - /* We have to allocate more space for RES. If any of the input - parameters are identical to RES, defer deallocation of the old - space. */ + if (xn == mn && mpn_cmp (xp, mp, mn) >= 0) + mpn_sub_n (xp, xp, mp, mn); + goto finishup; + } - if (rp == mp || rp == bp) + while (c != 0) + { + mpn_sqr_n (tp, xp, xn); + tn = 2 * xn; tn -= tp[tn - 1] == 0; + if (tn < mn) { - free_me = rp; - free_me_size = res->_mp_alloc; + MPN_COPY (xp, tp, tn); + xn = tn; } else - (*_mp_free_func) (rp, res->_mp_alloc * BYTES_PER_MP_LIMB); + { + mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn); + xn = mn; + } - rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB); - res->_mp_alloc = size; - res->_mp_d = rp; + if ((mp_limb_signed_t) e < 0) + { + mpn_mul (tp, xp, xn, bp, bn); + tn = xn + bn; tn -= tp[tn - 1] == 0; + if (tn < mn) + { + MPN_COPY (xp, tp, tn); + xn = tn; + } + else + { + mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn); + xn = mn; + } + } + e <<= 1; + c--; } - else + + finishup: + /* We shifted m left m_zero_cnt steps. Adjust the result by reducing + it with the original MOD. */ + if (m_zero_cnt != 0) { - /* Make BASE, EXP and MOD not overlap with RES. */ - if (rp == bp) + mp_limb_t cy; + cy = mpn_lshift (tp, xp, xn, m_zero_cnt); + tp[xn] = cy; xn += cy != 0; + + if (xn < mn) { - /* RES and BASE are identical. Allocate temp. space for BASE. */ - bp = (mp_ptr) TMP_ALLOC (bsize * BYTES_PER_MP_LIMB); - MPN_COPY (bp, rp, bsize); + MPN_COPY (xp, tp, xn); } - if (rp == mp) + else { - /* RES and MOD are identical. Allocate temporary space for MOD. */ - mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB); - MPN_COPY (mp, rp, msize); + mpn_tdiv_qr (qp, xp, 0L, tp, xn, mp, mn); + xn = mn; } + mpn_rshift (xp, xp, xn, m_zero_cnt); } + MPN_NORMALIZE (xp, xn); - MPN_COPY (rp, bp, bsize); - rsize = bsize; - - { - mp_ptr xp = (mp_ptr) TMP_ALLOC (2 * (msize + 1) * BYTES_PER_MP_LIMB); - int c; - mp_limb_t e; - mp_limb_t carry_limb; - - negative_result = (exp & 1) && base->_mp_size < 0; - - e = exp; - count_leading_zeros (c, e); - e = (e << c) << 1; /* shift the exp bits to the left, lose msb */ - c = BITS_PER_MP_LIMB - 1 - c; - - /* Main loop. - - Make the result be pointed to alternately by XP and RP. This - helps us avoid block copying, which would otherwise be necessary - with the overlap restrictions of mpn_divmod. With 50% probability - the result after this loop will be in the area originally pointed - by RP (==RES->_mp_d), and with 50% probability in the area originally - pointed to by XP. */ - - while (c != 0) - { - mp_ptr tp; - mp_size_t xsize; - - mpn_mul_n (xp, rp, rp, rsize); - xsize = 2 * rsize; - xsize -= xp[xsize - 1] == 0; - if (xsize > msize) - { - mpn_divmod (xp + msize, xp, xsize, mp, msize); - xsize = msize; - } - - tp = rp; rp = xp; xp = tp; - rsize = xsize; - - if ((mp_limb_signed_t) e < 0) - { - mpn_mul (xp, rp, rsize, bp, bsize); - xsize = rsize + bsize; - xsize -= xp[xsize - 1] == 0; - if (xsize > msize) - { - mpn_divmod (xp + msize, xp, xsize, mp, msize); - xsize = msize; - } - - tp = rp; rp = xp; xp = tp; - rsize = xsize; - } - e <<= 1; - c--; - } - - /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT - steps. Adjust the result by reducing it with the original MOD. - - Also make sure the result is put in RES->_mp_d (where it already - might be, see above). */ - - if (mod_shift_cnt != 0) - { - carry_limb = mpn_lshift (res->_mp_d, rp, rsize, mod_shift_cnt); - rp = res->_mp_d; - if (carry_limb != 0) - { - rp[rsize] = carry_limb; - rsize++; - } - } - else - { - MPN_COPY (res->_mp_d, rp, rsize); - rp = res->_mp_d; - } - - if (rsize >= msize) - { - mpn_divmod (rp + msize, rp, rsize, mp, msize); - rsize = msize; - } - - /* Remove any leading zero words from the result. */ - if (mod_shift_cnt != 0) - mpn_rshift (rp, rp, rsize, mod_shift_cnt); - MPN_NORMALIZE (rp, rsize); - } - - if (negative_result && rsize != 0) + if ((el & 1) != 0 && SIZ(b) < 0 && xn != 0) { - if (mod_shift_cnt != 0) - mpn_rshift (mp, mp, msize, mod_shift_cnt); - mpn_sub (rp, mp, msize, rp, rsize); - rsize = msize; - MPN_NORMALIZE (rp, rsize); + mp = PTR(m); /* want original, unnormalized m */ + mpn_sub (xp, mp, mn, xp, xn); + xn = mn; + MPN_NORMALIZE (xp, xn); } - res->_mp_size = rsize; + MPZ_REALLOC (r, xn); + SIZ (r) = xn; + MPN_COPY (PTR(r), xp, xn); - if (free_me != NULL) - (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB); TMP_FREE (marker); }