ghc/rts/gmp/mpz/powm_ui.c

   1 /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
   2
   3 Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation,
   4 Inc.
   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
   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
  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
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  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
  23 #include <stdio.h> /* for NULL */
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27
  28 void
  29 #if __STDC__
  30 mpz_powm_ui (mpz_ptr res, mpz_srcptr base, unsigned long int exp, mpz_srcptr mod)
  31 #else
  32 mpz_powm_ui (res, base, exp, mod)
  33      mpz_ptr res;
  34      mpz_srcptr base;
  35      unsigned long int exp;
  36      mpz_srcptr mod;
  37 #endif
  38 {
  39   mp_ptr rp, mp, bp;
  40   mp_size_t msize, bsize, rsize;
  41   mp_size_t size;
  42   int mod_shift_cnt;
  43   int negative_result;
  44   mp_limb_t *free_me = NULL;
  45   size_t free_me_size;
  46   TMP_DECL (marker);
  47
  48   msize = ABS (mod->_mp_size);
  49   size = 2 * msize;
  50
  51   rp = res->_mp_d;
  52
  53   if (msize == 0)
  54     DIVIDE_BY_ZERO;
  55
  56   if (exp == 0)
  57     {
  58       /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
  59          depending on if MOD equals 1.  */
  60       res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1;
  61       rp[0] = 1;
  62       return;
  63     }
  64
  65   TMP_MARK (marker);
  66
  67   /* Normalize MOD (i.e. make its most significant bit set) as required by
  68      mpn_divmod.  This will make the intermediate values in the calculation
  69      slightly larger, but the correct result is obtained after a final
  70      reduction using the original MOD value.  */
  71
  72   mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB);
  73   count_leading_zeros (mod_shift_cnt, mod->_mp_d[msize - 1]);
  74   if (mod_shift_cnt != 0)
  75     mpn_lshift (mp, mod->_mp_d, msize, mod_shift_cnt);
  76   else
  77     MPN_COPY (mp, mod->_mp_d, msize);
  78
  79   bsize = ABS (base->_mp_size);
  80   if (bsize > msize)
  81     {
  82       /* The base is larger than the module.  Reduce it.  */
  83
  84       /* Allocate (BSIZE + 1) with space for remainder and quotient.
  85          (The quotient is (bsize - msize + 1) limbs.)  */
  86       bp = (mp_ptr) TMP_ALLOC ((bsize + 1) * BYTES_PER_MP_LIMB);
  87       MPN_COPY (bp, base->_mp_d, bsize);
  88       /* We don't care about the quotient, store it above the remainder,
  89          at BP + MSIZE.  */
  90       mpn_divmod (bp + msize, bp, bsize, mp, msize);
  91       bsize = msize;
  92       /* Canonicalize the base, since we are going to multiply with it
  93          quite a few times.  */
  94       MPN_NORMALIZE (bp, bsize);
  95     }
  96   else
  97     bp = base->_mp_d;
  98
  99   if (bsize == 0)
 100     {
 101       res->_mp_size = 0;
 102       TMP_FREE (marker);
 103       return;
 104     }
 105
 106   if (res->_mp_alloc < size)
 107     {
 108       /* We have to allocate more space for RES.  If any of the input
 109          parameters are identical to RES, defer deallocation of the old
 110          space.  */
 111
 112       if (rp == mp || rp == bp)
 113         {
 114           free_me = rp;
 115           free_me_size = res->_mp_alloc;
 116         }
 117       else
 118         (*_mp_free_func) (rp, res->_mp_alloc * BYTES_PER_MP_LIMB);
 119
 120       rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB);
 121       res->_mp_alloc = size;
 122       res->_mp_d = rp;
 123     }
 124   else
 125     {
 126       /* Make BASE, EXP and MOD not overlap with RES.  */
 127       if (rp == bp)
 128         {
 129           /* RES and BASE are identical.  Allocate temp. space for BASE.  */
 130           bp = (mp_ptr) TMP_ALLOC (bsize * BYTES_PER_MP_LIMB);
 131           MPN_COPY (bp, rp, bsize);
 132         }
 133       if (rp == mp)
 134         {
 135           /* RES and MOD are identical.  Allocate temporary space for MOD.  */
 136           mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB);
 137           MPN_COPY (mp, rp, msize);
 138         }
 139     }
 140
 141   MPN_COPY (rp, bp, bsize);
 142   rsize = bsize;
 143
 144   {
 145     mp_ptr xp = (mp_ptr) TMP_ALLOC (2 * (msize + 1) * BYTES_PER_MP_LIMB);
 146     int c;
 147     mp_limb_t e;
 148     mp_limb_t carry_limb;
 149
 150     negative_result = (exp & 1) && base->_mp_size < 0;
 151
 152     e = exp;
 153     count_leading_zeros (c, e);
 154     e = (e << c) << 1;          /* shift the exp bits to the left, lose msb */
 155     c = BITS_PER_MP_LIMB - 1 - c;
 156
 157     /* Main loop.
 158
 159        Make the result be pointed to alternately by XP and RP.  This
 160        helps us avoid block copying, which would otherwise be necessary
 161        with the overlap restrictions of mpn_divmod.  With 50% probability
 162        the result after this loop will be in the area originally pointed
 163        by RP (==RES->_mp_d), and with 50% probability in the area originally
 164        pointed to by XP.  */
 165
 166     while (c != 0)
 167       {
 168         mp_ptr tp;
 169         mp_size_t xsize;
 170
 171         mpn_mul_n (xp, rp, rp, rsize);
 172         xsize = 2 * rsize;
 173         xsize -= xp[xsize - 1] == 0;
 174         if (xsize > msize)
 175           {
 176             mpn_divmod (xp + msize, xp, xsize, mp, msize);
 177             xsize = msize;
 178           }
 179
 180         tp = rp; rp = xp; xp = tp;
 181         rsize = xsize;
 182
 183         if ((mp_limb_signed_t) e < 0)
 184           {
 185             mpn_mul (xp, rp, rsize, bp, bsize);
 186             xsize = rsize + bsize;
 187             xsize -= xp[xsize - 1] == 0;
 188             if (xsize > msize)
 189               {
 190                 mpn_divmod (xp + msize, xp, xsize, mp, msize);
 191                 xsize = msize;
 192               }
 193
 194             tp = rp; rp = xp; xp = tp;
 195             rsize = xsize;
 196           }
 197         e <<= 1;
 198         c--;
 199       }
 200
 201     /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
 202        steps.  Adjust the result by reducing it with the original MOD.
 203
 204        Also make sure the result is put in RES->_mp_d (where it already
 205        might be, see above).  */
 206
 207     if (mod_shift_cnt != 0)
 208       {
 209         carry_limb = mpn_lshift (res->_mp_d, rp, rsize, mod_shift_cnt);
 210         rp = res->_mp_d;
 211         if (carry_limb != 0)
 212           {
 213             rp[rsize] = carry_limb;
 214             rsize++;
 215           }
 216       }
 217     else
 218       {
 219         MPN_COPY (res->_mp_d, rp, rsize);
 220         rp = res->_mp_d;
 221       }
 222
 223     if (rsize >= msize)
 224       {
 225         mpn_divmod (rp + msize, rp, rsize, mp, msize);
 226         rsize = msize;
 227       }
 228
 229     /* Remove any leading zero words from the result.  */
 230     if (mod_shift_cnt != 0)
 231       mpn_rshift (rp, rp, rsize, mod_shift_cnt);
 232     MPN_NORMALIZE (rp, rsize);
 233   }
 234
 235   if (negative_result && rsize != 0)
 236     {
 237       if (mod_shift_cnt != 0)
 238         mpn_rshift (mp, mp, msize, mod_shift_cnt);
 239       mpn_sub (rp, mp, msize, rp, rsize);
 240       rsize = msize;
 241       MPN_NORMALIZE (rp, rsize);
 242     }
 243   res->_mp_size = rsize;
 244
 245   if (free_me != NULL)
 246     (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
 247   TMP_FREE (marker);
 248 }