pb_exptmod.c

   1 /* LibTomPoly, Polynomial Basis Math -- Tom St Denis
   2  *
   3  * LibTomPoly is a public domain library that provides
   4  * polynomial basis arithmetic support.  It relies on
   5  * LibTomMath for large integer support.
   6  *
   7  * This library is free for all purposes without any
   8  * express guarantee that it works.
   9  *
  10  * Tom St Denis, tomstdenis@iahu.ca, http://poly.libtomcrypt.org
  11  */
  12 #include <tompoly.h>
  13
  14 #ifdef MP_LOW_MEM
  15    #define TAB_SIZE 32
  16 #else
  17    #define TAB_SIZE 256
  18 #endif
  19
  20 int pb_exptmod (pb_poly * G, mp_int * X, pb_poly * P, pb_poly * Y)
  21 {
  22   pb_poly  M[TAB_SIZE], res;
  23   mp_digit buf;
  24   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
  25
  26   /* find window size */
  27   x = mp_count_bits (X);
  28   if (x <= 7) {
  29     winsize = 2;
  30   } else if (x <= 36) {
  31     winsize = 3;
  32   } else if (x <= 140) {
  33     winsize = 4;
  34   } else if (x <= 450) {
  35     winsize = 5;
  36   } else if (x <= 1303) {
  37     winsize = 6;
  38   } else if (x <= 3529) {
  39     winsize = 7;
  40   } else {
  41     winsize = 8;
  42   }
  43
  44 #ifdef MP_LOW_MEM
  45     if (winsize > 5) {
  46        winsize = 5;
  47     }
  48 #endif
  49
  50   /* init M array */
  51   /* init first cell */
  52   if ((err = pb_init(&M[1], &(Y->characteristic))) != MP_OKAY) {
  53      return err;
  54   }
  55
  56   /* now init the second half of the array */
  57   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
  58     if ((err = pb_init(&M[x], &(Y->characteristic))) != MP_OKAY) {
  59       for (y = 1<<(winsize-1); y < x; y++) {
  60         pb_clear (&M[y]);
  61       }
  62       pb_clear(&M[1]);
  63       return err;
  64     }
  65   }
  66
  67   /* create M table
  68    *
  69    * The M table contains powers of the base,
  70    * e.g. M[x] = G**x mod P
  71    *
  72    * The first half of the table is not
  73    * computed though accept for M[0] and M[1]
  74    */
  75
  76   if (X->sign == MP_ZPOS) {
  77      if ((err = pb_mod (G, P, &M[1])) != MP_OKAY)                             { goto __M; }
  78   } else {
  79      if ((err = pb_invmod(G, P, &M[1])) != MP_OKAY)                           { goto __M; }
  80   }
  81
  82   /* compute the value at M[1<<(winsize-1)] by squaring
  83    * M[1] (winsize-1) times
  84    */
  85   if ((err = pb_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY)             { goto __M; }
  86
  87   for (x = 0; x < (winsize - 1); x++) {
  88     if ((err = pb_mulmod (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)],
  89                           P, &M[1 << (winsize - 1)])) != MP_OKAY)             { goto __M; }
  90   }
  91
  92   /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
  93    * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
  94    */
  95   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
  96     if ((err = pb_mulmod (&M[x - 1], &M[1], P, &M[x])) != MP_OKAY)            { goto __M; }
  97   }
  98
  99   /* setup result */
 100   if ((err = pb_init (&res, &(Y->characteristic))) != MP_OKAY)                { goto __M; }
 101   mp_set (&(res.terms[0]), 1); res.used = 1;
 102
 103   /* set initial mode and bit cnt */
 104   mode   = 0;
 105   bitcnt = 1;
 106   buf    = 0;
 107   digidx = X->used - 1;
 108   bitcpy = 0;
 109   bitbuf = 0;
 110
 111   for (;;) {
 112     /* grab next digit as required */
 113     if (--bitcnt == 0) {
 114       /* if digidx == -1 we are out of digits */
 115       if (digidx == -1) {
 116         break;
 117       }
 118       /* read next digit and reset the bitcnt */
 119       buf    = X->dp[digidx--];
 120       bitcnt = (int) DIGIT_BIT;
 121     }
 122
 123     /* grab the next msb from the exponent */
 124     y     = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
 125     buf <<= (mp_digit)1;
 126
 127     /* if the bit is zero and mode == 0 then we ignore it
 128      * These represent the leading zero bits before the first 1 bit
 129      * in the exponent.  Technically this opt is not required but it
 130      * does lower the # of trivial squaring/reductions used
 131      */
 132     if (mode == 0 && y == 0) {
 133       continue;
 134     }
 135
 136     /* if the bit is zero and mode == 1 then we square */
 137     if (mode == 1 && y == 0) {
 138       if ((err = pb_mulmod (&res, &res, P, &res)) != MP_OKAY)                 { goto __RES; }
 139       continue;
 140     }
 141
 142     /* else we add it to the window */
 143     bitbuf |= (y << (winsize - ++bitcpy));
 144     mode    = 2;
 145
 146     if (bitcpy == winsize) {
 147       /* ok window is filled so square as required and multiply  */
 148       /* square first */
 149       for (x = 0; x < winsize; x++) {
 150         if ((err = pb_mulmod (&res, &res, P, &res)) != MP_OKAY)               { goto __RES; }
 151       }
 152
 153       /* then multiply */
 154       if ((err = pb_mulmod (&res, &M[bitbuf], P, &res)) != MP_OKAY)           { goto __RES; }
 155
 156       /* empty window and reset */
 157       bitcpy = 0;
 158       bitbuf = 0;
 159       mode   = 1;
 160     }
 161   }
 162
 163   /* if bits remain then square/multiply */
 164   if (mode == 2 && bitcpy > 0) {
 165     /* square then multiply if the bit is set */
 166     for (x = 0; x < bitcpy; x++) {
 167       if ((err = pb_mulmod (&res, &res, P, &res)) != MP_OKAY)                 { goto __RES; }
 168
 169       bitbuf <<= 1;
 170       if ((bitbuf & (1 << winsize)) != 0) {
 171         /* then multiply */
 172         if ((err = pb_mulmod (&res, &M[1], P, &res)) != MP_OKAY)              { goto __RES; }
 173       }
 174     }
 175   }
 176
 177   pb_exch (&res, Y);
 178   err = MP_OKAY;
 179 __RES:pb_clear (&res);
 180 __M:
 181   pb_clear(&M[1]);
 182   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
 183     pb_clear (&M[x]);
 184   }
 185   return err;
 186 }