libs/libtommath/bn_s_mp_exptmod.c

   1 #include <tommath.h>
   2 #ifdef BN_S_MP_EXPTMOD_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   4  *
   5  * LibTomMath is a library that provides multiple-precision
   6  * integer arithmetic as well as number theoretic functionality.
   7  *
   8  * The library was designed directly after the MPI library by
   9  * Michael Fromberger but has been written from scratch with
  10  * additional optimizations in place.
  11  *
  12  * The library is free for all purposes without any express
  13  * guarantee it works.
  14  *
  15  * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16  */
  17 #ifdef MP_LOW_MEM
  18    #define TAB_SIZE 32
  19 #else
  20    #define TAB_SIZE 256
  21 #endif
  22
  23 int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
  24 {
  25   mp_int  M[TAB_SIZE], res, mu;
  26   mp_digit buf;
  27   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
  28   int (*redux)(mp_int*,mp_int*,mp_int*);
  29
  30   /* find window size */
  31   x = mp_count_bits (X);
  32   if (x <= 7) {
  33     winsize = 2;
  34   } else if (x <= 36) {
  35     winsize = 3;
  36   } else if (x <= 140) {
  37     winsize = 4;
  38   } else if (x <= 450) {
  39     winsize = 5;
  40   } else if (x <= 1303) {
  41     winsize = 6;
  42   } else if (x <= 3529) {
  43     winsize = 7;
  44   } else {
  45     winsize = 8;
  46   }
  47
  48 #ifdef MP_LOW_MEM
  49     if (winsize > 5) {
  50        winsize = 5;
  51     }
  52 #endif
  53
  54   /* init M array */
  55   /* init first cell */
  56   if ((err = mp_init(&M[1])) != MP_OKAY) {
  57      return err;
  58   }
  59
  60   /* now init the second half of the array */
  61   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
  62     if ((err = mp_init(&M[x])) != MP_OKAY) {
  63       for (y = 1<<(winsize-1); y < x; y++) {
  64         mp_clear (&M[y]);
  65       }
  66       mp_clear(&M[1]);
  67       return err;
  68     }
  69   }
  70
  71   /* create mu, used for Barrett reduction */
  72   if ((err = mp_init (&mu)) != MP_OKAY) {
  73     goto LBL_M;
  74   }
  75
  76   if (redmode == 0) {
  77      if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
  78         goto LBL_MU;
  79      }
  80      redux = mp_reduce;
  81   } else {
  82      if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
  83         goto LBL_MU;
  84      }
  85      redux = mp_reduce_2k_l;
  86   }
  87
  88   /* create M table
  89    *
  90    * The M table contains powers of the base,
  91    * e.g. M[x] = G**x mod P
  92    *
  93    * The first half of the table is not
  94    * computed though accept for M[0] and M[1]
  95    */
  96   if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
  97     goto LBL_MU;
  98   }
  99
 100   /* compute the value at M[1<<(winsize-1)] by squaring
 101    * M[1] (winsize-1) times
 102    */
 103   if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
 104     goto LBL_MU;
 105   }
 106
 107   for (x = 0; x < (winsize - 1); x++) {
 108     /* square it */
 109     if ((err = mp_sqr (&M[1 << (winsize - 1)],
 110                        &M[1 << (winsize - 1)])) != MP_OKAY) {
 111       goto LBL_MU;
 112     }
 113
 114     /* reduce modulo P */
 115     if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
 116       goto LBL_MU;
 117     }
 118   }
 119
 120   /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
 121    * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
 122    */
 123   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
 124     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
 125       goto LBL_MU;
 126     }
 127     if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
 128       goto LBL_MU;
 129     }
 130   }
 131
 132   /* setup result */
 133   if ((err = mp_init (&res)) != MP_OKAY) {
 134     goto LBL_MU;
 135   }
 136   mp_set (&res, 1);
 137
 138   /* set initial mode and bit cnt */
 139   mode   = 0;
 140   bitcnt = 1;
 141   buf    = 0;
 142   digidx = X->used - 1;
 143   bitcpy = 0;
 144   bitbuf = 0;
 145
 146   for (;;) {
 147     /* grab next digit as required */
 148     if (--bitcnt == 0) {
 149       /* if digidx == -1 we are out of digits */
 150       if (digidx == -1) {
 151         break;
 152       }
 153       /* read next digit and reset the bitcnt */
 154       buf    = X->dp[digidx--];
 155       bitcnt = (int) DIGIT_BIT;
 156     }
 157
 158     /* grab the next msb from the exponent */
 159     y     = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
 160     buf <<= (mp_digit)1;
 161
 162     /* if the bit is zero and mode == 0 then we ignore it
 163      * These represent the leading zero bits before the first 1 bit
 164      * in the exponent.  Technically this opt is not required but it
 165      * does lower the # of trivial squaring/reductions used
 166      */
 167     if (mode == 0 && y == 0) {
 168       continue;
 169     }
 170
 171     /* if the bit is zero and mode == 1 then we square */
 172     if (mode == 1 && y == 0) {
 173       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
 174         goto LBL_RES;
 175       }
 176       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 177         goto LBL_RES;
 178       }
 179       continue;
 180     }
 181
 182     /* else we add it to the window */
 183     bitbuf |= (y << (winsize - ++bitcpy));
 184     mode    = 2;
 185
 186     if (bitcpy == winsize) {
 187       /* ok window is filled so square as required and multiply  */
 188       /* square first */
 189       for (x = 0; x < winsize; x++) {
 190         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
 191           goto LBL_RES;
 192         }
 193         if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 194           goto LBL_RES;
 195         }
 196       }
 197
 198       /* then multiply */
 199       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
 200         goto LBL_RES;
 201       }
 202       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 203         goto LBL_RES;
 204       }
 205
 206       /* empty window and reset */
 207       bitcpy = 0;
 208       bitbuf = 0;
 209       mode   = 1;
 210     }
 211   }
 212
 213   /* if bits remain then square/multiply */
 214   if (mode == 2 && bitcpy > 0) {
 215     /* square then multiply if the bit is set */
 216     for (x = 0; x < bitcpy; x++) {
 217       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
 218         goto LBL_RES;
 219       }
 220       if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 221         goto LBL_RES;
 222       }
 223
 224       bitbuf <<= 1;
 225       if ((bitbuf & (1 << winsize)) != 0) {
 226         /* then multiply */
 227         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
 228           goto LBL_RES;
 229         }
 230         if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 231           goto LBL_RES;
 232         }
 233       }
 234     }
 235   }
 236
 237   mp_exch (&res, Y);
 238   err = MP_OKAY;
 239 LBL_RES:mp_clear (&res);
 240 LBL_MU:mp_clear (&mu);
 241 LBL_M:
 242   mp_clear(&M[1]);
 243   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
 244     mp_clear (&M[x]);
 245   }
 246   return err;
 247 }
 248 #endif
 249
 250 /* $Source: /cvs/libtom/libtommath/bn_s_mp_exptmod.c,v $ */
 251 /* $Revision: 1.4 $ */
 252 /* $Date: 2006/03/31 14:18:44 $ */