/* fix the sign */  neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;  x.sign = y.sign = MP_ZPOS;  /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */  norm = mp_count_bits(&y) % DIGIT_BIT;  if (norm < (int)(DIGIT_BIT-1)) {     norm = (DIGIT_BIT-1) - norm;     if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {       goto LBL_Y;     }     if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {       goto LBL_Y;     }  } else {     norm = 0;  }  /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */  n = x.used - 1;  t = y.used - 1;  /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */  if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */    goto LBL_Y;  }  while (mp_cmp (&x, &y) != MP_LT) {    ++(q.dp[n - t]);    if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {      goto LBL_Y;    }  }  /* reset y by shifting it back down */  mp_rshd (&y, n - t);  /* step 3. for i from n down to (t + 1) */  for (i = n; i >= (t + 1); i--) {    if (i > x.used) {      continue;    }    /* step 3.1 if xi == yt then set q{i-t-1} to b-1,      * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */    if (x.dp[i] == y.dp[t]) {      q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);    } else {      mp_word tmp;      tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);      tmp |= ((mp_word) x.dp[i - 1]);      tmp /= ((mp_word) y.dp[t]);      if (tmp > (mp_word) MP_MASK)        tmp = MP_MASK;      q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));    }    /* while (q{i-t-1} * (yt * b + y{t-1})) >              xi * b**2 + xi-1 * b + xi-2             do q{i-t-1} -= 1;     */    q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;    do {      q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;      /* find left hand */      mp_zero (&t1);      t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];      t1.dp[1] = y.dp[t];      t1.used = 2;      if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {        goto LBL_Y;      }      /* find right hand */      t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];      t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];      t2.dp[2] = x.dp[i];      t2.used = 3;    } while (mp_cmp_mag(&t1, &t2) == MP_GT);    /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */    if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {      goto LBL_Y;    }    if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {      goto LBL_Y;    }    if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {      goto LBL_Y;    }    /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */    if (x.sign == MP_NEG) {      if ((res = mp_copy (&y, &t1)) != MP_OKAY) {        goto LBL_Y;      }      if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {        goto LBL_Y;      }      if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {        goto LBL_Y;      }      q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;    }  }  /* now q is the quotient and x is the remainder    * [which we have to normalize]    */    /* get sign before writing to c */  x.sign = x.used == 0 ? MP_ZPOS : a->sign;  if (c != NULL) {    mp_clamp (&q);    mp_exch (&q, c);    c->sign = neg;  }  if (d != NULL) {    mp_div_2d (&x, norm, &x, NULL);    mp_exch (&x, d);  }  res = MP_OKAY;LBL_Y:mp_clear (&y);LBL_X:mp_clear (&x);LBL_T2:mp_clear (&t2);LBL_T1:mp_clear (&t1);LBL_Q:mp_clear (&q);  return res;}#endif#endif/* End: bn_mp_div.c *//* Start: bn_mp_div_2.c */#include <tommath.h>#ifdef BN_MP_DIV_2_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 *//* b = a/2 */int mp_div_2(mp_int * a, mp_int * b){  int     x, res, oldused;  /* copy */  if (b->alloc < a->used) {    if ((res = mp_grow (b, a->used)) != MP_OKAY) {      return res;    }  }  oldused = b->used;  b->used = a->used;  {    register mp_digit r, rr, *tmpa, *tmpb;    /* source alias */    tmpa = a->dp + b->used - 1;    /* dest alias */    tmpb = b->dp + b->used - 1;    /* carry */    r = 0;    for (x = b->used - 1; x >= 0; x--) {      /* get the carry for the next iteration */      rr = *tmpa & 1;      /* shift the current digit, add in carry and store */      *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));      /* forward carry to next iteration */      r = rr;    }    /* zero excess digits */    tmpb = b->dp + b->used;    for (x = b->used; x < oldused; x++) {      *tmpb++ = 0;    }  }  b->sign = a->sign;  mp_clamp (b);  return MP_OKAY;}#endif/* End: bn_mp_div_2.c *//* Start: bn_mp_div_2d.c */#include <tommath.h>#ifdef BN_MP_DIV_2D_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 *//* shift right by a certain bit count (store quotient in c, optional remainder in d) */int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d){  mp_digit D, r, rr;  int     x, res;  mp_int  t;  /* if the shift count is <= 0 then we do no work */  if (b <= 0) {    res = mp_copy (a, c);    if (d != NULL) {      mp_zero (d);    }    return res;  }  if ((res = mp_init (&t)) != MP_OKAY) {    return res;  }  /* get the remainder */  if (d != NULL) {    if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {      mp_clear (&t);      return res;    }  }  /* copy */  if ((res = mp_copy (a, c)) != MP_OKAY) {    mp_clear (&t);    return res;  }  /* shift by as many digits in the bit count */  if (b >= (int)DIGIT_BIT) {    mp_rshd (c, b / DIGIT_BIT);  }  /* shift any bit count < DIGIT_BIT */  D = (mp_digit) (b % DIGIT_BIT);  if (D != 0) {    register mp_digit *tmpc, mask, shift;    /* mask */    mask = (((mp_digit)1) << D) - 1;    /* shift for lsb */    shift = DIGIT_BIT - D;    /* alias */    tmpc = c->dp + (c->used - 1);    /* carry */    r = 0;    for (x = c->used - 1; x >= 0; x--) {      /* get the lower  bits of this word in a temp */      rr = *tmpc & mask;      /* shift the current word and mix in the carry bits from the previous word */      *tmpc = (*tmpc >> D) | (r << shift);      --tmpc;      /* set the carry to the carry bits of the current word found above */      r = rr;    }  }  mp_clamp (c);  if (d != NULL) {    mp_exch (&t, d);  }  mp_clear (&t);  return MP_OKAY;}#endif/* End: bn_mp_div_2d.c *//* Start: bn_mp_div_3.c */#include <tommath.h>#ifdef BN_MP_DIV_3_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 *//* divide by three (based on routine from MPI and the GMP manual) */intmp_div_3 (mp_int * a, mp_int *c, mp_digit * d){  mp_int   q;  mp_word  w, t;  mp_digit b;  int      res, ix;    /* b = 2**DIGIT_BIT / 3 */  b = (((mp_word)1) << ((mp_word)DIGIT_BIT)) / ((mp_word)3);  if ((res = mp_init_size(&q, a->used)) != MP_OKAY) {     return res;  }    q.used = a->used;  q.sign = a->sign;  w = 0;  for (ix = a->used - 1; ix >= 0; ix--) {     w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);     if (w >= 3) {        /* multiply w by [1/3] */        t = (w * ((mp_word)b)) >> ((mp_word)DIGIT_BIT);        /* now subtract 3 * [w/3] from w, to get the remainder */        w -= t+t+t;        /* fixup the remainder as required since         * the optimization is not exact.         */        while (w >= 3) {           t += 1;           w -= 3;        }      } else {        t = 0;      }      q.dp[ix] = (mp_digit)t;  }  /* [optional] store the remainder */  if (d != NULL) {     *d = (mp_digit)w;  }  /* [optional] store the quotient */  if (c != NULL) {     mp_clamp(&q);     mp_exch(&q, c);  }  mp_clear(&q);    return res;}#endif/* End: bn_mp_div_3.c *//* Start: bn_mp_div_d.c */#include <tommath.h>#ifdef BN_MP_DIV_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 */static int s_is_power_of_two(mp_digit b, int *p){   int x;   for (x = 1; x < DIGIT_BIT; x++) {      if (b == (((mp_digit)1)<<x)) {         *p = x;         return 1;      }   }   return 0;}/* single digit division (based on routine from MPI) */int mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d){  mp_int  q;  mp_word w;  mp_digit t;  int     res, ix;  /* cannot divide by zero */  if (b == 0) {     return MP_VAL;  }  /* quick outs */  if (b == 1 || mp_iszero(a) == 1) {     if (d != NULL) {        *d = 0;     }     if (c != NULL) {        return mp_copy(a, c);     }     return MP_OKAY;  }  /* power of two ? */  if (s_is_power_of_two(b, &ix) == 1) {     if (d != NULL) {        *d = a->dp[0] & ((((mp_digit)1)<<ix) - 1);     }     if (c != NULL) {        return mp_div_2d(a, ix, c, NULL);     }     return MP_OKAY;  }#ifdef BN_MP_DIV_3_C  /* three? */  if (b == 3) {     return mp_div_3(a, c, d);  }#endif  /* no easy answer [c'est la vie].  Just division */  if ((res = mp_init_size(&q, a->used)) != MP_OKAY) {     return res;  }    q.used = a->used;  q.sign = a->sign;  w = 0;  for (ix = a->used - 1; ix >= 0; ix--) {     w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);          if (w >= b) {        t = (mp_digit)(w / b);        w -= ((mp_word)t) * ((mp_word)b);      } else {        t = 0;      }      q.dp[ix] = (mp_digit)t;  }    if (d != NULL) {     *d = (mp_digit)w;  }    if (c != NULL) {     mp_clamp(&q);     mp_exch(&q, c);  }  mp_clear(&q);    return res;}#endif/* End: bn_mp_div_d.c *//* Start: bn_mp_dr_is_modulus.c */#include <tommath.h>#ifdef BN_MP_DR_IS_MODULUS_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 *//* determines if a number is a valid DR modulus */int mp_dr_is_modulus(mp_int *a){   int ix;   /* must be at least two digits */   if (a->used < 2) {      return 0;   }   /* must be of the form b**k - a [a <= b] so all    * but the first digit must be equal to -1 (mod b).    */   for (ix = 1; ix < a->used; ix++) {       if (a->dp[ix] != MP_MASK) {          return 0;       }   }   return 1;}#endif/* End: bn_mp_dr_is_modulus.c *//* Start: bn_mp_dr_reduce.c */#include <tommath.h>#ifdef BN_MP_DR_REDUCE_C/* LibTomMath, multiple-precision integer library -- Tom St Denis