ghc/runtime/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, 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 void
  26 #ifdef __STDC__
  27 mpz_powm_ui (MP_INT *res, const MP_INT *base, unsigned long int exp,
  28              const MP_INT *mod)
  29 #else
  30 mpz_powm_ui (res, base, exp, mod)
  31      MP_INT *res;
  32      const MP_INT *base;
  33      unsigned long int exp;
  34      const MP_INT *mod;
  35 #endif
  36 {
  37   mp_ptr rp, mp, bp;
  38   mp_size msize, bsize, rsize;
  39   mp_size size;
  40   int mod_shift_cnt;
  41   int negative_result;
  42   mp_limb *free_me = NULL;
  43   size_t free_me_size;
  44
  45   msize = ABS (mod->size);
  46   size = 2 * msize;
  47
  48   rp = res->d;
  49
  50   /* Normalize MOD (i.e. make its most significant bit set) as required by
  51      mpn_div.  This will make the intermediate values in the calculation
  52      slightly larger, but the correct result is obtained after a final
  53      reduction using the original MOD value.  */
  54
  55   mp = (mp_ptr) alloca (msize * BYTES_PER_MP_LIMB);
  56   count_leading_zeros (mod_shift_cnt, mod->d[msize - 1]);
  57   if (mod_shift_cnt != 0)
  58     (void) mpn_lshift (mp, mod->d, msize, mod_shift_cnt);
  59   else
  60     MPN_COPY (mp, mod->d, msize);
  61
  62   bsize = ABS (base->size);
  63   if (bsize > msize)
  64     {
  65       /* The base is larger than the module.  Reduce it.  */
  66
  67       /* Allocate (BSIZE + 1) with space for remainder and quotient.
  68          (The quotient is (bsize - msize + 1) limbs.)  */
  69       bp = (mp_ptr) alloca ((bsize + 1) * BYTES_PER_MP_LIMB);
  70       MPN_COPY (bp, base->d, bsize);
  71       /* We don't care about the quotient, store it above the remainder,
  72          at BP + MSIZE.  */
  73       mpn_div (bp + msize, bp, bsize, mp, msize);
  74       bsize = msize;
  75       while (bsize > 0 && bp[bsize - 1] == 0)
  76         bsize--;
  77     }
  78   else
  79     {
  80       bp = base->d;
  81       bsize = ABS (base->size);
  82     }
  83
  84   if (res->alloc < size)
  85     {
  86       /* We have to allocate more space for RES.  If any of the input
  87          parameters are identical to RES, defer deallocation of the old
  88          space.  */
  89
  90       if (rp == mp || rp == bp)
  91         {
  92           free_me = rp;
  93           free_me_size = res->alloc;
  94         }
  95       else
  96         (*_mp_free_func) (rp, res->alloc * BYTES_PER_MP_LIMB);
  97
  98       rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB);
  99       res->alloc = size;
 100       res->d = rp;
 101     }
 102   else
 103     {
 104       /* Make BASE, EXP and MOD not overlap with RES.  */
 105       if (rp == bp)
 106         {
 107           /* RES and BASE are identical.  Allocate temp. space for BASE.  */
 108           bp = (mp_ptr) alloca (bsize * BYTES_PER_MP_LIMB);
 109           MPN_COPY (bp, rp, bsize);
 110         }
 111       if (rp == mp)
 112         {
 113           /* RES and MOD are identical.  Allocate temporary space for MOD.  */
 114           mp = (mp_ptr) alloca (msize * BYTES_PER_MP_LIMB);
 115           MPN_COPY (mp, rp, msize);
 116         }
 117     }
 118
 119   if (exp == 0)
 120     {
 121       rp[0] = 1;
 122       res->size = 1;
 123       return;
 124     }
 125
 126   MPN_COPY (rp, bp, bsize);
 127   rsize = bsize;
 128
 129   {
 130     mp_size xsize;
 131     mp_ptr dummyp = (mp_ptr) alloca ((msize + 1) * BYTES_PER_MP_LIMB);
 132     mp_ptr xp = (mp_ptr) alloca (2 * (msize + 1) * BYTES_PER_MP_LIMB);
 133     int c;
 134     mp_limb e;
 135     mp_limb carry_limb;
 136
 137     negative_result = (exp & 1) && base->size < 0;
 138
 139     e = exp;
 140     count_leading_zeros (c, e);
 141     e <<= (c + 1);              /* shift the exp bits to the left, loose msb */
 142     c = BITS_PER_MP_LIMB - 1 - c;
 143
 144     /* Main loop.
 145
 146        Make the result be pointed to alternately by XP and RP.  This
 147        helps us avoid block copying, which would otherwise be necessary
 148        with the overlap restrictions of mpn_div.  With 50% probability
 149        the result after this loop will be in the area originally pointed
 150        by RP (==RES->D), and with 50% probability in the area originally
 151        pointed to by XP.  */
 152
 153     while (c != 0)
 154       {
 155         mp_ptr tp;
 156         mp_size tsize;
 157
 158         xsize = mpn_mul (xp, rp, rsize, rp, rsize);
 159         mpn_div (dummyp, xp, xsize, mp, msize);
 160
 161         /* Remove any leading zero words from the result.  */
 162         if (xsize > msize)
 163           xsize = msize;
 164         while (xsize > 0 && xp[xsize - 1] == 0)
 165           xsize--;
 166
 167         tp = rp; rp = xp; xp = tp;
 168         tsize = rsize; rsize = xsize; xsize = tsize;
 169
 170         if ((mp_limb_signed) e < 0)
 171           {
 172             if (rsize > bsize)
 173               xsize = mpn_mul (xp, rp, rsize, bp, bsize);
 174             else
 175               xsize = mpn_mul (xp, bp, bsize, rp, rsize);
 176             mpn_div (dummyp, xp, xsize, mp, msize);
 177
 178             /* Remove any leading zero words from the result.  */
 179             if (xsize > msize)
 180               xsize = msize;
 181             while (xsize > 0 && xp[xsize - 1] == 0)
 182               xsize--;
 183
 184             tp = rp; rp = xp; xp = tp;
 185             tsize = rsize; rsize = xsize; xsize = tsize;
 186           }
 187         e <<= 1;
 188         c--;
 189       }
 190
 191     /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT
 192        steps.  Adjust the result by reducing it with the original MOD.
 193
 194        Also make sure the result is put in RES->D (where it already
 195        might be, see above).  */
 196
 197     carry_limb = mpn_lshift (res->d, rp, rsize, mod_shift_cnt);
 198     rp = res->d;
 199     if (carry_limb != 0)
 200       {
 201         rp[rsize] = carry_limb;
 202         rsize++;
 203       }
 204     mpn_div (dummyp, rp, rsize, mp, msize);
 205     /* Remove any leading zero words from the result.  */
 206     if (rsize > msize)
 207       rsize = msize;
 208     while (rsize > 0 && rp[rsize - 1] == 0)
 209       rsize--;
 210     rsize = mpn_rshift (rp, rp, rsize, mod_shift_cnt);
 211   }
 212
 213   res->size = negative_result >= 0 ?  rsize : -rsize;
 214
 215   if (free_me != NULL)
 216     (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
 217
 218   alloca (0);
 219 }