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