* guarantee it works. * * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org */#include <ltc_tommath.h>/* reduce "x" in place modulo "n" using the Diminished Radix algorithm. * * Based on algorithm from the paper * * "Generating Efficient Primes for Discrete Log Cryptosystems" *                 Chae Hoon Lim, Pil Loong Lee, *          POSTECH Information Research Laboratories * * The modulus must be of a special format [see manual] * * Has been modified to use algorithm 7.10 from the LTM book instead * * Input x must be in the range 0 <= x <= (n-1)**2 */intmp_dr_reduce (mp_int * x, mp_int * n, mp_digit k){  int      err, i, m;  mp_word  r;  mp_digit mu, *tmpx1, *tmpx2;  /* m = digits in modulus */  m = n->used;  /* ensure that "x" has at least 2m digits */  if (x->alloc < m + m) {    if ((err = mp_grow (x, m + m)) != MP_OKAY) {      return err;    }  }/* top of loop, this is where the code resumes if * another reduction pass is required. */top:  /* aliases for digits */  /* alias for lower half of x */  tmpx1 = x->dp;  /* alias for upper half of x, or x/B**m */  tmpx2 = x->dp + m;  /* set carry to zero */  mu = 0;  /* compute (x mod B**m) + k * [x/B**m] inline and inplace */  for (i = 0; i < m; i++) {      r         = ((mp_word)*tmpx2++) * ((mp_word)k) + *tmpx1 + mu;      *tmpx1++  = (mp_digit)(r & MP_MASK);      mu        = (mp_digit)(r >> ((mp_word)DIGIT_BIT));  }  /* set final carry */  *tmpx1++ = mu;  /* zero words above m */  for (i = m + 1; i < x->used; i++) {      *tmpx1++ = 0;  }  /* clamp, sub and return */  mp_clamp (x);  /* if x >= n then subtract and reduce again   * Each successive "recursion" makes the input smaller and smaller.   */  if (mp_cmp_mag (x, n) != MP_LT) {    s_mp_sub(x, n, x);    goto top;  }  return MP_OKAY;}/* End: bn_mp_dr_reduce.c *//* Start: bn_mp_dr_setup.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, tomstdenis@iahu.ca, http://math.libtomcrypt.org */#include <ltc_tommath.h>/* determines the setup value */void mp_dr_setup(mp_int *a, mp_digit *d){   /* the casts are required if DIGIT_BIT is one less than    * the number of bits in a mp_digit [e.g. DIGIT_BIT==31]    */   *d = (mp_digit)((((mp_word)1) << ((mp_word)DIGIT_BIT)) -         ((mp_word)a->dp[0]));}/* End: bn_mp_dr_setup.c *//* Start: bn_mp_exch.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, tomstdenis@iahu.ca, http://math.libtomcrypt.org */#include <ltc_tommath.h>/* swap the elements of two integers, for cases where you can't simply swap the  * mp_int pointers around */voidmp_exch (mp_int * a, mp_int * b){  mp_int  t;  t  = *a;  *a = *b;  *b = t;}/* End: bn_mp_exch.c *//* Start: bn_mp_expt_d.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, tomstdenis@iahu.ca, http://math.libtomcrypt.org */#include <ltc_tommath.h>/* calculate c = a**b  using a square-multiply algorithm */int mp_expt_d (mp_int * a, mp_digit b, mp_int * c){  int     res, x;  mp_int  g;  if ((res = mp_init_copy (&g, a)) != MP_OKAY) {    return res;  }  /* set initial result */  mp_set (c, 1);  for (x = 0; x < (int) DIGIT_BIT; x++) {    /* square */    if ((res = mp_sqr (c, c)) != MP_OKAY) {      mp_clear (&g);      return res;    }    /* if the bit is set multiply */    if ((b & (mp_digit) (((mp_digit)1) << (DIGIT_BIT - 1))) != 0) {      if ((res = mp_mul (c, &g, c)) != MP_OKAY) {         mp_clear (&g);         return res;      }    }    /* shift to next bit */    b <<= 1;  }  mp_clear (&g);  return MP_OKAY;}/* End: bn_mp_expt_d.c *//* Start: bn_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, tomstdenis@iahu.ca, http://math.libtomcrypt.org */#include <ltc_tommath.h>/* this is a shell function that calls either the normal or Montgomery * exptmod functions.  Originally the call to the montgomery code was * embedded in the normal function but that wasted alot of stack space * for nothing (since 99% of the time the Montgomery code would be called) */int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y){  int dr;  /* modulus P must be positive */  if (P->sign == MP_NEG) {     return MP_VAL;  }  /* if exponent X is negative we have to recurse */  if (X->sign == MP_NEG) {     mp_int tmpG, tmpX;     int err;     /* first compute 1/G mod P */     if ((err = mp_init(&tmpG)) != MP_OKAY) {        return err;     }     if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {        mp_clear(&tmpG);        return err;     }     /* now get |X| */     if ((err = mp_init(&tmpX)) != MP_OKAY) {        mp_clear(&tmpG);        return err;     }     if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {        mp_clear_multi(&tmpG, &tmpX, NULL);        return err;     }     /* and now compute (1/G)**|X| instead of G**X [X < 0] */     err = mp_exptmod(&tmpG, &tmpX, P, Y);     mp_clear_multi(&tmpG, &tmpX, NULL);     return err;  }  /* is it a DR modulus? */  dr = mp_dr_is_modulus(P);  /* if not, is it a uDR modulus? */  if (dr == 0) {     dr = mp_reduce_is_2k(P) << 1;  }      /* if the modulus is odd or dr != 0 use the fast method */  if (mp_isodd (P) == 1 || dr !=  0) {    return mp_exptmod_fast (G, X, P, Y, dr);  } else {    /* otherwise use the generic Barrett reduction technique */    return s_mp_exptmod (G, X, P, Y);  }}/* End: bn_mp_exptmod.c *//* Start: bn_mp_exptmod_fast.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, tomstdenis@iahu.ca, http://math.libtomcrypt.org */#include <ltc_tommath.h>/* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85 * * Uses a left-to-right k-ary sliding window to compute the modular exponentiation. * The value of k changes based on the size of the exponent. * * Uses Montgomery or Diminished Radix reduction [whichever appropriate] */#ifdef MP_LOW_MEM   #define TAB_SIZE 32#else   #define TAB_SIZE 256#endifintmp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode){  mp_int  M[TAB_SIZE], res;  mp_digit buf, mp;  int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;  /* use a pointer to the reduction algorithm.  This allows us to use   * one of many reduction algorithms without modding the guts of   * the code with if statements everywhere.   */  int     (*redux)(mp_int*,mp_int*,mp_digit);  /* 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;    }  }  /* determine and setup reduction code */  if (redmode == 0) {     /* now setup montgomery  */     if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {        goto __M;     }     /* automatically pick the comba one if available (saves quite a few calls/ifs) */     if (((P->used * 2 + 1) < MP_WARRAY) &&          P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {        redux = fast_mp_montgomery_reduce;     } else {        /* use slower baseline Montgomery method */        redux = mp_montgomery_reduce;     }  } else if (redmode == 1) {     /* setup DR reduction for moduli of the form B**k - b */     mp_dr_setup(P, &mp);     redux = mp_dr_reduce;  } else {     /* setup DR reduction for moduli of the form 2**k - b */     if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {        goto __M;     }     redux = mp_reduce_2k;  }  /* setup result */  if ((err = mp_init (&res)) != MP_OKAY) {    goto __M;  }  /* create M table   *   * The M table contains powers of the input 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 (redmode == 0) {     /* now we need R mod m */     if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {       goto __RES;     }     /* now set M[1] to G * R mod m */     if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {       goto __RES;     }  } else {     mp_set(&res, 1);     if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {        goto __RES;     }  }  /* 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 __RES;  }  for (x = 0; x < (winsize - 1); x++) {    if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {      goto __RES;    }    if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {      goto __RES;    }  }  /* create upper table */  for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {    if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {      goto __RES;    }    if ((err = redux (&M[x], P, mp)) != MP_OKAY) {      goto __RES;    }  }  /* 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 so break */      if (digidx == -1) {        break;      }      /* read next digit and reset bitcnt */      buf    = X->dp[digidx--];      bitcnt = (int)DIGIT_BIT;    }    /* grab the next msb from the exponent */    y     = (mp_digit)(buf >> (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 __RES;      }      if ((err = redux (&res, P, mp)) != MP_OKAY) {        goto __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 __RES;        }        if ((err = redux (&res, P, mp)) != MP_OKAY) {          goto __RES;        }      }      /* then multiply */      if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {        goto __RES;      }      if ((err = redux (&res, P, mp)) != MP_OKAY) {        goto __RES;      }      /* empty window and reset */      bitcpy = 0;      bitbuf = 0;      mode   = 1;    }  }  /* if bits remain then square/multiply */  if (mode == 2 && bitcpy > 0) {