[BACK]Return to powm.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2-ssh-2 / mpz

Annotation of OpenXM/src/kan96xx/gmp-2.0.2-ssh-2/mpz/powm.c, Revision 1.1.1.1

1.1       takayama    1: /* mpz_powm(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
                      2:
                      3: Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
                      4:
                      5: This file is part of the GNU MP Library.
                      6:
                      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.
                     11:
                     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.
                     16:
                     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. */
                     21:
                     22: #include "gmp.h"
                     23: #include "gmp-impl.h"
                     24: #include "longlong.h"
                     25:
                     26: #ifndef BERKELEY_MP
                     27: void
                     28: #if __STDC__
                     29: mpz_powm (mpz_ptr res, mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod)
                     30: #else
                     31: mpz_powm (res, base, exp, mod)
                     32:      mpz_ptr res;
                     33:      mpz_srcptr base;
                     34:      mpz_srcptr exp;
                     35:      mpz_srcptr mod;
                     36: #endif
                     37: #else /* BERKELEY_MP */
                     38: void
                     39: #if __STDC__
                     40: pow (mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod, mpz_ptr res)
                     41: #else
                     42: pow (base, exp, mod, res)
                     43:      mpz_srcptr base;
                     44:      mpz_srcptr exp;
                     45:      mpz_srcptr mod;
                     46:      mpz_ptr res;
                     47: #endif
                     48: #endif /* BERKELEY_MP */
                     49: {
                     50:   mp_ptr rp, ep, mp, bp;
                     51:   mp_size_t esize, msize, bsize, rsize;
                     52:   mp_size_t size;
                     53:   int mod_shift_cnt;
                     54:   int negative_result;
                     55:   mp_limb_t *free_me = NULL;
                     56:   size_t free_me_size;
                     57:   TMP_DECL (marker);
                     58:
                     59:   esize = ABS (exp->_mp_size);
                     60:   msize = ABS (mod->_mp_size);
                     61:   size = 2 * msize;
                     62:
                     63:   rp = res->_mp_d;
                     64:   ep = exp->_mp_d;
                     65:
                     66:   if (msize == 0)
                     67:     msize = 1 / msize;         /* provoke a signal */
                     68:
                     69:   if (esize == 0)
                     70:     {
                     71:       /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
                     72:         depending on if MOD equals 1.  */
                     73:       rp[0] = 1;
                     74:       res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1;
                     75:       return;
                     76:     }
                     77:
                     78:   TMP_MARK (marker);
                     79:
                     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.  */
                     84:
                     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);
                     89:   else
                     90:     MPN_COPY (mp, mod->_mp_d, msize);
                     91:
                     92:   bsize = ABS (base->_mp_size);
                     93:   if (bsize > msize)
                     94:     {
                     95:       /* The base is larger than the module.  Reduce it.  */
                     96:
                     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,
                    102:         at BP + MSIZE.  */
                    103:       mpn_divmod (bp + msize, bp, bsize, mp, msize);
                    104:       bsize = msize;
                    105:       /* Canonicalize the base, since we are going to multiply with it
                    106:         quite a few times.  */
                    107:       MPN_NORMALIZE (bp, bsize);
                    108:     }
                    109:   else
                    110:     bp = base->_mp_d;
                    111:
                    112:   if (bsize == 0)
                    113:     {
                    114:       res->_mp_size = 0;
                    115:       TMP_FREE (marker);
                    116:       return;
                    117:     }
                    118:
                    119:   if (res->_mp_alloc < size)
                    120:     {
                    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
                    123:         space.  */
                    124:
                    125:       if (rp == ep || rp == mp || rp == bp)
                    126:        {
                    127:          free_me = rp;
                    128:          free_me_size = res->_mp_alloc;
                    129:        }
                    130:       else
                    131:        (*_mp_free_func) (rp, res->_mp_alloc * BYTES_PER_MP_LIMB);
                    132:
                    133:       rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB);
                    134:       res->_mp_alloc = size;
                    135:       res->_mp_d = rp;
                    136:     }
                    137:   else
                    138:     {
                    139:       /* Make BASE, EXP and MOD not overlap with RES.  */
                    140:       if (rp == bp)
                    141:        {
                    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);
                    145:        }
                    146:       if (rp == ep)
                    147:        {
                    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);
                    151:        }
                    152:       if (rp == mp)
                    153:        {
                    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);
                    157:        }
                    158:     }
                    159:
                    160:   MPN_COPY (rp, bp, bsize);
                    161:   rsize = bsize;
                    162:
                    163:   {
                    164:     mp_size_t i;
                    165:     mp_ptr xp = (mp_ptr) TMP_ALLOC (2 * (msize + 1) * BYTES_PER_MP_LIMB);
                    166:     int c;
                    167:     mp_limb_t e;
                    168:     mp_limb_t carry_limb;
                    169:
                    170:     negative_result = (ep[0] & 1) && base->_mp_size < 0;
                    171:
                    172:     i = esize - 1;
                    173:     e = ep[i];
                    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;
                    177:
                    178:     /* Main loop.
                    179:
                    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
                    185:        pointed to by XP.  */
                    186:
                    187:     for (;;)
                    188:       {
                    189:        while (c != 0)
                    190:          {
                    191:            mp_ptr tp;
                    192:            mp_size_t xsize;
                    193:
                    194:            mpn_mul_n (xp, rp, rp, rsize);
                    195:            xsize = 2 * rsize;
                    196:            if (xsize > msize)
                    197:              {
                    198:                mpn_divmod (xp + msize, xp, xsize, mp, msize);
                    199:                xsize = msize;
                    200:              }
                    201:
                    202:            tp = rp; rp = xp; xp = tp;
                    203:            rsize = xsize;
                    204:
                    205:            if ((mp_limb_signed_t) e < 0)
                    206:              {
                    207:                mpn_mul (xp, rp, rsize, bp, bsize);
                    208:                xsize = rsize + bsize;
                    209:                if (xsize > msize)
                    210:                  {
                    211:                    mpn_divmod (xp + msize, xp, xsize, mp, msize);
                    212:                    xsize = msize;
                    213:                  }
                    214:
                    215:                tp = rp; rp = xp; xp = tp;
                    216:                rsize = xsize;
                    217:              }
                    218:            e <<= 1;
                    219:            c--;
                    220:          }
                    221:
                    222:        i--;
                    223:        if (i < 0)
                    224:          break;
                    225:        e = ep[i];
                    226:        c = BITS_PER_MP_LIMB;
                    227:       }
                    228:
                    229:     /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
                    230:        steps.  Adjust the result by reducing it with the original MOD.
                    231:
                    232:        Also make sure the result is put in RES->_mp_d (where it already
                    233:        might be, see above).  */
                    234:
                    235:     if (mod_shift_cnt != 0)
                    236:       {
                    237:        carry_limb = mpn_lshift (res->_mp_d, rp, rsize, mod_shift_cnt);
                    238:        rp = res->_mp_d;
                    239:        if (carry_limb != 0)
                    240:          {
                    241:            rp[rsize] = carry_limb;
                    242:            rsize++;
                    243:          }
                    244:       }
                    245:     else
                    246:       {
                    247:        MPN_COPY (res->_mp_d, rp, rsize);
                    248:        rp = res->_mp_d;
                    249:       }
                    250:
                    251:     if (rsize >= msize)
                    252:       {
                    253:        mpn_divmod (rp + msize, rp, rsize, mp, msize);
                    254:        rsize = msize;
                    255:       }
                    256:
                    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);
                    261:   }
                    262:
                    263:   if (negative_result && rsize != 0)
                    264:     {
                    265:       if (mod_shift_cnt != 0)
                    266:        mpn_rshift (mp, mp, msize, mod_shift_cnt);
                    267:       mpn_sub (rp, mp, msize, rp, rsize);
                    268:       rsize = msize;
                    269:       MPN_NORMALIZE (rp, rsize);
                    270:     }
                    271:   res->_mp_size = rsize;
                    272:
                    273:   if (free_me != NULL)
                    274:     (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
                    275:   TMP_FREE (marker);
                    276: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>