Import of libtommath 0.35

From ltm-0.35.tar.bz2 SHA1 of 3f193dbae9351e92d02530994fa18236f7fde01c

--HG--
branch : libtommath-orig
extra : convert_revision : 2b4b13ac88b2a81e5c86ba868c92c6a452630e02

LICENSE

+LibTomMath is hereby released into the Public Domain.  
+
+-- Tom St Denis
+

TODO

+things for book in order of importance...
+
+- Fix up pseudo-code [only] for combas that are not consistent with source
+- Start in chapter 3 [basics] and work up...
+   - re-write to prose [less abrupt]
+   - clean up pseudo code [spacing]
+   - more examples where appropriate and figures
+
+Goal:
+   - Get sync done by mid January [roughly 8-12 hours work]
+   - Finish ch3-6 by end of January [roughly 12-16 hours of work]
+   - Finish ch7-end by mid Feb [roughly 20-24 hours of work].
+
+Goal isn't "first edition" but merely cleaner to read.
+
+

bn.ilg

+This is makeindex, version 2.14 [02-Oct-2002] (kpathsea + Thai support).
+Scanning input file bn.idx....done (79 entries accepted, 0 rejected).
+Sorting entries....done (511 comparisons).
+Generating output file bn.ind....done (82 lines written, 0 warnings).
+Output written in bn.ind.
+Transcript written in bn.ilg.

bn.ind

