contrib/libgmp/mpz/powm.c

   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 }