/* Start: bn_error.c */#include <tommath.h>#ifdef BN_ERROR_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 */static const struct {     int code;     char *msg;} msgs[] = {     { MP_OKAY, "Successful" },     { MP_MEM,  "Out of heap" },     { MP_VAL,  "Value out of range" }};/* return a char * string for a given code */char *mp_error_to_string(int code){   int x;   /* scan the lookup table for the given message */   for (x = 0; x < (int)(sizeof(msgs) / sizeof(msgs[0])); x++) {       if (msgs[x].code == code) {          return msgs[x].msg;       }   }   /* generic reply for invalid code */   return "Invalid error code";}#endif/* End: bn_error.c *//* Start: bn_fast_mp_invmod.c */#include <tommath.h>#ifdef BN_FAST_MP_INVMOD_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 *//* computes the modular inverse via binary extended euclidean algorithm,  * that is c = 1/a mod b  * * Based on slow invmod except this is optimized for the case where b is  * odd as per HAC Note 14.64 on pp. 610 */intfast_mp_invmod (mp_int * a, mp_int * b, mp_int * c){  mp_int  x, y, u, v, B, D;  int     res, neg;  /* 2. [modified] b must be odd   */  if (mp_iseven (b) == 1) {    return MP_VAL;  }  /* init all our temps */  if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {     return res;  }  /* x == modulus, y == value to invert */  if ((res = mp_copy (b, &x)) != MP_OKAY) {    goto LBL_ERR;  }  /* we need y = |a| */  if ((res = mp_abs (a, &y)) != MP_OKAY) {    goto LBL_ERR;  }  /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */  if ((res = mp_copy (&x, &u)) != MP_OKAY) {    goto LBL_ERR;  }  if ((res = mp_copy (&y, &v)) != MP_OKAY) {    goto LBL_ERR;  }  mp_set (&D, 1);top:  /* 4.  while u is even do */  while (mp_iseven (&u) == 1) {    /* 4.1 u = u/2 */    if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {      goto LBL_ERR;    }    /* 4.2 if B is odd then */    if (mp_isodd (&B) == 1) {      if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {        goto LBL_ERR;      }    }    /* B = B/2 */    if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {      goto LBL_ERR;    }  }  /* 5.  while v is even do */  while (mp_iseven (&v) == 1) {    /* 5.1 v = v/2 */    if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {      goto LBL_ERR;    }    /* 5.2 if D is odd then */    if (mp_isodd (&D) == 1) {      /* D = (D-x)/2 */      if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {        goto LBL_ERR;      }    }    /* D = D/2 */    if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {      goto LBL_ERR;    }  }  /* 6.  if u >= v then */  if (mp_cmp (&u, &v) != MP_LT) {    /* u = u - v, B = B - D */    if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {      goto LBL_ERR;    }    if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {      goto LBL_ERR;    }  } else {    /* v - v - u, D = D - B */    if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {      goto LBL_ERR;    }    if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {      goto LBL_ERR;    }  }  /* if not zero goto step 4 */  if (mp_iszero (&u) == 0) {    goto top;  }  /* now a = C, b = D, gcd == g*v */  /* if v != 1 then there is no inverse */  if (mp_cmp_d (&v, 1) != MP_EQ) {    res = MP_VAL;    goto LBL_ERR;  }  /* b is now the inverse */  neg = a->sign;  while (D.sign == MP_NEG) {    if ((res = mp_add (&D, b, &D)) != MP_OKAY) {      goto LBL_ERR;    }  }  mp_exch (&D, c);  c->sign = neg;  res = MP_OKAY;LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);  return res;}#endif/* End: bn_fast_mp_invmod.c *//* Start: bn_fast_mp_montgomery_reduce.c */#include <tommath.h>#ifdef BN_FAST_MP_MONTGOMERY_REDUCE_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 *//* computes xR**-1 == x (mod N) via Montgomery Reduction * * This is an optimized implementation of montgomery_reduce * which uses the comba method to quickly calculate the columns of the * reduction. * * Based on Algorithm 14.32 on pp.601 of HAC.*/intfast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho){  int     ix, res, olduse;  mp_word W[MP_WARRAY];  /* get old used count */  olduse = x->used;  /* grow a as required */  if (x->alloc < n->used + 1) {    if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) {      return res;    }  }  /* first we have to get the digits of the input into   * an array of double precision words W[...]   */  {    register mp_word *_W;    register mp_digit *tmpx;    /* alias for the W[] array */    _W   = W;    /* alias for the digits of  x*/    tmpx = x->dp;    /* copy the digits of a into W[0..a->used-1] */    for (ix = 0; ix < x->used; ix++) {      *_W++ = *tmpx++;    }    /* zero the high words of W[a->used..m->used*2] */    for (; ix < n->used * 2 + 1; ix++) {      *_W++ = 0;    }  }  /* now we proceed to zero successive digits   * from the least significant upwards   */  for (ix = 0; ix < n->used; ix++) {    /* mu = ai * m' mod b     *     * We avoid a double precision multiplication (which isn't required)     * by casting the value down to a mp_digit.  Note this requires     * that W[ix-1] have  the carry cleared (see after the inner loop)     */    register mp_digit mu;    mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);    /* a = a + mu * m * b**i     *     * This is computed in place and on the fly.  The multiplication     * by b**i is handled by offseting which columns the results     * are added to.     *     * Note the comba method normally doesn't handle carries in the     * inner loop In this case we fix the carry from the previous     * column since the Montgomery reduction requires digits of the     * result (so far) [see above] to work.  This is     * handled by fixing up one carry after the inner loop.  The     * carry fixups are done in order so after these loops the     * first m->used words of W[] have the carries fixed     */    {      register int iy;      register mp_digit *tmpn;      register mp_word *_W;      /* alias for the digits of the modulus */      tmpn = n->dp;      /* Alias for the columns set by an offset of ix */      _W = W + ix;      /* inner loop */      for (iy = 0; iy < n->used; iy++) {          *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);      }    }    /* now fix carry for next digit, W[ix+1] */    W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);  }  /* now we have to propagate the carries and   * shift the words downward [all those least   * significant digits we zeroed].   */  {    register mp_digit *tmpx;    register mp_word *_W, *_W1;    /* nox fix rest of carries */    /* alias for current word */    _W1 = W + ix;    /* alias for next word, where the carry goes */    _W = W + ++ix;    for (; ix <= n->used * 2 + 1; ix++) {      *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);    }    /* copy out, A = A/b**n     *     * The result is A/b**n but instead of converting from an     * array of mp_word to mp_digit than calling mp_rshd     * we just copy them in the right order     */    /* alias for destination word */    tmpx = x->dp;    /* alias for shifted double precision result */    _W = W + n->used;    for (ix = 0; ix < n->used + 1; ix++) {      *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));    }    /* zero oldused digits, if the input a was larger than     * m->used+1 we'll have to clear the digits     */    for (; ix < olduse; ix++) {      *tmpx++ = 0;    }  }  /* set the max used and clamp */  x->used = n->used + 1;  mp_clamp (x);  /* if A >= m then A = A - m */  if (mp_cmp_mag (x, n) != MP_LT) {    return s_mp_sub (x, n, x);  }  return MP_OKAY;}#endif/* End: bn_fast_mp_montgomery_reduce.c *//* Start: bn_fast_s_mp_mul_digs.c */#include <tommath.h>#ifdef BN_FAST_S_MP_MUL_DIGS_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 *//* Fast (comba) multiplier * * This is the fast column-array [comba] multiplier.  It is  * designed to compute the columns of the product first  * then handle the carries afterwards.  This has the effect  * of making the nested loops that compute the columns very * simple and schedulable on super-scalar processors. * * This has been modified to produce a variable number of  * digits of output so if say only a half-product is required  * you don't have to compute the upper half (a feature  * required for fast Barrett reduction). * * Based on Algorithm 14.12 on pp.595 of HAC. * */intfast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs){  int     olduse, res, pa, ix, iz;  mp_digit W[MP_WARRAY];  register mp_word  _W;  /* grow the destination as required */  if (c->alloc < digs) {    if ((res = mp_grow (c, digs)) != MP_OKAY) {      return res;    }  }  /* number of output digits to produce */  pa = MIN(digs, a->used + b->used);  /* clear the carry */  _W = 0;  for (ix = 0; ix < pa; ix++) {       int      tx, ty;      int      iy;      mp_digit *tmpx, *tmpy;      /* get offsets into the two bignums */      ty = MIN(b->used-1, ix);      tx = ix - ty;      /* setup temp aliases */      tmpx = a->dp + tx;      tmpy = b->dp + ty;      /* this is the number of times the loop will iterrate, essentially its          while (tx++ < a->used && ty-- >= 0) { ... }       */      iy = MIN(a->used-tx, ty+1);      /* execute loop */      for (iz = 0; iz < iy; ++iz) {         _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);      }      /* store term */      W[ix] = ((mp_digit)_W) & MP_MASK;      /* make next carry */      _W = _W >> ((mp_word)DIGIT_BIT);  }  /* store final carry */  W[ix] = _W;  /* setup dest */  olduse  = c->used;  c->used = digs;  {    register mp_digit *tmpc;    tmpc = c->dp;    for (ix = 0; ix < digs; ix++) {      /* now extract the previous digit [below the carry] */      *tmpc++ = W[ix];    }    /* clear unused digits [that existed in the old copy of c] */    for (; ix < olduse; ix++) {      *tmpc++ = 0;    }  }  mp_clamp (c);  return MP_OKAY;}#endif/* End: bn_fast_s_mp_mul_digs.c *//* Start: bn_fast_s_mp_mul_high_digs.c */#include <tommath.h>#ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_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 *//* this is a modified version of fast_s_mul_digs that only produces * output digits *above* digs.  See the comments for fast_s_mul_digs * to see how it works. * * This is used in the Barrett reduction since for one of the multiplications * only the higher digits were needed.  This essentially halves the work. * * Based on Algorithm 14.12 on pp.595 of HAC. */intfast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs){  int     olduse, res, pa, ix, iz;  mp_digit W[MP_WARRAY];  mp_word  _W;  /* grow the destination as required */  pa = a->used + b->used;  if (c->alloc < pa) {    if ((res = mp_grow (c, pa)) != MP_OKAY) {      return res;    }  }  /* number of output digits to produce */  pa = a->used + b->used;  _W = 0;  for (ix = digs; ix < pa; ix++) {       int      tx, ty, iy;      mp_digit *tmpx, *tmpy;      /* get offsets into the two bignums */      ty = MIN(b->used-1, ix);      tx = ix - ty;      /* setup temp aliases */      tmpx = a->dp + tx;      tmpy = b->dp + ty;      /* this is the number of times the loop will iterrate, essentially its          while (tx++ < a->used && ty-- >= 0) { ... }       */      iy = MIN(a->used-tx, ty+1);      /* execute loop */      for (iz = 0; iz < iy; iz++) {         _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);      }      /* store term */