+\begin{theindex}
+
+  \item mp\_add, \hyperpage{29}
+  \item mp\_add\_d, \hyperpage{52}
+  \item mp\_and, \hyperpage{29}
+  \item mp\_clear, \hyperpage{11}
+  \item mp\_clear\_multi, \hyperpage{12}
+  \item mp\_cmp, \hyperpage{24}
+  \item mp\_cmp\_d, \hyperpage{25}
+  \item mp\_cmp\_mag, \hyperpage{23}
+  \item mp\_div, \hyperpage{30}
+  \item mp\_div\_2, \hyperpage{26}
+  \item mp\_div\_2d, \hyperpage{28}
+  \item mp\_div\_d, \hyperpage{52}
+  \item mp\_dr\_reduce, \hyperpage{40}
+  \item mp\_dr\_setup, \hyperpage{40}
+  \item MP\_EQ, \hyperpage{22}
+  \item mp\_error\_to\_string, \hyperpage{10}
+  \item mp\_expt\_d, \hyperpage{43}
+  \item mp\_exptmod, \hyperpage{43}
+  \item mp\_exteuclid, \hyperpage{51}
+  \item mp\_gcd, \hyperpage{51}
+  \item mp\_get\_int, \hyperpage{20}
+  \item mp\_grow, \hyperpage{16}
+  \item MP\_GT, \hyperpage{22}
+  \item mp\_init, \hyperpage{11}
+  \item mp\_init\_copy, \hyperpage{13}
+  \item mp\_init\_multi, \hyperpage{12}
+  \item mp\_init\_set, \hyperpage{21}
+  \item mp\_init\_set\_int, \hyperpage{21}
+  \item mp\_init\_size, \hyperpage{14}
+  \item mp\_int, \hyperpage{10}
+  \item mp\_invmod, \hyperpage{52}
+  \item mp\_jacobi, \hyperpage{52}
+  \item mp\_lcm, \hyperpage{51}
+  \item mp\_lshd, \hyperpage{28}
+  \item MP\_LT, \hyperpage{22}
+  \item MP\_MEM, \hyperpage{9}
+  \item mp\_mod, \hyperpage{35}
+  \item mp\_mod\_d, \hyperpage{52}
+  \item mp\_montgomery\_calc\_normalization, \hyperpage{38}
+  \item mp\_montgomery\_reduce, \hyperpage{37}
+  \item mp\_montgomery\_setup, \hyperpage{37}
+  \item mp\_mul, \hyperpage{31}
+  \item mp\_mul\_2, \hyperpage{26}
+  \item mp\_mul\_2d, \hyperpage{28}
+  \item mp\_mul\_d, \hyperpage{52}
+  \item mp\_n\_root, \hyperpage{44}
+  \item mp\_neg, \hyperpage{29}
+  \item MP\_NO, \hyperpage{9}
+  \item MP\_OKAY, \hyperpage{9}
+  \item mp\_or, \hyperpage{29}
+  \item mp\_prime\_fermat, \hyperpage{45}
+  \item mp\_prime\_is\_divisible, \hyperpage{45}
+  \item mp\_prime\_is\_prime, \hyperpage{46}
+  \item mp\_prime\_miller\_rabin, \hyperpage{45}
+  \item mp\_prime\_next\_prime, \hyperpage{46}
+  \item mp\_prime\_rabin\_miller\_trials, \hyperpage{46}
+  \item mp\_prime\_random, \hyperpage{47}
+  \item mp\_prime\_random\_ex, \hyperpage{47}
+  \item mp\_radix\_size, \hyperpage{49}
+  \item mp\_read\_radix, \hyperpage{49}
+  \item mp\_read\_unsigned\_bin, \hyperpage{50}
+  \item mp\_reduce, \hyperpage{36}
+  \item mp\_reduce\_2k, \hyperpage{41}
+  \item mp\_reduce\_2k\_setup, \hyperpage{41}
+  \item mp\_reduce\_setup, \hyperpage{36}
+  \item mp\_rshd, \hyperpage{28}
+  \item mp\_set, \hyperpage{19}
+  \item mp\_set\_int, \hyperpage{20}
+  \item mp\_shrink, \hyperpage{15}
+  \item mp\_sqr, \hyperpage{33}
+  \item mp\_sub, \hyperpage{29}
+  \item mp\_sub\_d, \hyperpage{52}
+  \item mp\_to\_unsigned\_bin, \hyperpage{50}
+  \item mp\_toradix, \hyperpage{49}
+  \item mp\_unsigned\_bin\_size, \hyperpage{50}
+  \item MP\_VAL, \hyperpage{9}
+  \item mp\_xor, \hyperpage{29}
+  \item MP\_YES, \hyperpage{9}
+
+\end{theindex}

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

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
+ */
+int fast_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_mod (a, b, &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

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.
+*/
+int fast_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

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.
+ *
+ */
+int fast_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 
+         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] = (mp_digit)(_W & MP_MASK);
+
+  /* setup dest */
+  olduse  = c->used;
+  c->used = pa;
+
+  {
+    register mp_digit *tmpc;
+    tmpc = c->dp;
+    for (ix = 0; ix < pa+1; 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

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.
+ */
+int fast_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 */
+      W[ix] = ((mp_digit)_W) & MP_MASK;
+
+      /* make next carry */
+      _W = _W >> ((mp_word)DIGIT_BIT);
+  }
+  
+  /* store final carry */
+  W[ix] = (mp_digit)(_W & MP_MASK);
+
+  /* setup dest */
+  olduse  = c->used;
+  c->used = pa;
+
+  {
+    register mp_digit *tmpc;
+
+    tmpc = c->dp + digs;
+    for (ix = digs; ix <= pa; 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

bn_fast_s_mp_sqr.c

+#include <tommath.h>
+#ifdef BN_FAST_S_MP_SQR_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
+ */
+
+/* the jist of squaring...
+ * you do like mult except the offset of the tmpx [one that 
+ * starts closer to zero] can't equal the offset of tmpy.  
+ * So basically you set up iy like before then you min it with
+ * (ty-tx) so that it never happens.  You double all those 
+ * you add in the inner loop
+
+After that loop you do the squares and add them in.
+*/
+
+int fast_s_mp_sqr (mp_int * a, mp_int * b)
+{
+  int       olduse, res, pa, ix, iz;
+  mp_digit   W[MP_WARRAY], *tmpx;
+  mp_word   W1;
+
+  /* grow the destination as required */
+  pa = a->used + a->used;
+  if (b->alloc < pa) {
+    if ((res = mp_grow (b, pa)) != MP_OKAY) {
+      return res;
+    }
+  }
+
+  /* number of output digits to produce */
+  W1 = 0;
+  for (ix = 0; ix < pa; ix++) { 
+      int      tx, ty, iy;
+      mp_word  _W;
+      mp_digit *tmpy;
+
+      /* clear counter */
+      _W = 0;
+
+      /* get offsets into the two bignums */
+      ty = MIN(a->used-1, ix);
+      tx = ix - ty;
+
+      /* setup temp aliases */
+      tmpx = a->dp + tx;
+      tmpy = a->dp + ty;
+
+      /* this is the number of times the loop will iterrate, essentially
+         while (tx++ < a->used && ty-- >= 0) { ... }
+       */
+      iy = MIN(a->used-tx, ty+1);
+
+      /* now for squaring tx can never equal ty 
+       * we halve the distance since they approach at a rate of 2x
+       * and we have to round because odd cases need to be executed
+       */
+      iy = MIN(iy, (ty-tx+1)>>1);
+
+      /* execute loop */
+      for (iz = 0; iz < iy; iz++) {
+         _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
+      }
+
+      /* double the inner product and add carry */
+      _W = _W + _W + W1;
+
+      /* even columns have the square term in them */
+      if ((ix&1) == 0) {
+         _W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]);
+      }
+
+      /* store it */
+      W[ix] = (mp_digit)(_W & MP_MASK);
+
+      /* make next carry */
+      W1 = _W >> ((mp_word)DIGIT_BIT);
+  }
+
+  /* setup dest */
+  olduse  = b->used;
+  b->used = a->used+a->used;
+
+  {
+    mp_digit *tmpb;
+    tmpb = b->dp;
+    for (ix = 0; ix < pa; ix++) {
+      *tmpb++ = W[ix] & MP_MASK;
+    }
+
+    /* clear unused digits [that existed in the old copy of c] */
+    for (; ix < olduse; ix++) {
+      *tmpb++ = 0;
+    }
+  }
+  mp_clamp (b);
+  return MP_OKAY;
+}
+#endif

bn_mp_2expt.c

+#include <tommath.h>
+#ifdef BN_MP_2EXPT_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 a = 2**b 
+ *
+ * Simple algorithm which zeroes the int, grows it then just sets one bit
+ * as required.
+ */
+int
+mp_2expt (mp_int * a, int b)
+{
+  int     res;
+
+  /* zero a as per default */
+  mp_zero (a);
+
+  /* grow a to accomodate the single bit */
+  if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
+    return res;
+  }
+
+  /* set the used count of where the bit will go */
+  a->used = b / DIGIT_BIT + 1;
+
+  /* put the single bit in its place */
+  a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT);
+
+  return MP_OKAY;
+}
+#endif

