255 lines
		
	
	
		
			6.3 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
			
		
		
	
	
			255 lines
		
	
	
		
			6.3 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
#include <tommath_private.h>
 | 
						|
#ifdef BN_S_MP_EXPTMOD_C
 | 
						|
/* LibTomMath, multiple-precision integer library -- Tom St Denis
 | 
						|
 *
 | 
						|
 * LibTomMath is a library that provides multiple-precision
 | 
						|
 * integer arithmetic as well as number theoretic functionality.
 | 
						|
 *
 | 
						|
 * The library was designed directly after the MPI library by
 | 
						|
 * Michael Fromberger but has been written from scratch with
 | 
						|
 * additional optimizations in place.
 | 
						|
 *
 | 
						|
 * The library is free for all purposes without any express
 | 
						|
 * guarantee it works.
 | 
						|
 *
 | 
						|
 * Tom St Denis, tstdenis82@gmail.com, http://libtom.org
 | 
						|
 */
 | 
						|
#ifdef MP_LOW_MEM
 | 
						|
#   define TAB_SIZE 32
 | 
						|
#else
 | 
						|
#   define TAB_SIZE 256
 | 
						|
#endif
 | 
						|
 | 
						|
int s_mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
 | 
						|
{
 | 
						|
   mp_int  M[TAB_SIZE], res, mu;
 | 
						|
   mp_digit buf;
 | 
						|
   int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
 | 
						|
   int (*redux)(mp_int *, const mp_int *, const mp_int *);
 | 
						|
 | 
						|
   /* find window size */
 | 
						|
   x = mp_count_bits(X);
 | 
						|
   if (x <= 7) {
 | 
						|
      winsize = 2;
 | 
						|
   } else if (x <= 36) {
 | 
						|
      winsize = 3;
 | 
						|
   } else if (x <= 140) {
 | 
						|
      winsize = 4;
 | 
						|
   } else if (x <= 450) {
 | 
						|
      winsize = 5;
 | 
						|
   } else if (x <= 1303) {
 | 
						|
      winsize = 6;
 | 
						|
   } else if (x <= 3529) {
 | 
						|
      winsize = 7;
 | 
						|
   } else {
 | 
						|
      winsize = 8;
 | 
						|
   }
 | 
						|
 | 
						|
#ifdef MP_LOW_MEM
 | 
						|
   if (winsize > 5) {
 | 
						|
      winsize = 5;
 | 
						|
   }
 | 
						|
#endif
 | 
						|
 | 
						|
   /* init M array */
 | 
						|
   /* init first cell */
 | 
						|
   if ((err = mp_init(&M[1])) != MP_OKAY) {
 | 
						|
      return err;
 | 
						|
   }
 | 
						|
 | 
						|
   /* now init the second half of the array */
 | 
						|
   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
 | 
						|
      if ((err = mp_init(&M[x])) != MP_OKAY) {
 | 
						|
         for (y = 1<<(winsize-1); y < x; y++) {
 | 
						|
            mp_clear(&M[y]);
 | 
						|
         }
 | 
						|
         mp_clear(&M[1]);
 | 
						|
         return err;
 | 
						|
      }
 | 
						|
   }
 | 
						|
 | 
						|
   /* create mu, used for Barrett reduction */
 | 
						|
   if ((err = mp_init(&mu)) != MP_OKAY) {
 | 
						|
      goto LBL_M;
 | 
						|
   }
 | 
						|
 | 
						|
   if (redmode == 0) {
 | 
						|
      if ((err = mp_reduce_setup(&mu, P)) != MP_OKAY) {
 | 
						|
         goto LBL_MU;
 | 
						|
      }
 | 
						|
      redux = mp_reduce;
 | 
						|
   } else {
 | 
						|
      if ((err = mp_reduce_2k_setup_l(P, &mu)) != MP_OKAY) {
 | 
						|
         goto LBL_MU;
 | 
						|
      }
 | 
						|
      redux = mp_reduce_2k_l;
 | 
						|
   }
 | 
						|
 | 
						|
   /* create M table
 | 
						|
    *
 | 
						|
    * The M table contains powers of the base,
 | 
						|
    * e.g. M[x] = G**x mod P
 | 
						|
    *
 | 
						|
    * The first half of the table is not
 | 
						|
    * computed though accept for M[0] and M[1]
 | 
						|
    */
 | 
						|
   if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
 | 
						|
      goto LBL_MU;
 | 
						|
   }
 | 
						|
 | 
						|
   /* compute the value at M[1<<(winsize-1)] by squaring
 | 
						|
    * M[1] (winsize-1) times
 | 
						|
    */
 | 
						|
   if ((err = mp_copy(&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
 | 
						|
      goto LBL_MU;
 | 
						|
   }
 | 
						|
 | 
						|
   for (x = 0; x < (winsize - 1); x++) {
 | 
						|
      /* square it */
 | 
						|
      if ((err = mp_sqr(&M[1 << (winsize - 1)],
 | 
						|
                        &M[1 << (winsize - 1)])) != MP_OKAY) {
 | 
						|
         goto LBL_MU;
 | 
						|
      }
 | 
						|
 | 
						|
      /* reduce modulo P */
 | 
						|
      if ((err = redux(&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
 | 
						|
         goto LBL_MU;
 | 
						|
      }
 | 
						|
   }
 | 
						|
 | 
						|
   /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
 | 
						|
    * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
 | 
						|
    */
 | 
						|
   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
 | 
						|
      if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
 | 
						|
         goto LBL_MU;
 | 
						|
      }
 | 
						|
      if ((err = redux(&M[x], P, &mu)) != MP_OKAY) {
 | 
						|
         goto LBL_MU;
 | 
						|
      }
 | 
						|
   }
 | 
						|
 | 
						|
   /* setup result */
 | 
						|
   if ((err = mp_init(&res)) != MP_OKAY) {
 | 
						|
      goto LBL_MU;
 | 
						|
   }
 | 
						|
   mp_set(&res, 1);
 | 
						|
 | 
						|
   /* set initial mode and bit cnt */
 | 
						|
   mode   = 0;
 | 
						|
   bitcnt = 1;
 | 
						|
   buf    = 0;
 | 
						|
   digidx = X->used - 1;
 | 
						|
   bitcpy = 0;
 | 
						|
   bitbuf = 0;
 | 
						|
 | 
						|
   for (;;) {
 | 
						|
      /* grab next digit as required */
 | 
						|
      if (--bitcnt == 0) {
 | 
						|
         /* if digidx == -1 we are out of digits */
 | 
						|
         if (digidx == -1) {
 | 
						|
            break;
 | 
						|
         }
 | 
						|
         /* read next digit and reset the bitcnt */
 | 
						|
         buf    = X->dp[digidx--];
 | 
						|
         bitcnt = (int)DIGIT_BIT;
 | 
						|
      }
 | 
						|
 | 
						|
      /* grab the next msb from the exponent */
 | 
						|
      y     = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
 | 
						|
      buf <<= (mp_digit)1;
 | 
						|
 | 
						|
      /* if the bit is zero and mode == 0 then we ignore it
 | 
						|
       * These represent the leading zero bits before the first 1 bit
 | 
						|
       * in the exponent.  Technically this opt is not required but it
 | 
						|
       * does lower the # of trivial squaring/reductions used
 | 
						|
       */
 | 
						|
      if ((mode == 0) && (y == 0)) {
 | 
						|
         continue;
 | 
						|
      }
 | 
						|
 | 
						|
      /* if the bit is zero and mode == 1 then we square */
 | 
						|
      if ((mode == 1) && (y == 0)) {
 | 
						|
         if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
 | 
						|
            goto LBL_RES;
 | 
						|
         }
 | 
						|
         if ((err = redux(&res, P, &mu)) != MP_OKAY) {
 | 
						|
            goto LBL_RES;
 | 
						|
         }
 | 
						|
         continue;
 | 
						|
      }
 | 
						|
 | 
						|
      /* else we add it to the window */
 | 
						|
      bitbuf |= (y << (winsize - ++bitcpy));
 | 
						|
      mode    = 2;
 | 
						|
 | 
						|
      if (bitcpy == winsize) {
 | 
						|
         /* ok window is filled so square as required and multiply  */
 | 
						|
         /* square first */
 | 
						|
         for (x = 0; x < winsize; x++) {
 | 
						|
            if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
 | 
						|
               goto LBL_RES;
 | 
						|
            }
 | 
						|
            if ((err = redux(&res, P, &mu)) != MP_OKAY) {
 | 
						|
               goto LBL_RES;
 | 
						|
            }
 | 
						|
         }
 | 
						|
 | 
						|
         /* then multiply */
 | 
						|
         if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
 | 
						|
            goto LBL_RES;
 | 
						|
         }
 | 
						|
         if ((err = redux(&res, P, &mu)) != MP_OKAY) {
 | 
						|
            goto LBL_RES;
 | 
						|
         }
 | 
						|
 | 
						|
         /* empty window and reset */
 | 
						|
         bitcpy = 0;
 | 
						|
         bitbuf = 0;
 | 
						|
         mode   = 1;
 | 
						|
      }
 | 
						|
   }
 | 
						|
 | 
						|
   /* if bits remain then square/multiply */
 | 
						|
   if ((mode == 2) && (bitcpy > 0)) {
 | 
						|
      /* square then multiply if the bit is set */
 | 
						|
      for (x = 0; x < bitcpy; x++) {
 | 
						|
         if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
 | 
						|
            goto LBL_RES;
 | 
						|
         }
 | 
						|
         if ((err = redux(&res, P, &mu)) != MP_OKAY) {
 | 
						|
            goto LBL_RES;
 | 
						|
         }
 | 
						|
 | 
						|
         bitbuf <<= 1;
 | 
						|
         if ((bitbuf & (1 << winsize)) != 0) {
 | 
						|
            /* then multiply */
 | 
						|
            if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
 | 
						|
               goto LBL_RES;
 | 
						|
            }
 | 
						|
            if ((err = redux(&res, P, &mu)) != MP_OKAY) {
 | 
						|
               goto LBL_RES;
 | 
						|
            }
 | 
						|
         }
 | 
						|
      }
 | 
						|
   }
 | 
						|
 | 
						|
   mp_exch(&res, Y);
 | 
						|
   err = MP_OKAY;
 | 
						|
LBL_RES:
 | 
						|
   mp_clear(&res);
 | 
						|
LBL_MU:
 | 
						|
   mp_clear(&mu);
 | 
						|
LBL_M:
 | 
						|
   mp_clear(&M[1]);
 | 
						|
   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
 | 
						|
      mp_clear(&M[x]);
 | 
						|
   }
 | 
						|
   return err;
 | 
						|
}
 | 
						|
#endif
 | 
						|
 | 
						|
/* ref:         $Format:%D$ */
 | 
						|
/* git commit:  $Format:%H$ */
 | 
						|
/* commit time: $Format:%ai$ */
 |