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