bn_mp_abs.c

+#include <tommath.h>
+#ifdef BN_MP_ABS_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| 
+ *
+ * Simple function copies the input and fixes the sign to positive
+ */
+int
+mp_abs (mp_int * a, mp_int * b)
+{
+  int     res;
+
+  /* copy a to b */
+  if (a != b) {
+     if ((res = mp_copy (a, b)) != MP_OKAY) {
+       return res;
+     }
+  }
+
+  /* force the sign of b to positive */
+  b->sign = MP_ZPOS;
+
+  return MP_OKAY;
+}
+#endif

bn_mp_add.c

+#include <tommath.h>
+#ifdef BN_MP_ADD_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
+ */
+
+/* high level addition (handles signs) */
+int mp_add (mp_int * a, mp_int * b, mp_int * c)
+{
+  int     sa, sb, res;
+
+  /* get sign of both inputs */
+  sa = a->sign;
+  sb = b->sign;
+
+  /* handle two cases, not four */
+  if (sa == sb) {
+    /* both positive or both negative */
+    /* add their magnitudes, copy the sign */
+    c->sign = sa;
+    res = s_mp_add (a, b, c);
+  } else {
+    /* one positive, the other negative */
+    /* subtract the one with the greater magnitude from */
+    /* the one of the lesser magnitude.  The result gets */
+    /* the sign of the one with the greater magnitude. */
+    if (mp_cmp_mag (a, b) == MP_LT) {
+      c->sign = sb;
+      res = s_mp_sub (b, a, c);
+    } else {
+      c->sign = sa;
+      res = s_mp_sub (a, b, c);
+    }
+  }
+  return res;
+}
+
+#endif

bn_mp_add_d.c

