fp_exptmod.c

00001 /* TomsFastMath, a fast ISO C bignum library.
00002  * 
00003  * This project is meant to fill in where LibTomMath
00004  * falls short.  That is speed ;-)
00005  *
00006  * This project is public domain and free for all purposes.
00007  * 
00008  * Tom St Denis, tomstdenis@gmail.com
00009  */
00010 #include <tfm.h>
00011 
00012 #ifdef TFM_TIMING_RESISTANT
00013 
00014 /* timing resistant montgomery ladder based exptmod 
00015 
00016    Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
00017 */
00018 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
00019 {
00020   fp_int   R[2];
00021   fp_digit buf, mp;
00022   int      err, bitcnt, digidx, y;
00023 
00024   /* now setup montgomery  */
00025   if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
00026      return err;
00027   }
00028 
00029   fp_init(&R[0]);   
00030   fp_init(&R[1]);   
00031    
00032   /* now we need R mod m */
00033   fp_montgomery_calc_normalization (&R[0], P);
00034 
00035   /* now set R[0][1] to G * R mod m */
00036   if (fp_cmp_mag(P, G) != FP_GT) {
00037      /* G > P so we reduce it first */
00038      fp_mod(G, P, &R[1]);
00039   } else {
00040      fp_copy(G, &R[1]);
00041   }
00042   fp_mulmod (&R[1], &R[0], P, &R[1]);
00043 
00044   /* for j = t-1 downto 0 do
00045         r_!k = R0*R1; r_k = r_k^2
00046   */
00047   
00048   /* set initial mode and bit cnt */
00049   bitcnt = 1;
00050   buf    = 0;
00051   digidx = X->used - 1;
00052 
00053   for (;;) {
00054     /* grab next digit as required */
00055     if (--bitcnt == 0) {
00056       /* if digidx == -1 we are out of digits so break */
00057       if (digidx == -1) {
00058         break;
00059       }
00060       /* read next digit and reset bitcnt */
00061       buf    = X->dp[digidx--];
00062       bitcnt = (int)DIGIT_BIT;
00063     }
00064 
00065     /* grab the next msb from the exponent */
00066     y     = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
00067     buf <<= (fp_digit)1;
00068 
00069     /* do ops */
00070     fp_mul(&R[0], &R[1], &R[y^1]); fp_montgomery_reduce(&R[y^1], P, mp);
00071     fp_sqr(&R[y], &R[y]);          fp_montgomery_reduce(&R[y], P, mp);
00072   }
00073 
00074    fp_montgomery_reduce(&R[0], P, mp);
00075    fp_copy(&R[0], Y);
00076    return FP_OKAY;
00077 }   
00078 
00079 #else
00080 
00081 /* y = g**x (mod b) 
00082  * Some restrictions... x must be positive and < b
00083  */
00084 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
00085 {
00086   fp_int   M[64], res;
00087   fp_digit buf, mp;
00088   int      err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
00089 
00090   /* find window size */
00091   x = fp_count_bits (X);
00092   if (x <= 21) {
00093     winsize = 1;
00094   } else if (x <= 36) {
00095     winsize = 3;
00096   } else if (x <= 140) {
00097     winsize = 4;
00098   } else if (x <= 450) {
00099     winsize = 5;
00100   } else {
00101     winsize = 6;
00102   } 
00103 
00104   /* init M array */
00105   memset(M, 0, sizeof(M)); 
00106 
00107   /* now setup montgomery  */
00108   if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
00109      return err;
00110   }
00111 
00112   /* setup result */
00113   fp_init(&res);
00114 
00115   /* create M table
00116    *
00117    * The M table contains powers of the input base, e.g. M[x] = G^x mod P
00118    *
00119    * The first half of the table is not computed though accept for M[0] and M[1]
00120    */
00121 
00122    /* now we need R mod m */
00123    fp_montgomery_calc_normalization (&res, P);
00124 
00125    /* now set M[1] to G * R mod m */
00126    if (fp_cmp_mag(P, G) != FP_GT) {
00127       /* G > P so we reduce it first */
00128       fp_mod(G, P, &M[1]);
00129    } else {
00130       fp_copy(G, &M[1]);
00131    }
00132    fp_mulmod (&M[1], &res, P, &M[1]);
00133 
00134   /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
00135   fp_copy (&M[1], &M[1 << (winsize - 1)]);
00136   for (x = 0; x < (winsize - 1); x++) {
00137     fp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)]);
00138     fp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp);
00139   }
00140 
00141   /* create upper table */
00142   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
00143     fp_mul(&M[x - 1], &M[1], &M[x]);
00144     fp_montgomery_reduce(&M[x], P, mp);
00145   }
00146 
00147   /* set initial mode and bit cnt */
00148   mode   = 0;
00149   bitcnt = 1;
00150   buf    = 0;
00151   digidx = X->used - 1;
00152   bitcpy = 0;
00153   bitbuf = 0;
00154 
00155   for (;;) {
00156     /* grab next digit as required */
00157     if (--bitcnt == 0) {
00158       /* if digidx == -1 we are out of digits so break */
00159       if (digidx == -1) {
00160         break;
00161       }
00162       /* read next digit and reset bitcnt */
00163       buf    = X->dp[digidx--];
00164       bitcnt = (int)DIGIT_BIT;
00165     }
00166 
00167     /* grab the next msb from the exponent */
00168     y     = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
00169     buf <<= (fp_digit)1;
00170 
00171     /* if the bit is zero and mode == 0 then we ignore it
00172      * These represent the leading zero bits before the first 1 bit
00173      * in the exponent.  Technically this opt is not required but it
00174      * does lower the # of trivial squaring/reductions used
00175      */
00176     if (mode == 0 && y == 0) {
00177       continue;
00178     }
00179 
00180     /* if the bit is zero and mode == 1 then we square */
00181     if (mode == 1 && y == 0) {
00182       fp_sqr(&res, &res);
00183       fp_montgomery_reduce(&res, P, mp);
00184       continue;
00185     }
00186 
00187     /* else we add it to the window */
00188     bitbuf |= (y << (winsize - ++bitcpy));
00189     mode    = 2;
00190 
00191     if (bitcpy == winsize) {
00192       /* ok window is filled so square as required and multiply  */
00193       /* square first */
00194       for (x = 0; x < winsize; x++) {
00195         fp_sqr(&res, &res);
00196         fp_montgomery_reduce(&res, P, mp);
00197       }
00198 
00199       /* then multiply */
00200       fp_mul(&res, &M[bitbuf], &res);
00201       fp_montgomery_reduce(&res, P, mp);
00202 
00203       /* empty window and reset */
00204       bitcpy = 0;
00205       bitbuf = 0;
00206       mode   = 1;
00207     }
00208   }
00209 
00210   /* if bits remain then square/multiply */
00211   if (mode == 2 && bitcpy > 0) {
00212     /* square then multiply if the bit is set */
00213     for (x = 0; x < bitcpy; x++) {
00214       fp_sqr(&res, &res);
00215       fp_montgomery_reduce(&res, P, mp);
00216 
00217       /* get next bit of the window */
00218       bitbuf <<= 1;
00219       if ((bitbuf & (1 << winsize)) != 0) {
00220         /* then multiply */
00221         fp_mul(&res, &M[1], &res);
00222         fp_montgomery_reduce(&res, P, mp);
00223       }
00224     }
00225   }
00226 
00227   /* fixup result if Montgomery reduction is used
00228    * recall that any value in a Montgomery system is
00229    * actually multiplied by R mod n.  So we have
00230    * to reduce one more time to cancel out the factor
00231    * of R.
00232    */
00233   fp_montgomery_reduce(&res, P, mp);
00234 
00235   /* swap res with Y */
00236   fp_copy (&res, Y);
00237   return FP_OKAY;
00238 }
00239 
00240 #endif
00241 
00242 
00243 int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
00244 {
00245    fp_int tmp;
00246    int    err;
00247    
00248 #ifdef TFM_CHECK
00249    /* prevent overflows */
00250    if (P->used > (FP_SIZE/2)) {
00251       return FP_VAL;
00252    }
00253 #endif
00254 
00255    /* is X negative?  */
00256    if (X->sign == FP_NEG) {
00257       /* yes, copy G and invmod it */
00258       fp_copy(G, &tmp);
00259       if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
00260          return err;
00261       }
00262       X->sign = FP_ZPOS;
00263       err =  _fp_exptmod(&tmp, X, P, Y);
00264       if (X != Y) {
00265          X->sign = FP_NEG;
00266       }
00267       return err;
00268    } else {
00269       /* Positive exponent so just exptmod */
00270       return _fp_exptmod(G, X, P, Y);
00271    }
00272 }
00273 
00274 /* $Source: /cvs/libtom/tomsfastmath/src/exptmod/fp_exptmod.c,v $ */
00275 /* $Revision: 1.1 $ */
00276 /* $Date: 2006/12/31 21:25:53 $ */