contrib/gmp/mpz/powm.c

   1 /* mpz_powm(res,base,exp,mod) -- Set R to (U^E) mod M.
   2
   3    Contributed to the GNU project by Torbjorn Granlund.
   4
   5 Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2005, 2008, 2009
   6 Free Software Foundation, Inc.
   7
   8 This file is part of the GNU MP Library.
   9
  10 The GNU MP Library is free software; you can redistribute it and/or modify
  11 it under the terms of the GNU Lesser General Public License as published by
  12 the Free Software Foundation; either version 3 of the License, or (at your
  13 option) any later version.
  14
  15 The GNU MP Library is distributed in the hope that it will be useful, but
  16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  18 License for more details.
  19
  20 You should have received a copy of the GNU Lesser General Public License
  21 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  22
  23
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27 #ifdef BERKELEY_MP
  28 #include "mp.h"
  29 #endif
  30
  31
  32 /* TODO
  33
  34  * Improve handling of buffers.  It is pretty ugly now.
  35
  36  * For even moduli, we compute a binvert of its odd part both here and in
  37    mpn_powm.  How can we avoid this recomputation?
  38 */
  39
  40 /*
  41   b ^ e mod m   res
  42   0   0     0    ?
  43   0   e     0    ?
  44   0   0     m    ?
  45   0   e     m    0
  46   b   0     0    ?
  47   b   e     0    ?
  48   b   0     m    1 mod m
  49   b   e     m    b^e mod m
  50 */
  51
  52 #define HANDLE_NEGATIVE_EXPONENT 1
  53
  54 void
  55 #ifndef BERKELEY_MP
  56 mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m)
  57 #else /* BERKELEY_MP */
  58 pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r)
  59 #endif /* BERKELEY_MP */
  60 {
  61   mp_size_t n, nodd, ncnt;
  62   int cnt;
  63   mp_ptr rp, tp;
  64   mp_srcptr bp, ep, mp;
  65   mp_size_t rn, bn, es, en, itch;
  66   TMP_DECL;
  67
  68   n = ABSIZ(m);
  69   if (n == 0)
  70     DIVIDE_BY_ZERO;
  71
  72   mp = PTR(m);
  73
  74   TMP_MARK;
  75
  76   es = SIZ(e);
  77   if (UNLIKELY (es <= 0))
  78     {
  79       mpz_t new_b;
  80       if (es == 0)
  81         {
  82           /* b^0 mod m,  b is anything and m is non-zero.
  83              Result is 1 mod m, i.e., 1 or 0 depending on if m = 1.  */
  84           SIZ(r) = n != 1 || mp[0] != 1;
  85           PTR(r)[0] = 1;
  86           TMP_FREE;     /* we haven't really allocated anything here */
  87           return;
  88         }
  89 #if HANDLE_NEGATIVE_EXPONENT
  90       MPZ_TMP_INIT (new_b, n + 1);
  91
  92       if (! mpz_invert (new_b, b, m))
  93         DIVIDE_BY_ZERO;
  94       b = new_b;
  95       es = -es;
  96 #else
  97       DIVIDE_BY_ZERO;
  98 #endif
  99     }
 100   en = es;
 101
 102   bn = ABSIZ(b);
 103
 104   if (UNLIKELY (bn == 0))
 105     {
 106       SIZ(r) = 0;
 107       TMP_FREE;
 108       return;
 109     }
 110
 111   ep = PTR(e);
 112
 113   /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case.  */
 114   if (UNLIKELY (en == 1 && ep[0] == 1))
 115     {
 116       rp = TMP_ALLOC_LIMBS (n);
 117       bp = PTR(b);
 118       if (bn >= n)
 119         {
 120           mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1);
 121           mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n);
 122           rn = n;
 123           MPN_NORMALIZE (rp, rn);
 124
 125           if (SIZ(b) < 0 && rn != 0)
 126             {
 127               mpn_sub (rp, mp, n, rp, rn);
 128               rn = n;
 129               MPN_NORMALIZE (rp, rn);
 130             }
 131         }
 132       else
 133         {
 134           if (SIZ(b) < 0)
 135             {
 136               mpn_sub (rp, mp, n, bp, bn);
 137               rn = n;
 138               rn -= (rp[rn - 1] == 0);
 139             }
 140           else
 141             {
 142               MPN_COPY (rp, bp, bn);
 143               rn = bn;
 144             }
 145         }
 146       goto ret;
 147     }
 148
 149   /* Remove low zero limbs from M.  This loop will terminate for correctly
 150      represented mpz numbers.  */
 151   ncnt = 0;
 152   while (UNLIKELY (mp[0] == 0))
 153     {
 154       mp++;
 155       ncnt++;
 156     }
 157   nodd = n - ncnt;
 158   cnt = 0;
 159   if (mp[0] % 2 == 0)
 160     {
 161       mp_ptr new = TMP_ALLOC_LIMBS (nodd);
 162       count_trailing_zeros (cnt, mp[0]);
 163       mpn_rshift (new, mp, nodd, cnt);
 164       nodd -= new[nodd - 1] == 0;
 165       mp = new;
 166       ncnt++;
 167     }
 168
 169   if (ncnt != 0)
 170     {
 171       /* We will call both mpn_powm and mpn_powlo.  */
 172       /* rp needs n, mpn_powlo needs 4n, the 2 mpn_binvert might need more */
 173       mp_size_t n_largest_binvert = MAX (ncnt, nodd);
 174       mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert);
 175       itch = 3 * n + MAX (itch_binvert, 2 * n);
 176     }
 177   else
 178     {
 179       /* We will call just mpn_powm.  */
 180       mp_size_t itch_binvert = mpn_binvert_itch (nodd);
 181       itch = n + MAX (itch_binvert, 2 * n);
 182     }
 183   tp = TMP_ALLOC_LIMBS (itch);
 184
 185   rp = tp;  tp += n;
 186
 187   bp = PTR(b);
 188   mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp);
 189
 190   rn = n;
 191
 192   if (ncnt != 0)
 193     {
 194       mp_ptr r2, xp, yp, odd_inv_2exp;
 195       unsigned long t;
 196       int bcnt;
 197
 198       if (bn < ncnt)
 199         {
 200           mp_ptr new = TMP_ALLOC_LIMBS (ncnt);
 201           MPN_COPY (new, bp, bn);
 202           MPN_ZERO (new + bn, ncnt - bn);
 203           bp = new;
 204         }
 205
 206       r2 = tp;
 207
 208       if (bp[0] % 2 == 0)
 209         {
 210           if (en > 1)
 211             {
 212               MPN_ZERO (r2, ncnt);
 213               goto zero;
 214             }
 215
 216           ASSERT (en == 1);
 217           t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt;
 218
 219           /* Count number of low zero bits in B, up to 3.  */
 220           bcnt = (0x1213 >> ((bp[0] & 7) << 1)) & 0x3;
 221           /* Note that ep[0] * bcnt might overflow, but that just results
 222              in a missed optimization.  */
 223           if (ep[0] * bcnt >= t)
 224             {
 225               MPN_ZERO (r2, ncnt);
 226               goto zero;
 227             }
 228         }
 229
 230       mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt);
 231
 232     zero:
 233       if (nodd < ncnt)
 234         {
 235           mp_ptr new = TMP_ALLOC_LIMBS (ncnt);
 236           MPN_COPY (new, mp, nodd);
 237           MPN_ZERO (new + nodd, ncnt - nodd);
 238           mp = new;
 239         }
 240
 241       odd_inv_2exp = tp + n;
 242       mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n);
 243
 244       mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd);
 245
 246       xp = tp + 2 * n;
 247       mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt);
 248
 249       if (cnt != 0)
 250         xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1;
 251
 252       yp = tp;
 253       if (ncnt > nodd)
 254         mpn_mul (yp, xp, ncnt, mp, nodd);
 255       else
 256         mpn_mul (yp, mp, nodd, xp, ncnt);
 257
 258       mpn_add (rp, yp, n, rp, nodd);
 259
 260       ASSERT (nodd + ncnt >= n);
 261       ASSERT (nodd + ncnt <= n + 1);
 262     }
 263
 264   MPN_NORMALIZE (rp, rn);
 265
 266   if ((ep[0] & 1) && SIZ(b) < 0 && rn != 0)
 267     {
 268       mpn_sub (rp, PTR(m), n, rp, rn);
 269       rn = n;
 270       MPN_NORMALIZE (rp, rn);
 271     }
 272
 273  ret:
 274   MPZ_REALLOC (r, rn);
 275   SIZ(r) = rn;
 276   MPN_COPY (PTR(r), rp, rn);
 277
 278   TMP_FREE;
 279 }