+#include <tommath.h>
+#ifdef BN_MP_ADD_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
+ */
+
+/* single digit addition */
+int
+mp_add_d (mp_int * a, mp_digit b, mp_int * c)
+{
+  int     res, ix, oldused;
+  mp_digit *tmpa, *tmpc, mu;
+
+  /* grow c as required */
+  if (c->alloc < a->used + 1) {
+     if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
+        return res;
+     }
+  }
+
+  /* if a is negative and |a| >= b, call c = |a| - b */
+  if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) {
+     /* temporarily fix sign of a */
+     a->sign = MP_ZPOS;
+
+     /* c = |a| - b */
+     res = mp_sub_d(a, b, c);
+
+     /* fix sign  */
+     a->sign = c->sign = MP_NEG;
+
+     return res;
+  }
+
+  /* old number of used digits in c */
+  oldused = c->used;
+
+  /* sign always positive */
+  c->sign = MP_ZPOS;
+
+  /* source alias */
+  tmpa    = a->dp;
+
+  /* destination alias */
+  tmpc    = c->dp;
+
+  /* if a is positive */
+  if (a->sign == MP_ZPOS) {
+     /* add digit, after this we're propagating
+      * the carry.
+      */
+     *tmpc   = *tmpa++ + b;
+     mu      = *tmpc >> DIGIT_BIT;
+     *tmpc++ &= MP_MASK;
+
+     /* now handle rest of the digits */
+     for (ix = 1; ix < a->used; ix++) {
+        *tmpc   = *tmpa++ + mu;
+        mu      = *tmpc >> DIGIT_BIT;
+        *tmpc++ &= MP_MASK;
+     }
+     /* set final carry */
+     ix++;
+     *tmpc++  = mu;
+
+     /* setup size */
+     c->used = a->used + 1;
+  } else {
+     /* a was negative and |a| < b */
+     c->used  = 1;
+
+     /* the result is a single digit */
+     if (a->used == 1) {
+        *tmpc++  =  b - a->dp[0];
+     } else {
+        *tmpc++  =  b;
+     }
+
+     /* setup count so the clearing of oldused
+      * can fall through correctly
+      */
+     ix       = 1;
+  }
+
+  /* now zero to oldused */
+  while (ix++ < oldused) {
+     *tmpc++ = 0;
+  }
+  mp_clamp(c);
+
+  return MP_OKAY;
+}
+
+#endif

bn_mp_addmod.c

+#include <tommath.h>
+#ifdef BN_MP_ADDMOD_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
+ */
+
+/* d = a + b (mod c) */
+int
+mp_addmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
+{
+  int     res;
+  mp_int  t;
+
+  if ((res = mp_init (&t)) != MP_OKAY) {
+    return res;
+  }
+
+  if ((res = mp_add (a, b, &t)) != MP_OKAY) {
+    mp_clear (&t);
+    return res;
+  }
+  res = mp_mod (&t, c, d);
+  mp_clear (&t);
+  return res;
+}
+#endif

bn_mp_and.c

+#include <tommath.h>
+#ifdef BN_MP_AND_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
+ */
+
+/* AND two ints together */
+int
+mp_and (mp_int * a, mp_int * b, mp_int * c)
+{
+  int     res, ix, px;
+  mp_int  t, *x;
+
+  if (a->used > b->used) {
+    if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
+      return res;
+    }
+    px = b->used;
+    x = b;
+  } else {
+    if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
+      return res;
+    }
+    px = a->used;
+    x = a;
+  }
+
+  for (ix = 0; ix < px; ix++) {
+    t.dp[ix] &= x->dp[ix];
+  }
+
+  /* zero digits above the last from the smallest mp_int */
+  for (; ix < t.used; ix++) {
+    t.dp[ix] = 0;
+  }
+
+  mp_clamp (&t);
+  mp_exch (c, &t);
+  mp_clear (&t);
+  return MP_OKAY;
+}
+#endif

bn_mp_clamp.c

+#include <tommath.h>
+#ifdef BN_MP_CLAMP_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
+ */
+
+/* trim unused digits 
+ *
+ * This is used to ensure that leading zero digits are
+ * trimed and the leading "used" digit will be non-zero
+ * Typically very fast.  Also fixes the sign if there
+ * are no more leading digits
+ */
+void
+mp_clamp (mp_int * a)
+{
+  /* decrease used while the most significant digit is
+   * zero.
+   */
+  while (a->used > 0 && a->dp[a->used - 1] == 0) {
+    --(a->used);
+  }
+
+  /* reset the sign flag if used == 0 */
+  if (a->used == 0) {
+    a->sign = MP_ZPOS;
+  }
+}
+#endif

bn_mp_clear.c

+#include <tommath.h>
+#ifdef BN_MP_CLEAR_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
+ */
+
+/* clear one (frees)  */
+void
+mp_clear (mp_int * a)
+{
+  int i;
+
+  /* only do anything if a hasn't been freed previously */
+  if (a->dp != NULL) {
+    /* first zero the digits */
+    for (i = 0; i < a->used; i++) {
+        a->dp[i] = 0;
+    }
+
+    /* free ram */
+    XFREE(a->dp);
+
+    /* reset members to make debugging easier */
+    a->dp    = NULL;
+    a->alloc = a->used = 0;
+    a->sign  = MP_ZPOS;
+  }
+}
+#endif