ghc/runtime/gmp/mpz_powm.c

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