pd.if.org Git - zpackage/blob - tomsfastmath/src/exptmod/fp_exptmod.c

   1 /* TomsFastMath, a fast ISO C bignum library.
   2  *
   3  * This project is meant to fill in where LibTomMath
   4  * falls short.  That is speed ;-)
   5  *
   6  * This project is public domain and free for all purposes.
   7  *
   8  * Tom St Denis, tomstdenis@gmail.com
   9  */
  10 #include <tfm_private.h>
  11
  12 #ifdef TFM_TIMING_RESISTANT
  13
  14 /* timing resistant montgomery ladder based exptmod
  15
  16    Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
  17 */
  18 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
  19 {
  20   fp_int   R[2];
  21   fp_digit buf, mp;
  22   int      err, bitcnt, digidx, y;
  23
  24   /* now setup montgomery  */
  25   if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
  26      return err;
  27   }
  28
  29   fp_init(&R[0]);
  30   fp_init(&R[1]);
  31
  32   /* now we need R mod m */
  33   fp_montgomery_calc_normalization (&R[0], P);
  34
  35   /* now set R[0][1] to G * R mod m */
  36   if (fp_cmp_mag(P, G) != FP_GT) {
  37      /* G > P so we reduce it first */
  38      fp_mod(G, P, &R[1]);
  39   } else {
  40      fp_copy(G, &R[1]);
  41   }
  42   fp_mulmod (&R[1], &R[0], P, &R[1]);
  43
  44   /* for j = t-1 downto 0 do
  45         r_!k = R0*R1; r_k = r_k^2
  46   */
  47
  48   /* set initial mode and bit cnt */
  49   bitcnt = 1;
  50   buf    = 0;
  51   digidx = X->used - 1;
  52
  53   for (;;) {
  54     /* grab next digit as required */
  55     if (--bitcnt == 0) {
  56       /* if digidx == -1 we are out of digits so break */
  57       if (digidx == -1) {
  58         break;
  59       }
  60       /* read next digit and reset bitcnt */
  61       buf    = X->dp[digidx--];
  62       bitcnt = (int)DIGIT_BIT;
  63     }
  64
  65     /* grab the next msb from the exponent */
  66     y     = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
  67     buf <<= (fp_digit)1;
  68
  69     /* do ops */
  70     fp_mul(&R[0], &R[1], &R[y^1]); fp_montgomery_reduce(&R[y^1], P, mp);
  71     fp_sqr(&R[y], &R[y]);          fp_montgomery_reduce(&R[y], P, mp);
  72   }
  73
  74    fp_montgomery_reduce(&R[0], P, mp);
  75    fp_copy(&R[0], Y);
  76    return FP_OKAY;
  77 }
  78
  79 #else
  80
  81 /* y = g**x (mod b)
  82  * Some restrictions... x must be positive and < b
  83  */
  84 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
  85 {
  86   fp_int   M[64], res;
  87   fp_digit buf, mp;
  88   int      err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
  89
  90   /* find window size */
  91   x = fp_count_bits (X);
  92   if (x <= 21) {
  93     winsize = 1;
  94   } else if (x <= 36) {
  95     winsize = 3;
  96   } else if (x <= 140) {
  97     winsize = 4;
  98   } else if (x <= 450) {
  99     winsize = 5;
 100   } else {
 101     winsize = 6;
 102   }
 103
 104   /* init M array */
 105   memset(M, 0, sizeof(M));
 106
 107   /* now setup montgomery  */
 108   if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
 109      return err;
 110   }
 111
 112   /* setup result */
 113   fp_init(&res);
 114
 115   /* create M table
 116    *
 117    * The M table contains powers of the input base, e.g. M[x] = G^x mod P
 118    *
 119    * The first half of the table is not computed though accept for M[0] and M[1]
 120    */
 121
 122    /* now we need R mod m */
 123    fp_montgomery_calc_normalization (&res, P);
 124
 125    /* now set M[1] to G * R mod m */
 126    if (fp_cmp_mag(P, G) != FP_GT) {
 127       /* G > P so we reduce it first */
 128       fp_mod(G, P, &M[1]);
 129    } else {
 130       fp_copy(G, &M[1]);
 131    }
 132    fp_mulmod (&M[1], &res, P, &M[1]);
 133
 134   /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
 135   fp_copy (&M[1], &M[1 << (winsize - 1)]);
 136   for (x = 0; x < (winsize - 1); x++) {
 137     fp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)]);
 138     fp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp);
 139   }
 140
 141   /* create upper table */
 142   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
 143     fp_mul(&M[x - 1], &M[1], &M[x]);
 144     fp_montgomery_reduce(&M[x], P, mp);
 145   }
 146
 147   /* set initial mode and bit cnt */
 148   mode   = 0;
 149   bitcnt = 1;
 150   buf    = 0;
 151   digidx = X->used - 1;
 152   bitcpy = 0;
 153   bitbuf = 0;
 154
 155   for (;;) {
 156     /* grab next digit as required */
 157     if (--bitcnt == 0) {
 158       /* if digidx == -1 we are out of digits so break */
 159       if (digidx == -1) {
 160         break;
 161       }
 162       /* read next digit and reset bitcnt */
 163       buf    = X->dp[digidx--];
 164       bitcnt = (int)DIGIT_BIT;
 165     }
 166
 167     /* grab the next msb from the exponent */
 168     y     = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
 169     buf <<= (fp_digit)1;
 170
 171     /* if the bit is zero and mode == 0 then we ignore it
 172      * These represent the leading zero bits before the first 1 bit
 173      * in the exponent.  Technically this opt is not required but it
 174      * does lower the # of trivial squaring/reductions used
 175      */
 176     if (mode == 0 && y == 0) {
 177       continue;
 178     }
 179
 180     /* if the bit is zero and mode == 1 then we square */
 181     if (mode == 1 && y == 0) {
 182       fp_sqr(&res, &res);
 183       fp_montgomery_reduce(&res, P, mp);
 184       continue;
 185     }
 186
 187     /* else we add it to the window */
 188     bitbuf |= (y << (winsize - ++bitcpy));
 189     mode    = 2;
 190
 191     if (bitcpy == winsize) {
 192       /* ok window is filled so square as required and multiply  */
 193       /* square first */
 194       for (x = 0; x < winsize; x++) {
 195         fp_sqr(&res, &res);
 196         fp_montgomery_reduce(&res, P, mp);
 197       }
 198
 199       /* then multiply */
 200       fp_mul(&res, &M[bitbuf], &res);
 201       fp_montgomery_reduce(&res, P, mp);
 202
 203       /* empty window and reset */
 204       bitcpy = 0;
 205       bitbuf = 0;
 206       mode   = 1;
 207     }
 208   }
 209
 210   /* if bits remain then square/multiply */
 211   if (mode == 2 && bitcpy > 0) {
 212     /* square then multiply if the bit is set */
 213     for (x = 0; x < bitcpy; x++) {
 214       fp_sqr(&res, &res);
 215       fp_montgomery_reduce(&res, P, mp);
 216
 217       /* get next bit of the window */
 218       bitbuf <<= 1;
 219       if ((bitbuf & (1 << winsize)) != 0) {
 220         /* then multiply */
 221         fp_mul(&res, &M[1], &res);
 222         fp_montgomery_reduce(&res, P, mp);
 223       }
 224     }
 225   }
 226
 227   /* fixup result if Montgomery reduction is used
 228    * recall that any value in a Montgomery system is
 229    * actually multiplied by R mod n.  So we have
 230    * to reduce one more time to cancel out the factor
 231    * of R.
 232    */
 233   fp_montgomery_reduce(&res, P, mp);
 234
 235   /* swap res with Y */
 236   fp_copy (&res, Y);
 237   return FP_OKAY;
 238 }
 239
 240 #endif
 241
 242
 243 int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
 244 {
 245    fp_int tmp;
 246    int    err;
 247
 248 #ifdef TFM_CHECK
 249    /* prevent overflows */
 250    if (P->used > (FP_SIZE/2)) {
 251       return FP_VAL;
 252    }
 253 #endif
 254
 255    /* is X negative?  */
 256    if (X->sign == FP_NEG) {
 257       /* yes, copy G and invmod it */
 258       fp_copy(G, &tmp);
 259       if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
 260          return err;
 261       }
 262       X->sign = FP_ZPOS;
 263       err =  _fp_exptmod(&tmp, X, P, Y);
 264       if (X != Y) {
 265          X->sign = FP_NEG;
 266       }
 267       return err;
 268    } else {
 269       /* Positive exponent so just exptmod */
 270       return _fp_exptmod(G, X, P, Y);
 271    }
 272 }
 273
 274 /* $Source$ */
 275 /* $Revision$ */
 276 /* $Date$ */