added libtommath-0.34

bn.tex

@@ -49,7 +49,7 @@
 \begin{document}
 \frontmatter
 \pagestyle{empty}
-\title{LibTomMath User Manual \\ v0.33}
+\title{LibTomMath User Manual \\ v0.34}
 \author{Tom St Denis \\ tomstdenis@iahu.ca}
 \maketitle
 This text, the library and the accompanying textbook are all hereby placed in the public domain.  This book has been 
@@ -263,12 +263,12 @@ are the pros and cons of LibTomMath by comparing it to the math routines from Gn
 \begin{center}
 \begin{tabular}{|l|c|c|l|}
 \hline \textbf{Criteria} & \textbf{Pro} & \textbf{Con} & \textbf{Notes} \\
-\hline Few lines of code per file & X & & GnuPG $ = 300.9$, LibTomMath  $ = 76.04$ \\
+\hline Few lines of code per file & X & & GnuPG $ = 300.9$, LibTomMath  $ = 71.97$ \\
 \hline Commented function prototypes & X && GnuPG function names are cryptic. \\
 \hline Speed && X & LibTomMath is slower.  \\
 \hline Totally free & X & & GPL has unfavourable restrictions.\\
 \hline Large function base & X & & GnuPG is barebones. \\
-\hline Four modular reduction algorithms & X & & Faster modular exponentiation. \\
+\hline Five modular reduction algorithms & X & & Faster modular exponentiation for a variety of moduli. \\
 \hline Portable & X & & GnuPG requires configuration to build. \\
 \hline
 \end{tabular}
@@ -284,9 +284,12 @@ would require when working with large integers.
 So it may feel tempting to just rip the math code out of GnuPG (or GnuMP where it was taken from originally) in your
 own application but I think there are reasons not to.  While LibTomMath is slower than libraries such as GnuMP it is
 not normally significantly slower.  On x86 machines the difference is normally a factor of two when performing modular
-exponentiations.
+exponentiations.  It depends largely on the processor, compiler and the moduli being used.
 
-Essentially the only time you wouldn't use LibTomMath is when blazing speed is the primary concern.
+Essentially the only time you wouldn't use LibTomMath is when blazing speed is the primary concern.  However,
+on the other side of the coin LibTomMath offers you a totally free (public domain) well structured math library
+that is very flexible, complete and performs well in resource contrained environments.  Fast RSA for example can
+be performed with as little as 8KB of ram for data (again depending on build options).  
 
 \chapter{Getting Started with LibTomMath}
 \section{Building Programs}
@@ -809,7 +812,7 @@ mp\_int variables based on their digits only.
 
 \index{mp\_cmp\_mag}
 \begin{alltt}
-int mp_cmp(mp_int * a, mp_int * b);
+int mp_cmp_mag(mp_int * a, mp_int * b);
 \end{alltt}
 This will compare $a$ to $b$ placing $a$ to the left of $b$.  This function cannot fail and will return one of the
 three compare codes listed in figure \ref{fig:CMP}.
@@ -1220,12 +1223,13 @@ int mp_sqr (mp_int * a, mp_int * b);
 \end{alltt}
 
 Will square $a$ and store it in $b$.  Like the case of multiplication there are four different squaring
-algorithms all which can be called from mp\_sqr().  It is ideal to use mp\_sqr over mp\_mul when squaring terms.
+algorithms all which can be called from mp\_sqr().  It is ideal to use mp\_sqr over mp\_mul when squaring terms because
+of the speed difference.  
 
 \section{Tuning Polynomial Basis Routines}
 
 Both of the Toom-Cook and Karatsuba multiplication algorithms are faster than the traditional $O(n^2)$ approach that
-the Comba and baseline algorithms use.  At $O(n^{1.464973})$ and $O(n^{1.584962})$ running times respectfully they require 
+the Comba and baseline algorithms use.  At $O(n^{1.464973})$ and $O(n^{1.584962})$ running times respectively they require 
 considerably less work.  For example, a 10000-digit multiplication would take roughly 724,000 single precision
 multiplications with Toom-Cook or 100,000,000 single precision multiplications with the standard Comba (a factor
 of 138).
@@ -1297,14 +1301,14 @@ of $b$.  This algorithm accepts an input $a$ of any range and is not limited by
 \section{Barrett Reduction}
 
 Barrett reduction is a generic optimized reduction algorithm that requires pre--computation to achieve
-a decent speedup over straight division.  First a $mu$ value must be precomputed with the following function.
+a decent speedup over straight division.  First a $\mu$ value must be precomputed with the following function.
 
 \index{mp\_reduce\_setup}
 \begin{alltt}
 int mp_reduce_setup(mp_int *a, mp_int *b);
 \end{alltt}
 
-Given a modulus in $b$ this produces the required $mu$ value in $a$.  For any given modulus this only has to
+Given a modulus in $b$ this produces the required $\mu$ value in $a$.  For any given modulus this only has to
 be computed once.  Modular reduction can now be performed with the following.
 
 \index{mp\_reduce}
@@ -1312,7 +1316,7 @@ be computed once.  Modular reduction can now be performed with the following.
 int mp_reduce(mp_int *a, mp_int *b, mp_int *c);
 \end{alltt}
 
-This will reduce $a$ in place modulo $b$ with the precomputed $mu$ value in $c$.  $a$ must be in the range
+This will reduce $a$ in place modulo $b$ with the precomputed $\mu$ value in $c$.  $a$ must be in the range
 $0 \le a < b^2$.
 
 \begin{alltt}
@@ -1578,7 +1582,8 @@ will return $-2$.
 This algorithm uses the ``Newton Approximation'' method and will converge on the correct root fairly quickly.  Since
 the algorithm requires raising $a$ to the power of $b$ it is not ideal to attempt to find roots for large
 values of $b$.  If particularly large roots are required then a factor method could be used instead.  For example,
-$a^{1/16}$ is equivalent to $\left (a^{1/4} \right)^{1/4}$.
+$a^{1/16}$ is equivalent to $\left (a^{1/4} \right)^{1/4}$ or simply 
+$\left ( \left ( \left ( a^{1/2} \right )^{1/2} \right )^{1/2} \right )^{1/2}$
 
 \chapter{Prime Numbers}
 \section{Trial Division}

bn_fast_mp_invmod.c

@@ -21,8 +21,7 @@
  * 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
+int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
-fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
 {
   mp_int  x, y, u, v, B, D;
   int     res, neg;

bn_fast_mp_montgomery_reduce.c

@@ -23,8 +23,7 @@
  *
  * Based on Algorithm 14.32 on pp.601 of HAC.
 */
-int
+int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
-fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
 {
   int     ix, res, olduse;
   mp_word W[MP_WARRAY];

bn_fast_s_mp_mul_digs.c

@@ -31,8 +31,7 @@
  * Based on Algorithm 14.12 on pp.595 of HAC.
  *
  */
-int
+int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
-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];
@@ -81,7 +80,7 @@ fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
   }
 
   /* store final carry */
-  W[ix] = _W;
+  W[ix] = _W & MP_MASK;
 
   /* setup dest */
   olduse  = c->used;

bn_fast_s_mp_mul_high_digs.c

@@ -24,8 +24,7 @@
  *
  * Based on Algorithm 14.12 on pp.595 of HAC.
  */
-int
+int fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
-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];
@@ -72,7 +71,7 @@ fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
   }
   
   /* store final carry */
-  W[ix] = _W;
+  W[ix] = _W & MP_MASK;
 
   /* setup dest */
   olduse  = c->used;

bn_fast_s_mp_sqr.c

@@ -101,7 +101,7 @@ int fast_s_mp_sqr (mp_int * a, mp_int * b)
       }
 
       /* store it */
-      W[ix] = _W;
+      W[ix] = _W & MP_MASK;
 
       /* make next carry */
       W1 = _W >> ((mp_word)DIGIT_BIT);

bn_mp_exptmod.c

@@ -65,21 +65,29 @@ int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 #endif
   }
 
+/* modified diminished radix reduction */
+#if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C)
+  if (mp_reduce_is_2k_l(P) == MP_YES) {
+     return s_mp_exptmod(G, X, P, Y, 1);
+  }
+#endif
+
 #ifdef BN_MP_DR_IS_MODULUS_C
   /* is it a DR modulus? */
   dr = mp_dr_is_modulus(P);
 #else
+  /* default to no */
   dr = 0;
 #endif
 
 #ifdef BN_MP_REDUCE_IS_2K_C
-  /* if not, is it a uDR modulus? */
+  /* if not, is it a unrestricted DR modulus? */
   if (dr == 0) {
      dr = mp_reduce_is_2k(P) << 1;
   }
 #endif
     
-  /* if the modulus is odd or dr != 0 use the fast method */
+  /* if the modulus is odd or dr != 0 use the montgomery method */
 #ifdef BN_MP_EXPTMOD_FAST_C
   if (mp_isodd (P) == 1 || dr !=  0) {
     return mp_exptmod_fast (G, X, P, Y, dr);
@@ -87,7 +95,7 @@ int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 #endif
 #ifdef BN_S_MP_EXPTMOD_C
     /* otherwise use the generic Barrett reduction technique */
-    return s_mp_exptmod (G, X, P, Y);
+    return s_mp_exptmod (G, X, P, Y, 0);
 #else
     /* no exptmod for evens */
     return MP_VAL;

bn_mp_exptmod_fast.c

@@ -29,8 +29,7 @@
    #define TAB_SIZE 256
 #endif
 
-int
+int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
-mp_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;

bn_mp_mul_d.c

@@ -57,8 +57,9 @@ mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
     u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
   }
 
-  /* store final carry [if any] */
+  /* store final carry [if any] and increment ix offset  */
   *tmpc++ = u;
+  ++ix;
 
   /* now zero digits above the top */
   while (ix++ < olduse) {

bn_mp_prime_random_ex.c

@@ -60,7 +60,7 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
 
    /* calc the maskOR_msb */
    maskOR_msb        = 0;
-   maskOR_msb_offset = (size - 2) >> 3;
+   maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
    if (flags & LTM_PRIME_2MSB_ON) {
       maskOR_msb     |= 1 << ((size - 2) & 7);
    } else if (flags & LTM_PRIME_2MSB_OFF) {
@@ -68,7 +68,7 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
    } 
 
    /* get the maskOR_lsb */
-   maskOR_lsb         = 0;
+   maskOR_lsb         = 1;
    if (flags & LTM_PRIME_BBS) {
       maskOR_lsb     |= 3;
    }

bn_mp_read_radix.c

@@ -16,7 +16,7 @@
  */
 
 /* read a string [ASCII] in a given radix */
-int mp_read_radix (mp_int * a, char *str, int radix)
+int mp_read_radix (mp_int * a, const char *str, int radix)
 {
   int     y, res, neg;
   char    ch;

bn_mp_reduce.c

@@ -19,8 +19,7 @@
  * precomputed via mp_reduce_setup.
  * From HAC pp.604 Algorithm 14.42
  */
-int
+int mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
-mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
 {
   mp_int  q;
   int     res, um = m->used;

bn_mp_reduce_2k.c

@@ -16,8 +16,7 @@
  */
 
 /* reduces a modulo n where n is of the form 2**p - d */
-int
+int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
-mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
 {
    mp_int q;
    int    p, res;

bn_mp_reduce_2k_l.c

@@ -0,0 +1,58 @@
+#include <tommath.h>
+#ifdef BN_MP_REDUCE_2K_L_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
+ */
+
+/* reduces a modulo n where n is of the form 2**p - d 
+   This differs from reduce_2k since "d" can be larger
+   than a single digit.
+*/
+int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
+{
+   mp_int q;
+   int    p, res;
+   
+   if ((res = mp_init(&q)) != MP_OKAY) {
+      return res;
+   }
+   
+   p = mp_count_bits(n);    
+top:
+   /* q = a/2**p, a = a mod 2**p */
+   if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
+      goto ERR;
+   }
+   
+   /* q = q * d */
+   if ((res = mp_mul(&q, d, &q)) != MP_OKAY) { 
+      goto ERR;
+   }
+   
+   /* a = a + q */
+   if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
+      goto ERR;
+   }
+   
+   if (mp_cmp_mag(a, n) != MP_LT) {
+      s_mp_sub(a, n, a);
+      goto top;
+   }
+   
+ERR:
+   mp_clear(&q);
+   return res;
+}
+
+#endif

bn_mp_reduce_2k_setup.c

@@ -16,8 +16,7 @@
  */
 
 /* determines the setup value */
-int 
+int mp_reduce_2k_setup(mp_int *a, mp_digit *d)
-mp_reduce_2k_setup(mp_int *a, mp_digit *d)
 {
    int res, p;
    mp_int tmp;

bn_mp_reduce_2k_setup_l.c

@@ -0,0 +1,40 @@
+#include <tommath.h>
+#ifdef BN_MP_REDUCE_2K_SETUP_L_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 the setup value */
+int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
+{
+   int    res;
+   mp_int tmp;
+   
+   if ((res = mp_init(&tmp)) != MP_OKAY) {
+      return res;
+   }
+   
+   if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
+      goto ERR;
+   }
+   
+   if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
+      goto ERR;
+   }
+   
+ERR:
+   mp_clear(&tmp);
+   return res;
+}
+#endif

bn_mp_reduce_is_2k.c

@@ -22,9 +22,9 @@ int mp_reduce_is_2k(mp_int *a)
    mp_digit iz;
    
    if (a->used == 0) {
-      return 0;
+      return MP_NO;
    } else if (a->used == 1) {
-      return 1;
+      return MP_YES;
    } else if (a->used > 1) {
       iy = mp_count_bits(a);
       iz = 1;
@@ -33,7 +33,7 @@ int mp_reduce_is_2k(mp_int *a)
       /* Test every bit from the second digit up, must be 1 */
       for (ix = DIGIT_BIT; ix < iy; ix++) {
           if ((a->dp[iw] & iz) == 0) {
-             return 0;
+             return MP_NO;
           }
           iz <<= 1;
           if (iz > (mp_digit)MP_MASK) {
@@ -42,7 +42,7 @@ int mp_reduce_is_2k(mp_int *a)
           }
       }
    }
-   return 1;
+   return MP_YES;
 }
 
 #endif

bn_mp_reduce_is_2k_l.c

@@ -0,0 +1,40 @@
+#include <tommath.h>
+#ifdef BN_MP_REDUCE_IS_2K_L_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 reduce_2k_l can be used */
+int mp_reduce_is_2k_l(mp_int *a)
+{
+   int ix, iy;
+   
+   if (a->used == 0) {
+      return MP_NO;
+   } else if (a->used == 1) {
+      return MP_YES;
+   } else if (a->used > 1) {
+      /* if more than half of the digits are -1 we're sold */
+      for (iy = ix = 0; ix < a->used; ix++) {
+          if (a->dp[ix] == MP_MASK) {
+              ++iy;
+          }
+      }
+      return (iy >= (a->used/2)) ? MP_YES : MP_NO;
+      
+   }
+   return MP_NO;
+}
+
+#endif

bn_mp_to_signed_bin.c

@@ -16,8 +16,7 @@
  */
 
 /* store in signed [big endian] format */
-int
+int mp_to_signed_bin (mp_int * a, unsigned char *b)
-mp_to_signed_bin (mp_int * a, unsigned char *b)
 {
   int     res;
 

bn_mp_to_signed_bin_n.c

@@ -0,0 +1,27 @@
+#include <tommath.h>
+#ifdef BN_MP_TO_SIGNED_BIN_N_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
+ */
+
+/* store in signed [big endian] format */
+int mp_to_signed_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen)
+{
+   if (*outlen < (unsigned long)mp_signed_bin_size(a)) {
+      return MP_VAL;
+   }
+   *outlen = mp_signed_bin_size(a);
+   return mp_to_signed_bin(a, b);
+}
+#endif

bn_mp_to_unsigned_bin.c

@@ -16,8 +16,7 @@
  */
 
 /* store in unsigned [big endian] format */
-int
+int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
-mp_to_unsigned_bin (mp_int * a, unsigned char *b)
 {
   int     x, res;
   mp_int  t;

bn_mp_to_unsigned_bin_n.c

@@ -0,0 +1,27 @@
+#include <tommath.h>
+#ifdef BN_MP_TO_UNSIGNED_BIN_N_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
+ */
+
+/* store in unsigned [big endian] format */
+int mp_to_unsigned_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen)
+{
+   if (*outlen < (unsigned long)mp_unsigned_bin_size(a)) {
+      return MP_VAL;
+   }
+   *outlen = mp_unsigned_bin_size(a);
+   return mp_to_unsigned_bin(a, b);
+}
+#endif

bn_mp_unsigned_bin_size.c

@@ -16,8 +16,7 @@
  */
 
 /* get the size for an unsigned equivalent */
-int
+int mp_unsigned_bin_size (mp_int * a)
-mp_unsigned_bin_size (mp_int * a)
 {
   int     size = mp_count_bits (a);
   return (size / 8 + ((size & 7) != 0 ? 1 : 0));

bn_s_mp_exptmod.c

@@ -21,11 +21,12 @@
    #define TAB_SIZE 256
 #endif
 
-int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
+int s_mp_exptmod (mp_int * G, mp_int * X, 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*,mp_int*,mp_int*);
 
   /* find window size */
   x = mp_count_bits (X);
@@ -72,9 +73,18 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
   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
    *
@@ -96,11 +106,14 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
   }
 
   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;
     }
-    if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
+
+    /* reduce modulo P */
+    if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
       goto LBL_MU;
     }
   }
@@ -112,7 +125,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
       goto LBL_MU;
     }
-    if ((err = mp_reduce (&M[x], P, &mu)) != MP_OKAY) {
+    if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
       goto LBL_MU;
     }
   }
@@ -161,7 +174,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
         goto LBL_RES;
       }
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
         goto LBL_RES;
       }
       continue;
@@ -178,7 +191,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
           goto LBL_RES;
         }
-        if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+        if ((err = redux (&res, P, &mu)) != MP_OKAY) {
           goto LBL_RES;
         }
       }
@@ -187,7 +200,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
         goto LBL_RES;
       }
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
         goto LBL_RES;
       }
 
@@ -205,7 +218,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
         goto LBL_RES;
       }
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
         goto LBL_RES;
       }
 
@@ -215,7 +228,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
           goto LBL_RES;
         }
-        if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+        if ((err = redux (&res, P, &mu)) != MP_OKAY) {
           goto LBL_RES;
         }
       }

bncore.c

@@ -20,11 +20,12 @@
  CPU                    /Compiler     /MUL CUTOFF/SQR CUTOFF
 -------------------------------------------------------------
  Intel P4 Northwood     /GCC v3.4.1   /        88/       128/LTM 0.32 ;-)
+ AMD Athlon64           /GCC v3.4.4   /        74/       124/LTM 0.34
  
 */
 
-int     KARATSUBA_MUL_CUTOFF = 88,      /* Min. number of digits before Karatsuba multiplication is used. */
+int     KARATSUBA_MUL_CUTOFF = 74,      /* Min. number of digits before Karatsuba multiplication is used. */
-        KARATSUBA_SQR_CUTOFF = 128,     /* Min. number of digits before Karatsuba squaring is used. */
+        KARATSUBA_SQR_CUTOFF = 124,     /* Min. number of digits before Karatsuba squaring is used. */
         
         TOOM_MUL_CUTOFF      = 350,      /* no optimal values of these are known yet so set em high */
         TOOM_SQR_CUTOFF      = 400; 

changes.txt

@@ -1,3 +1,15 @@
+February 12th, 2005
+v0.34  -- Fixed two more small errors in mp_prime_random_ex()
+       -- Fixed overflow in mp_mul_d() [Kevin Kenny]
+       -- Added mp_to_(un)signed_bin_n() functions which do bounds checking for ya [and report the size]
+       -- Added "large" diminished radix support.  Speeds up things like DSA where the moduli is of the form 2^k - P for some P < 2^(k/2) or so
+          Actually is faster than Montgomery on my AMD64 (and probably much faster on a P4)
+       -- Updated the manual a bit
+       -- Ok so I haven't done the textbook work yet... My current freelance gig has landed me in France till the 
+          end of Feb/05.  Once I get back I'll have tons of free time and I plan to go to town on the book.
+          As of this release the API will freeze.  At least until the book catches up with all the changes.  I welcome
+          bug reports but new algorithms will have to wait.
+
 December 23rd, 2004
 v0.33  -- Fixed "small" variant for mp_div() which would munge with negative dividends...
        -- Fixed bug in mp_prime_random_ex() which would set the most significant byte to zero when

demo/demo.c

@@ -12,6 +12,7 @@
 void ndraw(mp_int * a, char *name)
 {
    char buf[16000];
+
    printf("%s: ", name);
    mp_toradix(a, buf, 10);
    printf("%s\n", buf);
@@ -39,7 +40,9 @@ int lbit(void)
 int myrng(unsigned char *dst, int len, void *dat)
 {
    int x;
-   for (x = 0; x < len; x++) dst[x] = rand() & 0xFF;
+
+   for (x = 0; x < len; x++)
+      dst[x] = rand() & 0xFF;
    return len;
 }
 
@@ -49,8 +52,8 @@ int myrng(unsigned char *dst, int len, void *dat)
 int main(void)
 {
    mp_int a, b, c, d, e, f;
-   unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n, inv_n,
+   unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n,
-                 div2_n, mul2_n, add_d_n, sub_d_n, t;
+      gcd_n, lcm_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n, t;
    unsigned rr;
    int i, n, err, cnt, ix, old_kara_m, old_kara_s;
 
@@ -76,36 +79,37 @@ int main(void)
       }
    }
    mp_set_int(&a, 0);
-  if (mp_get_int(&a)!=0)
+   if (mp_get_int(&a) != 0) {
-  { printf("mp_get_int() bad result!\n");
+      printf("mp_get_int() bad result!\n");
       return 1;
    }
    mp_set_int(&a, 0xffffffff);
-  if (mp_get_int(&a)!=0xffffffff)
+   if (mp_get_int(&a) != 0xffffffff) {
-  { printf("mp_get_int() bad result!\n");
+      printf("mp_get_int() bad result!\n");
       return 1;
    }
-
    // test mp_sqrt
    printf("Testing: mp_sqrt\n");
    for (i = 0; i < 1000; ++i) {
-    printf("%6d\r", i); fflush(stdout);
+      printf("%6d\r", i);
+      fflush(stdout);
       n = (rand() & 15) + 1;
       mp_rand(&a, n);
-    if (mp_sqrt(&a,&b) != MP_OKAY)
+      if (mp_sqrt(&a, &b) != MP_OKAY) {
-    { printf("mp_sqrt() error!\n");
+	 printf("mp_sqrt() error!\n");
 	 return 1;
       }
       mp_n_root(&a, 2, &a);
-    if (mp_cmp_mag(&b,&a) != MP_EQ)
+      if (mp_cmp_mag(&b, &a) != MP_EQ) {
-    { printf("mp_sqrt() bad result!\n");
+	 printf("mp_sqrt() bad result!\n");
 	 return 1;
       }
    }
 
    printf("\nTesting: mp_is_square\n");
    for (i = 0; i < 1000; ++i) {
-    printf("%6d\r", i); fflush(stdout);
+      printf("%6d\r", i);
+      fflush(stdout);
 
       /* test mp_is_square false negatives */
       n = (rand() & 7) + 1;
@@ -136,8 +140,12 @@ int main(void)
 
    /* test for size */
    for (ix = 10; ix < 256; ix++) {
-       printf("Testing (not safe-prime): %9d bits    \r", ix); fflush(stdout);
+      printf("Testing (not safe-prime): %9d bits    \r", ix);
-       err = mp_prime_random_ex(&a, 8, ix, (rand()&1)?LTM_PRIME_2MSB_OFF:LTM_PRIME_2MSB_ON, myrng, NULL);
+      fflush(stdout);
+      err =
+	 mp_prime_random_ex(&a, 8, ix,
+			    (rand() & 1) ? LTM_PRIME_2MSB_OFF :
+			    LTM_PRIME_2MSB_ON, myrng, NULL);
       if (err != MP_OKAY) {
 	 printf("failed with err code %d\n", err);
 	 return EXIT_FAILURE;
@@ -149,8 +157,13 @@ int main(void)
    }
 
    for (ix = 16; ix < 256; ix++) {
-       printf("Testing (   safe-prime): %9d bits    \r", ix); fflush(stdout);
+      printf("Testing (   safe-prime): %9d bits    \r", ix);
-       err = mp_prime_random_ex(&a, 8, ix, ((rand()&1)?LTM_PRIME_2MSB_OFF:LTM_PRIME_2MSB_ON)|LTM_PRIME_SAFE, myrng, NULL);
+      fflush(stdout);
+      err =
+	 mp_prime_random_ex(&a, 8, ix,
+			    ((rand() & 1) ? LTM_PRIME_2MSB_OFF :
+			     LTM_PRIME_2MSB_ON) | LTM_PRIME_SAFE, myrng,
+			    NULL);
       if (err != MP_OKAY) {
 	 printf("failed with err code %d\n", err);
 	 return EXIT_FAILURE;
@@ -205,6 +218,7 @@ int main(void)
    printf("Testing mp_reduce_2k...\n");
    for (cnt = 3; cnt <= 128; ++cnt) {
       mp_digit tmp;
+
       mp_2expt(&a, cnt);
       mp_sub_d(&a, 2, &a);	/* a = 2**cnt - 2 */
 
@@ -214,7 +228,10 @@ int main(void)
       mp_reduce_2k_setup(&a, &tmp);
       printf("(%d)", tmp);
       for (ix = 0; ix < 1000; ix++) {
-           if (!(ix & 127)) {printf("."); fflush(stdout); }
+	 if (!(ix & 127)) {
+	    printf(".");
+	    fflush(stdout);
+	 }
 	 mp_rand(&b, (cnt / DIGIT_BIT + 1) * 2);
 	 mp_copy(&c, &b);
 	 mp_mod(&c, &a, &c);
@@ -232,7 +249,8 @@ int main(void)
    for (cnt = 0; cnt < 10000;) {
       mp_digit r1, r2;
 
-      if (!(++cnt & 127)) printf("%9d\r", cnt);
+      if (!(++cnt & 127))
+	 printf("%9d\r", cnt);
       mp_rand(&a, abs(rand()) % 128 + 1);
       mp_div(&a, &d, &b, &e);
       mp_div_3(&a, &c, &r2);
@@ -260,15 +278,20 @@ int main(void)
 
       rr = 0;
       do {
-         if (!(rr & 127)) { printf("%9lu\r", rr); fflush(stdout); }
+	 if (!(rr & 127)) {
-         mp_sqr(&b, &b); mp_add_d(&b, 1, &b);
+	    printf("%9lu\r", rr);
+	    fflush(stdout);
+	 }
+	 mp_sqr(&b, &b);
+	 mp_add_d(&b, 1, &b);
 	 mp_copy(&b, &c);
 
 	 mp_mod(&b, &a, &b);
 	 mp_dr_reduce(&c, &a, (((mp_digit) 1) << DIGIT_BIT) - a.dp[0]);
 
 	 if (mp_cmp(&b, &c) != MP_EQ) {
-            printf("Failed on trial %lu\n", rr); exit(-1);
+	    printf("Failed on trial %lu\n", rr);
+	    exit(-1);
 
 	 }
       } while (++rr < 500);
@@ -277,8 +300,58 @@ int main(void)
 
 #endif
 
+/* test the mp_reduce_2k_l code */
+#if 0
+#if 0
+/* first load P with 2^1024 - 0x2A434 B9FDEC95 D8F9D550 FFFFFFFF FFFFFFFF */
+   mp_2expt(&a, 1024);
+   mp_read_radix(&b, "2A434B9FDEC95D8F9D550FFFFFFFFFFFFFFFF", 16);
+   mp_sub(&a, &b, &a);
+#elif 1
+/*  p = 2^2048 - 0x1 00000000 00000000 00000000 00000000 4945DDBF 8EA2A91D 5776399B B83E188F  */
+   mp_2expt(&a, 2048);
+   mp_read_radix(&b,
+		 "1000000000000000000000000000000004945DDBF8EA2A91D5776399BB83E188F",
+		 16);
+   mp_sub(&a, &b, &a);
+#endif
+
+   mp_todecimal(&a, buf);
+   printf("p==%s\n", buf);
+/* now mp_reduce_is_2k_l() should return */
+   if (mp_reduce_is_2k_l(&a) != 1) {
+      printf("mp_reduce_is_2k_l() return 0, should be 1\n");
+      return EXIT_FAILURE;
+   }
+   mp_reduce_2k_setup_l(&a, &d);
+   /* now do a million square+1 to see if it varies */
+   mp_rand(&b, 64);
+   mp_mod(&b, &a, &b);
+   mp_copy(&b, &c);
+   printf("testing mp_reduce_2k_l...");
+   fflush(stdout);
+   for (cnt = 0; cnt < (1UL << 20); cnt++) {
+      mp_sqr(&b, &b);
+      mp_add_d(&b, 1, &b);
+      mp_reduce_2k_l(&b, &a, &d);
+      mp_sqr(&c, &c);
+      mp_add_d(&c, 1, &c);
+      mp_mod(&c, &a, &c);
+      if (mp_cmp(&b, &c) != MP_EQ) {
+	 printf("mp_reduce_2k_l() failed at step %lu\n", cnt);
+	 mp_tohex(&b, buf);
+	 printf("b == %s\n", buf);
+	 mp_tohex(&c, buf);
+	 printf("c == %s\n", buf);
+	 return EXIT_FAILURE;
+      }
+   }
+   printf("...Passed\n");
+#endif
+
    div2_n = mul2_n = inv_n = expt_n = lcm_n = gcd_n = add_n =
-   sub_n = mul_n = div_n = sqr_n = mul2d_n = div2d_n = cnt = add_d_n = sub_d_n= 0;
+      sub_n = mul_n = div_n = sqr_n = mul2d_n = div2d_n = cnt = add_d_n =
+      sub_d_n = 0;
 
    /* force KARA and TOOM to enable despite cutoffs */
    KARATSUBA_SQR_CUTOFF = KARATSUBA_MUL_CUTOFF = 110;
@@ -287,24 +360,51 @@ int main(void)
    for (;;) {
       /* randomly clear and re-init one variable, this has the affect of triming the alloc space */
       switch (abs(rand()) % 7) {
-           case 0:  mp_clear(&a); mp_init(&a); break;
+      case 0:
-           case 1:  mp_clear(&b); mp_init(&b); break;
+	 mp_clear(&a);
-           case 2:  mp_clear(&c); mp_init(&c); break;
+	 mp_init(&a);
-           case 3:  mp_clear(&d); mp_init(&d); break;
+	 break;
-           case 4:  mp_clear(&e); mp_init(&e); break;
+      case 1:
-           case 5:  mp_clear(&f); mp_init(&f); break;
+	 mp_clear(&b);
-           case 6:  break; /* don't clear any */
+	 mp_init(&b);
+	 break;
+      case 2:
+	 mp_clear(&c);
+	 mp_init(&c);
+	 break;
+      case 3:
+	 mp_clear(&d);
+	 mp_init(&d);
+	 break;
+      case 4:
+	 mp_clear(&e);
+	 mp_init(&e);
+	 break;
+      case 5:
+	 mp_clear(&f);
+	 mp_init(&f);
+	 break;
+      case 6:
+	 break;			/* don't clear any */
       }
 
 
-       printf("%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu ", add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n, expt_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n);
+      printf
+	 ("%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu ",
+	  add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n,
+	  expt_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n);
       fgets(cmd, 4095, stdin);
       cmd[strlen(cmd) - 1] = 0;
-       printf("%s  ]\r",cmd); fflush(stdout);
+      printf("%s  ]\r", cmd);
-       if (!strcmp(cmd, "mul2d")) { ++mul2d_n;
+      fflush(stdout);
-          fgets(buf, 4095, stdin); mp_read_radix(&a, buf, 64);
+      if (!strcmp(cmd, "mul2d")) {
-          fgets(buf, 4095, stdin); sscanf(buf, "%d", &rr);
+	 ++mul2d_n;
-          fgets(buf, 4095, stdin); mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 sscanf(buf, "%d", &rr);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 
 	 mp_mul_2d(&a, rr, &a);
 	 a.sign = b.sign;
@@ -314,29 +414,42 @@ int main(void)
 	    draw(&b);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "div2d")) { ++div2d_n;
+      } else if (!strcmp(cmd, "div2d")) {
-          fgets(buf, 4095, stdin); mp_read_radix(&a, buf, 64);
+	 ++div2d_n;
-          fgets(buf, 4095, stdin); sscanf(buf, "%d", &rr);
+	 fgets(buf, 4095, stdin);
-          fgets(buf, 4095, stdin); mp_read_radix(&b, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 sscanf(buf, "%d", &rr);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 
 	 mp_div_2d(&a, rr, &a, &e);
 	 a.sign = b.sign;
-          if (a.used == b.used && a.used == 0) { a.sign = b.sign = MP_ZPOS; }
+	 if (a.used == b.used && a.used == 0) {
+	    a.sign = b.sign = MP_ZPOS;
+	 }
 	 if (mp_cmp(&a, &b) != MP_EQ) {
 	    printf("div2d failed, rr == %d\n", rr);
 	    draw(&a);
 	    draw(&b);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "add")) { ++add_n;
+      } else if (!strcmp(cmd, "add")) {
-          fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++add_n;
-          fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-          fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_add(&d, &b, &d);
 	 if (mp_cmp(&c, &d) != MP_EQ) {
 	    printf("add %lu failure!\n", add_n);
-draw(&a);draw(&b);draw(&c);draw(&d);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
 	    return 0;
 	 }
 
@@ -365,104 +478,167 @@ draw(&a);draw(&b);draw(&c);draw(&d);
 	    return 0;
 	 }
 
-       } else if (!strcmp(cmd, "sub")) { ++sub_n;
+      } else if (!strcmp(cmd, "sub")) {
-          fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++sub_n;
-          fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-          fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_sub(&d, &b, &d);
 	 if (mp_cmp(&c, &d) != MP_EQ) {
 	    printf("sub %lu failure!\n", sub_n);
-draw(&a);draw(&b);draw(&c);draw(&d);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "mul")) { ++mul_n;
+      } else if (!strcmp(cmd, "mul")) {
-          fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++mul_n;
-          fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-          fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_mul(&d, &b, &d);
 	 if (mp_cmp(&c, &d) != MP_EQ) {
 	    printf("mul %lu failure!\n", mul_n);
-draw(&a);draw(&b);draw(&c);draw(&d);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "div")) { ++div_n;
+      } else if (!strcmp(cmd, "div")) {
-          fgets(buf, 4095, stdin); mp_read_radix(&a, buf, 64);
+	 ++div_n;
-          fgets(buf, 4095, stdin); mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-          fgets(buf, 4095, stdin); mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
-          fgets(buf, 4095, stdin); mp_read_radix(&d, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&d, buf, 64);
 
 	 mp_div(&a, &b, &e, &f);
 	 if (mp_cmp(&c, &e) != MP_EQ || mp_cmp(&d, &f) != MP_EQ) {
-             printf("div %lu %d, %d, failure!\n", div_n, mp_cmp(&c, &e), mp_cmp(&d, &f));
+	    printf("div %lu %d, %d, failure!\n", div_n, mp_cmp(&c, &e),
-draw(&a);draw(&b);draw(&c);draw(&d); draw(&e); draw(&f);
+		   mp_cmp(&d, &f));
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
+	    draw(&e);
+	    draw(&f);
 	    return 0;
 	 }
 
-       } else if (!strcmp(cmd, "sqr")) { ++sqr_n;
+      } else if (!strcmp(cmd, "sqr")) {
-          fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++sqr_n;
-          fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 	 mp_copy(&a, &c);
 	 mp_sqr(&c, &c);
 	 if (mp_cmp(&b, &c) != MP_EQ) {
 	    printf("sqr %lu failure!\n", sqr_n);
-draw(&a);draw(&b);draw(&c);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "gcd")) { ++gcd_n;
+      } else if (!strcmp(cmd, "gcd")) {
-          fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++gcd_n;
-          fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-          fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_gcd(&d, &b, &d);
 	 d.sign = c.sign;
 	 if (mp_cmp(&c, &d) != MP_EQ) {
 	    printf("gcd %lu failure!\n", gcd_n);
-draw(&a);draw(&b);draw(&c);draw(&d);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "lcm")) { ++lcm_n;
+      } else if (!strcmp(cmd, "lcm")) {
-             fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++lcm_n;
-             fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-             fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_lcm(&d, &b, &d);
 	 d.sign = c.sign;
 	 if (mp_cmp(&c, &d) != MP_EQ) {
 	    printf("lcm %lu failure!\n", lcm_n);
-   draw(&a);draw(&b);draw(&c);draw(&d);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "expt")) {  ++expt_n;
+      } else if (!strcmp(cmd, "expt")) {
-             fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++expt_n;
-             fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-             fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
-             fgets(buf, 4095, stdin);  mp_read_radix(&d, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&d, buf, 64);
 	 mp_copy(&a, &e);
 	 mp_exptmod(&e, &b, &c, &e);
 	 if (mp_cmp(&d, &e) != MP_EQ) {
 	    printf("expt %lu failure!\n", expt_n);
-   draw(&a);draw(&b);draw(&c);draw(&d); draw(&e);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
+	    draw(&e);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "invmod")) {  ++inv_n;
+      } else if (!strcmp(cmd, "invmod")) {
-             fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++inv_n;
-             fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
-             fgets(buf, 4095, stdin);  mp_read_radix(&c, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&c, buf, 64);
 	 mp_invmod(&a, &b, &d);
 	 mp_mulmod(&d, &a, &b, &e);
 	 if (mp_cmp_d(&e, 1) != MP_EQ) {
 	    printf("inv [wrong value from MPI?!] failure\n");
-                draw(&a);draw(&b);draw(&c);draw(&d);
+	    draw(&a);
+	    draw(&b);
+	    draw(&c);
+	    draw(&d);
 	    mp_gcd(&a, &b, &e);
 	    draw(&e);
 	    return 0;
 	 }
 
-       } else if (!strcmp(cmd, "div2")) { ++div2_n;
+      } else if (!strcmp(cmd, "div2")) {
-             fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++div2_n;
-             fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 	 mp_div_2(&a, &c);
 	 if (mp_cmp(&c, &b) != MP_EQ) {
 	    printf("div_2 %lu failure\n", div2_n);
@@ -471,9 +647,12 @@ draw(&a);draw(&b);draw(&c);draw(&d);
 	    draw(&c);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "mul2")) { ++mul2_n;
+      } else if (!strcmp(cmd, "mul2")) {
-             fgets(buf, 4095, stdin);  mp_read_radix(&a, buf, 64);
+	 ++mul2_n;
-             fgets(buf, 4095, stdin);  mp_read_radix(&b, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 	 mp_mul_2(&a, &c);
 	 if (mp_cmp(&c, &b) != MP_EQ) {
 	    printf("mul_2 %lu failure\n", mul2_n);
@@ -482,10 +661,14 @@ draw(&a);draw(&b);draw(&c);draw(&d);
 	    draw(&c);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "add_d")) { ++add_d_n;
+      } else if (!strcmp(cmd, "add_d")) {
-              fgets(buf, 4095, stdin); mp_read_radix(&a, buf, 64);
+	 ++add_d_n;
-              fgets(buf, 4095, stdin); sscanf(buf, "%d", &ix);
+	 fgets(buf, 4095, stdin);
-              fgets(buf, 4095, stdin); mp_read_radix(&b, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 sscanf(buf, "%d", &ix);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 	 mp_add_d(&a, ix, &c);
 	 if (mp_cmp(&b, &c) != MP_EQ) {
 	    printf("add_d %lu failure\n", add_d_n);
@@ -495,10 +678,14 @@ draw(&a);draw(&b);draw(&c);draw(&d);
 	    printf("d == %d\n", ix);
 	    return 0;
 	 }
-       } else if (!strcmp(cmd, "sub_d")) { ++sub_d_n;
+      } else if (!strcmp(cmd, "sub_d")) {
-              fgets(buf, 4095, stdin); mp_read_radix(&a, buf, 64);
+	 ++sub_d_n;
-              fgets(buf, 4095, stdin); sscanf(buf, "%d", &ix);
+	 fgets(buf, 4095, stdin);
-              fgets(buf, 4095, stdin); mp_read_radix(&b, buf, 64);
+	 mp_read_radix(&a, buf, 64);
+	 fgets(buf, 4095, stdin);
+	 sscanf(buf, "%d", &ix);
+	 fgets(buf, 4095, stdin);
+	 mp_read_radix(&b, buf, 64);
 	 mp_sub_d(&a, ix, &c);
 	 if (mp_cmp(&b, &c) != MP_EQ) {
 	    printf("sub_d %lu failure\n", sub_d_n);
@@ -512,4 +699,3 @@ draw(&a);draw(&b);draw(&c);draw(&d);
    }
    return 0;
 }
-

demo/timing.c

@@ -14,6 +14,7 @@ ulong64 _tt;
 void ndraw(mp_int * a, char *name)
 {
    char buf[4096];
+
    printf("%s: ", name);
    mp_toradix(a, buf, 64);
    printf("%s\n", buf);
@@ -44,13 +45,16 @@ static ulong64 TIMFUNC (void)
 #if defined __GNUC__
 #if defined(__i386__) || defined(__x86_64__)
    unsigned long long a;
-         __asm__ __volatile__ ("rdtsc\nmovl %%eax,%0\nmovl %%edx,4+%0\n"::"m"(a):"%eax","%edx");
+   __asm__ __volatile__("rdtsc\nmovl %%eax,%0\nmovl %%edx,4+%0\n"::
+			"m"(a):"%eax", "%edx");
    return a;
 #else /* gcc-IA64 version */
    unsigned long result;
    __asm__ __volatile__("mov %0=ar.itc":"=r"(result)::"memory");
+
    while (__builtin_expect((int) result == -1, 0))
       __asm__ __volatile__("mov %0=ar.itc":"=r"(result)::"memory");
+
    return result;
 #endif
 
@@ -77,7 +81,7 @@ static ulong64 TIMFUNC (void)
 int main(void)
 {
    ulong64 tt, gg, CLK_PER_SEC;
-   FILE *log, *logb, *logc;
+   FILE *log, *logb, *logc, *logd;
    mp_int a, b, c, d, e, f;
    int n, cnt, ix, old_kara_m, old_kara_s;
    unsigned rr;
@@ -100,7 +104,7 @@ int main(void)
    CLK_PER_SEC = TIMFUNC() - CLK_PER_SEC;
 
    printf("CLK_PER_SEC == %llu\n", CLK_PER_SEC);
-      
+   goto exptmod;
    log = fopen("logs/add.log", "w");
    for (cnt = 8; cnt <= 128; cnt += 8) {
       SLEEP;
@@ -112,10 +116,13 @@ int main(void)
 	 gg = TIMFUNC();
 	 DO(mp_add(&a, &b, &c));
 	 gg = (TIMFUNC() - gg) >> 1;
-            if (tt > gg) tt = gg;
+	 if (tt > gg)
+	    tt = gg;
       } while (++rr < 100000);
-         printf("Adding\t\t%4d-bit => %9llu/sec, %9llu cycles\n", mp_count_bits(&a), CLK_PER_SEC/tt, tt);
+      printf("Adding\t\t%4d-bit => %9llu/sec, %9llu cycles\n",
-         fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, tt); fflush(log);
+	     mp_count_bits(&a), CLK_PER_SEC / tt, tt);
+      fprintf(log, "%d %9llu\n", cnt * DIGIT_BIT, tt);
+      fflush(log);
    }
    fclose(log);
 
@@ -130,25 +137,29 @@ int main(void)
 	 gg = TIMFUNC();
 	 DO(mp_sub(&a, &b, &c));
 	 gg = (TIMFUNC() - gg) >> 1;
-            if (tt > gg) tt = gg;
+	 if (tt > gg)
+	    tt = gg;
       } while (++rr < 100000);
 
-         printf("Subtracting\t\t%4d-bit => %9llu/sec, %9llu cycles\n", mp_count_bits(&a), CLK_PER_SEC/tt, tt);
+      printf("Subtracting\t\t%4d-bit => %9llu/sec, %9llu cycles\n",
-         fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, tt);  fflush(log);
+	     mp_count_bits(&a), CLK_PER_SEC / tt, tt);
+      fprintf(log, "%d %9llu\n", cnt * DIGIT_BIT, tt);
+      fflush(log);
    }
    fclose(log);
 
    /* do mult/square twice, first without karatsuba and second with */
+ multtest:
    old_kara_m = KARATSUBA_MUL_CUTOFF;
    old_kara_s = KARATSUBA_SQR_CUTOFF;
-   for (ix = 0; ix < 1; ix++) {
+   for (ix = 0; ix < 2; ix++) {
       printf("With%s Karatsuba\n", (ix == 0) ? "out" : "");
 
       KARATSUBA_MUL_CUTOFF = (ix == 0) ? 9999 : old_kara_m;
       KARATSUBA_SQR_CUTOFF = (ix == 0) ? 9999 : old_kara_s;
 
       log = fopen((ix == 0) ? "logs/mult.log" : "logs/mult_kara.log", "w");
-      for (cnt = 4; cnt <= 288; cnt += 2) {
+      for (cnt = 4; cnt <= 10240 / DIGIT_BIT; cnt += 2) {
 	 SLEEP;
 	 mp_rand(&a, cnt);
 	 mp_rand(&b, cnt);
@@ -158,15 +169,18 @@ int main(void)
 	    gg = TIMFUNC();
 	    DO(mp_mul(&a, &b, &c));
 	    gg = (TIMFUNC() - gg) >> 1;
-            if (tt > gg) tt = gg;
+	    if (tt > gg)
+	       tt = gg;
 	 } while (++rr < 100);
-         printf("Multiplying\t%4d-bit => %9llu/sec, %9llu cycles\n", mp_count_bits(&a), CLK_PER_SEC/tt, tt);
+	 printf("Multiplying\t%4d-bit => %9llu/sec, %9llu cycles\n",
-         fprintf(log, "%d %9llu\n", mp_count_bits(&a), tt);  fflush(log);
+		mp_count_bits(&a), CLK_PER_SEC / tt, tt);
+	 fprintf(log, "%d %9llu\n", mp_count_bits(&a), tt);
+	 fflush(log);
       }
       fclose(log);
 
       log = fopen((ix == 0) ? "logs/sqr.log" : "logs/sqr_kara.log", "w");
-      for (cnt = 4; cnt <= 288; cnt += 2) {
+      for (cnt = 4; cnt <= 10240 / DIGIT_BIT; cnt += 2) {
 	 SLEEP;
 	 mp_rand(&a, cnt);
 	 rr = 0;
@@ -175,17 +189,25 @@ int main(void)
 	    gg = TIMFUNC();
 	    DO(mp_sqr(&a, &b));
 	    gg = (TIMFUNC() - gg) >> 1;
-            if (tt > gg) tt = gg;
+	    if (tt > gg)
+	       tt = gg;
 	 } while (++rr < 100);
-         printf("Squaring\t%4d-bit => %9llu/sec, %9llu cycles\n", mp_count_bits(&a), CLK_PER_SEC/tt, tt);
+	 printf("Squaring\t%4d-bit => %9llu/sec, %9llu cycles\n",
-         fprintf(log, "%d %9llu\n", mp_count_bits(&a), tt);  fflush(log);
+		mp_count_bits(&a), CLK_PER_SEC / tt, tt);
+	 fprintf(log, "%d %9llu\n", mp_count_bits(&a), tt);
+	 fflush(log);
       }
       fclose(log);
 
    }
+ exptmod:
 
    {
       char *primes[] = {
+	 /* 2K large moduli */
+	 "179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586239334100047359817950870678242457666208137217",
+	 "32317006071311007300714876688669951960444102669715484032130345427524655138867890893197201411522913463688717960921898019494119559150490921095088152386448283120630877367300996091750197750389652106796057638384067568276792218642619756161838094338476170470581645852036305042887575891541065808607552399123930385521914333389668342420684974786564569494856176035326322058077805659331026192708460314150258592864177116725943603718461857357598351152301645904403697613233287231227125684710820209725157101726931323469678542580656697935045997268352998638099733077152121140120031150424541696791951097529546801429027668869927491725169",
+	 "1044388881413152506691752710716624382579964249047383780384233483283953907971557456848826811934997558340890106714439262837987573438185793607263236087851365277945956976543709998340361590134383718314428070011855946226376318839397712745672334684344586617496807908705803704071284048740118609114467977783598029006686938976881787785946905630190260940599579453432823469303026696443059025015972399867714215541693835559885291486318237914434496734087811872639496475100189041349008417061675093668333850551032972088269550769983616369411933015213796825837188091833656751221318492846368125550225998300412344784862595674492194617023806505913245610825731835380087608622102834270197698202313169017678006675195485079921636419370285375124784014907159135459982790513399611551794271106831134090584272884279791554849782954323534517065223269061394905987693002122963395687782878948440616007412945674919823050571642377154816321380631045902916136926708342856440730447899971901781465763473223850267253059899795996090799469201774624817718449867455659250178329070473119433165550807568221846571746373296884912819520317457002440926616910874148385078411929804522981857338977648103126085902995208257421855249796721729039744118165938433694823325696642096892124547425283",
 	 /* 2K moduli mersenne primes */
 	 "6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151",
 	 "531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127",
@@ -216,6 +238,7 @@ int main(void)
       log = fopen("logs/expt.log", "w");
       logb = fopen("logs/expt_dr.log", "w");
       logc = fopen("logs/expt_2k.log", "w");
+      logd = fopen("logs/expt_2kl.log", "w");
       for (n = 0; primes[n]; n++) {
 	 SLEEP;
 	 mp_read_radix(&a, primes[n], 10);
@@ -234,7 +257,8 @@ int main(void)
 	    gg = TIMFUNC();
 	    DO(mp_exptmod(&c, &b, &a, &d));
 	    gg = (TIMFUNC() - gg) >> 1;
-            if (tt > gg) tt = gg;
+	    if (tt > gg)
+	       tt = gg;
 	 } while (++rr < 10);
 	 mp_sub_d(&a, 1, &e);
 	 mp_sub(&e, &b, &b);
@@ -245,13 +269,16 @@ int main(void)
 	    draw(&d);
 	    exit(0);
 	 }
-      printf("Exponentiating\t%4d-bit => %9llu/sec, %9llu cycles\n", mp_count_bits(&a), CLK_PER_SEC/tt, tt);
+	 printf("Exponentiating\t%4d-bit => %9llu/sec, %9llu cycles\n",
-      fprintf((n < 6) ? logc : (n < 13) ? logb : log, "%d %9llu\n", mp_count_bits(&a), tt);
+		mp_count_bits(&a), CLK_PER_SEC / tt, tt);
+	 fprintf(n < 4 ? logd : (n < 9) ? logc : (n < 16) ? logb : log,
+		 "%d %9llu\n", mp_count_bits(&a), tt);
       }
    }
    fclose(log);
    fclose(logb);
    fclose(logc);
+   fclose(logd);
 
    log = fopen("logs/invmod.log", "w");
    for (cnt = 4; cnt <= 128; cnt += 4) {
@@ -270,18 +297,19 @@ int main(void)
 	 gg = TIMFUNC();
 	 DO(mp_invmod(&b, &a, &c));
 	 gg = (TIMFUNC() - gg) >> 1;
-         if (tt > gg) tt = gg;
+	 if (tt > gg)
+	    tt = gg;
       } while (++rr < 1000);
       mp_mulmod(&b, &c, &a, &d);
       if (mp_cmp_d(&d, 1) != MP_EQ) {
 	 printf("Failed to invert\n");
 	 return 0;
       }
-      printf("Inverting mod\t%4d-bit => %9llu/sec, %9llu cycles\n", mp_count_bits(&a), CLK_PER_SEC/tt, tt);
+      printf("Inverting mod\t%4d-bit => %9llu/sec, %9llu cycles\n",
+	     mp_count_bits(&a), CLK_PER_SEC / tt, tt);
       fprintf(log, "%d %9llu\n", cnt * DIGIT_BIT, tt);
    }
    fclose(log);
 
    return 0;
 }
-

dep.pl

@@ -13,6 +13,8 @@ print CLASS "#if !(defined(LTM1) && defined(LTM2) && defined(LTM3))\n#if defined
 foreach my $filename (glob "bn*.c") {
    my $define = $filename;
 
+print "Processing $filename\n";
+
    # convert filename to upper case so we can use it as a define 
    $define =~ tr/[a-z]/[A-Z]/;
    $define =~ tr/\./_/;

etc/tune.c

@@ -10,13 +10,44 @@
  */
 #define TIMES (1UL<<14UL)
 
+/* RDTSC from Scott Duplichan */
+static ulong64 TIMFUNC (void)
+   {
+   #if defined __GNUC__
+      #if defined(__i386__) || defined(__x86_64__)
+         unsigned long long a;
+         __asm__ __volatile__ ("rdtsc\nmovl %%eax,%0\nmovl %%edx,4+%0\n"::"m"(a):"%eax","%edx");
+         return a;
+      #else /* gcc-IA64 version */
+         unsigned long result;
+         __asm__ __volatile__("mov %0=ar.itc" : "=r"(result) :: "memory");
+         while (__builtin_expect ((int) result == -1, 0))
+         __asm__ __volatile__("mov %0=ar.itc" : "=r"(result) :: "memory");
+         return result;
+      #endif
+
+   // Microsoft and Intel Windows compilers
+   #elif defined _M_IX86
+     __asm rdtsc
+   #elif defined _M_AMD64
+     return __rdtsc ();
+   #elif defined _M_IA64
+     #if defined __INTEL_COMPILER
+       #include <ia64intrin.h>
+     #endif
+      return __getReg (3116);
+   #else
+     #error need rdtsc function for this build
+   #endif
+   }
+
 
 #ifndef X86_TIMER
 
 /* generic ISO C timer */
 ulong64 LBL_T;
-void t_start(void) { LBL_T = clock(); }
+void t_start(void) { LBL_T = TIMFUNC(); }
-ulong64 t_read(void) { return clock() - LBL_T; }
+ulong64 t_read(void) { return TIMFUNC() - LBL_T; }
 
 #else
 extern void t_start(void);

logs/add.log

@@ -1,10 +1,10 @@
-480        88
+480        87
-960       113
+960       111
-1440       138
+1440       135
-1920       163
+1920       159
-2400       202
+2400       200
-2880       226
+2880       224
-3360       251
+3360       248
 3840       272
 4320       296
 4800       320

logs/expt.log

@@ -1,7 +1,7 @@
-513   1499509
+513   1489160
-769   3682671
+769   3688476
-1025   8098887
+1025   8162061
-2049  49332743
+2049  49260015
-2561  89647783
+2561  89579052
-3073 149440713
+3073 148797060
-4097 326135364
+4097 324449263

logs/expt_2k.log

@@ -1,6 +1,5 @@
-521   1423346
+607   2272809
-607   1841305
+1279   9557382
-1279   8375656
+2203  36250309
-2203  34104708
+3217  87666486
-3217  83830729
+4253 174168369
-4253 167916804

logs/expt_2kl.log

@@ -0,0 +1,4 @@
+1024   6954080
+2048  35993987
+4096 176068521
+521   1683720

logs/expt_dr.log

@@ -1,7 +1,7 @@
-532   1803110
+532   1989592
-784   3607375
+784   3898697
-1036   6089790
+1036   6519700
-1540  14739797
+1540  15676650
-2072  33251589
+2072  33128187
-3080  82794331
+3080  82963362
-4116 165212734
+4116 168358337

logs/mult.log

@@ -1,143 +1,84 @@
-271       580
+271       555
-390       861
+390       855
-511      1177
+508      1161
-630      1598
+631      1605
-749      2115
+749      2117
-871      2670
+871      2687
-991      3276
+991      3329
-1111      3987
+1108      4084
-1231      4722
+1231      4786
-1351      5474
+1351      5624
-1471      6281
+1470      6392
-1589      7126
+1586      7364
-1710      8114
+1710      8218
-1831      8988
+1830      9255
-1946     10038
+1951     10217
-2071     10995
+2067     11461
-2188     12286
+2191     12463
-2310     13152
+2308     13677
-2430     14480
+2430     14800
-2549     15521
+2551     16232
-2671     17171
+2671     17460
-2790     18081
+2791     18899
-2911     19754
+2902     20247
-3031     20809
+3028     21902
-3150     22849
+3151     23240
-3269     23757
+3267     24927
-3391     25772
+3391     26441
-3508     26832
+3511     28277
-3631     29304
+3631     29838
-3750     30149
+3749     31751
-3865     32581
+3869     33673
-3988     33644
+3989     35431
-4111     36565
+4111     37518
-4231     37309
+4231     39426
-4351     40152
+4349     41504
-4471     41188
+4471     43567
-4590     44658
+4591     45786
-4710     45256
+4711     47876
-4827     48538
+4831     50299
-4951     49490
+4951     52427
-5070     53472
+5071     54785
-5190     53902
+5189     57241
-5308     57619
+5307     59730
-5431     58509
+5431     62194
-5550     63044
+5551     64761
-5664     63333
+5670     67322
-5791     67542
+5789     70073
-5911     68279
+5907     72663
-6028     73477
+6030     75437
-6150     73475
+6151     78242
-6271     78189
+6268     81202
-6390     78842
+6389     83948
-6510     84691
+6509     86985
-6631     84444
+6631     89903
-6751     89721
+6747     93184
-6871     90186
+6869     96044
-6991     96665
+6991     99286
-7111     96119
+7109    102395
-7231    101937
+7229    105917
-7350    102212
+7351    108940
-7471    109439
+7470    112490
-7591    108491
+7589    115702
-7709    114965
+7711    119508
-7829    115025
+7831    122632
-7951    123002
+7951    126410
-8071    121630
+8071    129808
-8190    128725
+8190    133895
-8311    128536
+8311    137146
-8430    137298
+8431    141218
-8550    135568
+8549    144732
-8671    143265
+8667    149131
-8791    142793
+8790    152462
-8911    152432
+8911    156754
-9030    150202
+9030    160479
-9151    158616
+9149    165138
-9271    157848
+9271    168601
-9391    168374
+9391    173185
-9511    165651
+9511    176988
-9627    174775
+9627    181976
-9750    173375
+9751    185539
-9871    185067
+9870    190388
-9985    181845
+9991    194335
-10111    191708
+10110    199605
-10229    190239
+10228    203298
-10351    202585
-10467    198704
-10591    209193
-10711    207322
-10831    220842
-10950    215882
-11071    227761
-11191    225501
-11311    239669
-11430    234809
-11550    243511
-11671    255947
-11791    255243
-11906    267828
-12029    263437
-12149    276571
-12270    275579
-12390    288963
-12510    284001
-12631    298196
-12751    297018
-12869    310848
-12990    305369
-13111    319086
-13230    318940
-13349    333685
-13471    327495
-13588    343678
-13711    341817
-13831    357181
-13948    350440
-14071    367526
-14189    365330
-14311    381551
-14429    374149
-14549    392203
-14670    389764
-14791    406761
-14910    398652
-15026    417718
-15150    414733
-15269    432759
-15390   1037071
-15511   1053454
-15631   1069198
-15748   1086164
-15871   1112820
-15991   1129676
-16111   1145924
-16230   1163016
-16345   1179911
-16471   1197048
-16586   1214352
-16711   1232095
-16829   1249338
-16947   1266987
-17071   1284181
-17188   1302521
-17311   1320539

logs/mult_kara.log

@@ -1,33 +1,84 @@
-924     16686
+271       560
-1146     25334
+391       870
-1371     35304
+511      1159
-1591     47122
+631      1605
-1820     61500
+750      2111
-2044     75254
+871      2737
-2266     91732
+991      3361
-2492    111656
+1111      4054
-2716    129428
+1231      4778
-2937    147508
+1351      5600
-3164    167758
+1471      6404
-3388    188248
+1591      7323
-3612    210826
+1710      8255
-3836    233814
+1831      9239
-4059    256898
+1948     10257
-4284    280210
+2070     11397
-4508    310372
+2190     12531
-4731    333902
+2308     13665
-4955    376502
+2429     14870
-5179    402854
+2550     16175
-5404    432004
+2671     17539
-5626    459010
+2787     18879
-5849    491868
+2911     20350
-6076    520550
+3031     21807
-6300    547400
+3150     23415
-6524    575968
+3270     24897
-6747    608482
+3388     26567
-6971    642850
+3511     28205
-7196    673670
+3627     30076
-7419    710680
+3751     31744
-7644    743942
+3869     33657
-7868    780394
+3991     35425
-8092    817342
+4111     37522
+4229     39363
+4351     41503
+4470     43491
+4590     45827
+4711     47795
+4828     50166
+4951     52318
+5070     54911
+5191     57036
+5308     58237
+5431     60248
+5551     62678
+5671     64786
+5791     67294
+5908     69343
+6031     71607
+6151     74166
+6271     76590
+6391     78734
+6511     81175
+6631     83742
+6750     86403
+6868     88873
+6990     91150
+7110     94211
+7228     96922
+7351     99445
+7469    102216
+7589    104968
+7711    108113
+7827    110758
+7950    113714
+8071    116511
+8186    119643
+8310    122679
+8425    125581
+8551    128715
+8669    131778
+8788    135116
+8910    138138
+9031    141628
+9148    144754
+9268    148367
+9391    151551
+9511    155033
+9631    158652
+9751    162125
+9871    165248
+9988    168627
+10111    172427
+10231    176412

logs/sqr.log

@@ -1,143 +1,84 @@
-271       552
+265       562
-389       883
+389       882
-510      1191
+509      1207
-629      1572
+631      1572
-750      1996
+750      1990
-863      2428
+859      2433
-991      2891
+991      2894
-1108      3539
+1109      3555
-1231      4182
+1230      4228
-1351      4980
+1350      5018
-1471      5771
+1471      5805
-1590      6551
+1591      6579
-1711      7313
+1709      7415
-1830      8240
+1829      8329
-1951      9184
+1949      9225
-2070     10087
+2071     10139
-2191     11140
+2188     11239
-2311     12111
+2309     12178
-2431     13219
+2431     13212
-2550     14247
+2551     14294
-2669     15353
+2671     15551
-2791     16446
+2791     16512
-2911     17692
+2911     17718
-3029     18848
+3030     18876
-3151     20028
+3150     20259
-3268     21282
+3270     21374
-3391     22696
+3391     22650
-3511     23971
+3511     23948
-3631     25303
+3631     25493
-3751     26675
+3750     26756
-3871     28245
+3870     28225
-3990     29736
+3989     29705
-4111     31124
+4110     31409
-4229     32714
+4230     32834
-4347     34397
+4351     34327
-4471     35877
+4471     35818
-4587     37269
+4591     37636
-4710     39011
+4711     39228
-4831     40884
+4830     40868
-4950     42501
+4949     42393
-5070     44005
+5070     44541
-5191     46026
+5191     46269
-5310     48168
+5310     48162
-5431     49801
+5429     49728
-5551     51385
+5548     51985
-5671     53604
+5671     53948
-5787     55942
+5791     55885
-5910     57757
+5910     57584
-6031     59391
+6031     60082
-6151     61754
+6150     62239
-6271     64234
+6270     64309
-6390     66110
+6390     66014
-6511     67845
+6511     68766
-6627     70474
+6631     71012
-6751     73113
+6750     73172
-6871     75064
+6871     74952
-6990     76940
+6991     77909
-7111     79681
+7111     80371
-7230     82548
+7231     82666
-7351     84597
+7351     84531
-7471     86507
+7469     87698
-7591     89497
+7589     90318
-7711    225216
+7711    225384
-7831    232192
+7830    232428
-7951    239583
+7950    240009
-8071    247302
+8070    246522
-8191    255497
+8190    253662
-8308    261587
+8310    260961
-8431    271490
+8431    269253
-8550    279492
+8549    275743
-8671    286927
+8671    283769
-8790    294680
+8789    290811
-8910    302974
+8911    300034
-9030    311300
+9030    306873
-9150    318635
+9149    315085
-9271    326740
+9270    323944
-9390    335304
+9390    332390
-9511    344297
+9508    337519
-9630    352056
+9631    348986
-9748    358652
+9749    356904
-9870    369723
+9871    367013
-9991    379119
+9989    373831
-10111    386982
+10108    381033
-10231    396075
+10230    393475
-10349    404396
-10470    415375
-10590    424146
-10711    433390
-10829    442662
-10950    453238
-11071    462178
-11186    469811
-11311    482529
-11431    493214
-11550    503210
-11671    513486
-11791    524244
-11911    535277
-12031    544872
-12151    555695
-12271    566893
-12391    578385
-12510    588658
-12628    596914
-12751    611324
-12871    623437
-12991    633907
-13110    645605
-13231    657684
-13351    670037
-13471    680939
-13591    693047
-13710    705363
-13829    718178
-13949    727930
-14069    739641
-14190    754817
-14310    768192
-14431    779875
-14551    792655
-14667    802847
-14791    819806
-14911    831684
-15031    844936
-15151    858813
-15270    873037
-15387    882123
-15510    899117
-15631    913465
-15750    927989
-15870    940790
-15991    954948
-16110    969483
-16231    984544
-16350    997837
-16470   1012445
-16590   1027834
-16710   1043032
-16831   1056394
-16951   1071408
-17069   1097263
-17191   1113364
-17306   1123650

logs/sqr_kara.log

@@ -1,33 +1,84 @@
-922     11272
+271       560
-1148     16004
+388       878
-1370     21958
+511      1179
-1596     28684
+629      1625
-1817     37832
+751      1988
-2044     46386
+871      2423
-2262     56218
+989      2896
-2492     66388
+1111      3561
-2716     77478
+1231      4209
-2940     89380
+1350      5015
-3163    103680
+1470      5804
-3385    116274
+1591      6556
-3612    135334
+1709      7420
-3836    151332
+1831      8263
-4057    164938
+1951      9173
-4284    183178
+2070     10153
-4508    198864
+2191     11229
-4731    215222
+2310     12167
-4954    231986
+2431     13211
-5180    251660
+2550     14309
-5404    269414
+2671     15524
-5626    288454
+2788     16525
-5850    307806
+2910     17712
-6076    329458
+3028     18822
-6299    347726
+3148     20220
-6523    369864
+3271     21343
-6748    387832
+3391     22652
-6971    413010
+3511     23944
-7194    453310
+3630     25485
-7415    476936
+3750     26778
-7643    497118
+3868     28201
-7867    521394
+3990     29653
-8091    540224
+4111     31393
+4225     32841
+4350     34328
+4471     35786
+4590     37652
+4711     39245
+4830     40876
+4951     42433
+5068     44547
+5191     46321
+5311     48140
+5430     49727
+5550     52034
+5671     53954
+5791     55921
+5908     57597
+6031     60084
+6148     62226
+6270     64295
+6390     66045
+6511     68779
+6629     71003
+6751     73169
+6871     74992
+6991     77895
+7110     80376
+7231     82628
+7351     84468
+7470     87664
+7591     90284
+7711     91352
+7828     93995
+7950     96276
+8071     98691
+8190    101256
+8308    103631
+8431    105222
+8550    108343
+8671    110281
+8787    112764
+8911    115397
+9031    117690
+9151    120266
+9271    122715
+9391    124624
+9510    127937
+9630    130313
+9750    132914
+9871    136129
+9991    138517
+10108    141525
+10231    144225

logs/sub.log

@@ -1,16 +1,16 @@
-480        87
+480        94
-960       114
+960       116
-1440       139
+1440       140
-1920       159
+1920       164
-2400       204
+2400       205
-2880       228
+2880       229
-3360       250
+3360       253
-3840       273
+3840       277
-4320       300
+4320       299
 4800       321
-5280       348
+5280       345
-5760       370
+5760       371
-6240       393
+6240       395
-6720       420
+6720       419
-7200       444
+7200       441
-7680       466
+7680       465

makefile

@@ -3,7 +3,7 @@
 #Tom St Denis
 
 #version of library 
-VERSION=0.33
+VERSION=0.34
 
 CFLAGS  +=  -I./ -Wall -W -Wshadow -Wsign-compare
 
@@ -57,11 +57,13 @@ bn_mp_prime_is_prime.o bn_mp_prime_next_prime.o bn_mp_dr_reduce.o \
 bn_mp_dr_is_modulus.o bn_mp_dr_setup.o bn_mp_reduce_setup.o \
 bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_div_3.o bn_s_mp_exptmod.o \
 bn_mp_reduce_2k.o bn_mp_reduce_is_2k.o bn_mp_reduce_2k_setup.o \
+bn_mp_reduce_2k_l.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_2k_setup_l.o \
 bn_mp_radix_smap.o bn_mp_read_radix.o bn_mp_toradix.o bn_mp_radix_size.o \
 bn_mp_fread.o bn_mp_fwrite.o bn_mp_cnt_lsb.o bn_error.o \
 bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_exteuclid.o bn_mp_toradix_n.o \
 bn_mp_prime_random_ex.o bn_mp_get_int.o bn_mp_sqrt.o bn_mp_is_square.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o
+bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o \
+bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin_n.o
 
 libtommath.a:  $(OBJECTS)
 	$(AR) $(ARFLAGS) libtommath.a $(OBJECTS)

makefile.bcc

@@ -27,11 +27,13 @@ bn_mp_prime_is_prime.obj bn_mp_prime_next_prime.obj bn_mp_dr_reduce.obj \
 bn_mp_dr_is_modulus.obj bn_mp_dr_setup.obj bn_mp_reduce_setup.obj \
 bn_mp_toom_mul.obj bn_mp_toom_sqr.obj bn_mp_div_3.obj bn_s_mp_exptmod.obj \
 bn_mp_reduce_2k.obj bn_mp_reduce_is_2k.obj bn_mp_reduce_2k_setup.obj \
+bn_mp_reduce_2k_l.obj bn_mp_reduce_is_2k_l.obj bn_mp_reduce_2k_setup_l.obj \
 bn_mp_radix_smap.obj bn_mp_read_radix.obj bn_mp_toradix.obj bn_mp_radix_size.obj \
 bn_mp_fread.obj bn_mp_fwrite.obj bn_mp_cnt_lsb.obj bn_error.obj \
 bn_mp_init_multi.obj bn_mp_clear_multi.obj bn_mp_exteuclid.obj bn_mp_toradix_n.obj \
 bn_mp_prime_random_ex.obj bn_mp_get_int.obj bn_mp_sqrt.obj bn_mp_is_square.obj \
-bn_mp_init_set.obj bn_mp_init_set_int.obj bn_mp_invmod_slow.obj bn_mp_prime_rabin_miller_trials.obj
+bn_mp_init_set.obj bn_mp_init_set_int.obj bn_mp_invmod_slow.obj bn_mp_prime_rabin_miller_trials.obj \
+bn_mp_to_signed_bin_n.obj bn_mp_to_unsigned_bin_n.obj
 
 TARGET = libtommath.lib
 

makefile.cygwin_dll

@@ -32,11 +32,13 @@ bn_mp_prime_is_prime.o bn_mp_prime_next_prime.o bn_mp_dr_reduce.o \
 bn_mp_dr_is_modulus.o bn_mp_dr_setup.o bn_mp_reduce_setup.o \
 bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_div_3.o bn_s_mp_exptmod.o \
 bn_mp_reduce_2k.o bn_mp_reduce_is_2k.o bn_mp_reduce_2k_setup.o \
+bn_mp_reduce_2k_l.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_2k_setup_l.o \
 bn_mp_radix_smap.o bn_mp_read_radix.o bn_mp_toradix.o bn_mp_radix_size.o \
 bn_mp_fread.o bn_mp_fwrite.o bn_mp_cnt_lsb.o bn_error.o \
 bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_exteuclid.o bn_mp_toradix_n.o \
 bn_mp_prime_random_ex.o bn_mp_get_int.o bn_mp_sqrt.o bn_mp_is_square.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o
+bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o \
+bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin_n.o
 
 # make a Windows DLL via Cygwin
 windll:  $(OBJECTS)

makefile.icc

@@ -59,11 +59,13 @@ bn_mp_prime_is_prime.o bn_mp_prime_next_prime.o bn_mp_dr_reduce.o \
 bn_mp_dr_is_modulus.o bn_mp_dr_setup.o bn_mp_reduce_setup.o \
 bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_div_3.o bn_s_mp_exptmod.o \
 bn_mp_reduce_2k.o bn_mp_reduce_is_2k.o bn_mp_reduce_2k_setup.o \
+bn_mp_reduce_2k_l.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_2k_setup_l.o \
 bn_mp_radix_smap.o bn_mp_read_radix.o bn_mp_toradix.o bn_mp_radix_size.o \
 bn_mp_fread.o bn_mp_fwrite.o bn_mp_cnt_lsb.o bn_error.o \
 bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_exteuclid.o bn_mp_toradix_n.o \
 bn_mp_prime_random_ex.o bn_mp_get_int.o bn_mp_sqrt.o bn_mp_is_square.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o
+bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o \
+bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin_n.o
 
 libtommath.a:  $(OBJECTS)
 	$(AR) $(ARFLAGS) libtommath.a $(OBJECTS)

makefile.msvc

@@ -26,11 +26,13 @@ bn_mp_prime_is_prime.obj bn_mp_prime_next_prime.obj bn_mp_dr_reduce.obj \
 bn_mp_dr_is_modulus.obj bn_mp_dr_setup.obj bn_mp_reduce_setup.obj \
 bn_mp_toom_mul.obj bn_mp_toom_sqr.obj bn_mp_div_3.obj bn_s_mp_exptmod.obj \
 bn_mp_reduce_2k.obj bn_mp_reduce_is_2k.obj bn_mp_reduce_2k_setup.obj \
+bn_mp_reduce_2k_l.obj bn_mp_reduce_is_2k_l.obj bn_mp_reduce_2k_setup_l.obj \
 bn_mp_radix_smap.obj bn_mp_read_radix.obj bn_mp_toradix.obj bn_mp_radix_size.obj \
 bn_mp_fread.obj bn_mp_fwrite.obj bn_mp_cnt_lsb.obj bn_error.obj \
 bn_mp_init_multi.obj bn_mp_clear_multi.obj bn_mp_exteuclid.obj bn_mp_toradix_n.obj \
 bn_mp_prime_random_ex.obj bn_mp_get_int.obj bn_mp_sqrt.obj bn_mp_is_square.obj \
-bn_mp_init_set.obj bn_mp_init_set_int.obj bn_mp_invmod_slow.obj bn_mp_prime_rabin_miller_trials.obj
+bn_mp_init_set.obj bn_mp_init_set_int.obj bn_mp_invmod_slow.obj bn_mp_prime_rabin_miller_trials.obj \
+bn_mp_to_signed_bin_n.obj bn_mp_to_unsigned_bin_n.obj
 
 library: $(OBJECTS)
 	lib /out:tommath.lib $(OBJECTS)

makefile.shared

@@ -1,7 +1,7 @@
 #Makefile for GCC
 #
 #Tom St Denis
-VERSION=0:33
+VERSION=0:34
 
 CC = libtool --mode=compile gcc
 CFLAGS  +=  -I./ -Wall -W -Wshadow -Wsign-compare
@@ -53,11 +53,14 @@ bn_mp_prime_is_prime.o bn_mp_prime_next_prime.o bn_mp_dr_reduce.o \
 bn_mp_dr_is_modulus.o bn_mp_dr_setup.o bn_mp_reduce_setup.o \
 bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_div_3.o bn_s_mp_exptmod.o \
 bn_mp_reduce_2k.o bn_mp_reduce_is_2k.o bn_mp_reduce_2k_setup.o \
+bn_mp_reduce_2k_l.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_2k_setup_l.o \
 bn_mp_radix_smap.o bn_mp_read_radix.o bn_mp_toradix.o bn_mp_radix_size.o \
 bn_mp_fread.o bn_mp_fwrite.o bn_mp_cnt_lsb.o bn_error.o \
 bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_exteuclid.o bn_mp_toradix_n.o \
 bn_mp_prime_random_ex.o bn_mp_get_int.o bn_mp_sqrt.o bn_mp_is_square.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o
+bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o \
+bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin_n.o
+
 
 libtommath.la:  $(OBJECTS)
 	libtool --mode=link gcc *.lo -o libtommath.la -rpath $(LIBPATH) -version-info $(VERSION)

pre_gen/mpi.c

@@ -69,8 +69,7 @@ char *mp_error_to_string(int code)
  * 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
+int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
-fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
 {
   mp_int  x, y, u, v, B, D;
   int     res, neg;
@@ -220,8 +219,7 @@ LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);
  *
  * Based on Algorithm 14.32 on pp.601 of HAC.
 */
-int
+int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
-fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
 {
   int     ix, res, olduse;
   mp_word W[MP_WARRAY];
@@ -401,8 +399,7 @@ fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
  * Based on Algorithm 14.12 on pp.595 of HAC.
  *
  */
-int
+int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
-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];
@@ -451,7 +448,7 @@ fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
   }
 
   /* store final carry */
-  W[ix] = _W;
+  W[ix] = _W & MP_MASK;
 
   /* setup dest */
   olduse  = c->used;
@@ -504,8 +501,7 @@ fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
  *
  * Based on Algorithm 14.12 on pp.595 of HAC.
  */
-int
+int fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
-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];
@@ -552,7 +548,7 @@ fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
   }
   
   /* store final carry */
-  W[ix] = _W;
+  W[ix] = _W & MP_MASK;
 
   /* setup dest */
   olduse  = c->used;
@@ -683,7 +679,7 @@ int fast_s_mp_sqr (mp_int * a, mp_int * b)
       }
 
       /* store it */
-      W[ix] = _W;
+      W[ix] = _W & MP_MASK;
 
       /* make next carry */
       W1 = _W >> ((mp_word)DIGIT_BIT);
@@ -2467,21 +2463,29 @@ int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 #endif
   }
 
+/* modified diminished radix reduction */
+#if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C)
+  if (mp_reduce_is_2k_l(P) == MP_YES) {
+     return s_mp_exptmod(G, X, P, Y, 1);
+  }
+#endif
+
 #ifdef BN_MP_DR_IS_MODULUS_C
   /* is it a DR modulus? */
   dr = mp_dr_is_modulus(P);
 #else
+  /* default to no */
   dr = 0;
 #endif
 
 #ifdef BN_MP_REDUCE_IS_2K_C
-  /* if not, is it a uDR modulus? */
+  /* if not, is it a unrestricted DR modulus? */
   if (dr == 0) {
      dr = mp_reduce_is_2k(P) << 1;
   }
 #endif
     
-  /* if the modulus is odd or dr != 0 use the fast method */
+  /* if the modulus is odd or dr != 0 use the montgomery method */
 #ifdef BN_MP_EXPTMOD_FAST_C
   if (mp_isodd (P) == 1 || dr !=  0) {
     return mp_exptmod_fast (G, X, P, Y, dr);
@@ -2489,7 +2493,7 @@ int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 #endif
 #ifdef BN_S_MP_EXPTMOD_C
     /* otherwise use the generic Barrett reduction technique */
-    return s_mp_exptmod (G, X, P, Y);
+    return s_mp_exptmod (G, X, P, Y, 0);
 #else
     /* no exptmod for evens */
     return MP_VAL;
@@ -2535,8 +2539,7 @@ int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
    #define TAB_SIZE 256
 #endif
 
-int
+int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
-mp_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;
@@ -4989,8 +4992,9 @@ mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
     u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
   }
 
-  /* store final carry [if any] */
+  /* store final carry [if any] and increment ix offset  */
   *tmpc++ = u;
+  ++ix;
 
   /* now zero digits above the top */
   while (ix++ < olduse) {
@@ -5847,7 +5851,7 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
 
    /* calc the maskOR_msb */
    maskOR_msb        = 0;
-   maskOR_msb_offset = (size - 2) >> 3;
+   maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
    if (flags & LTM_PRIME_2MSB_ON) {
       maskOR_msb     |= 1 << ((size - 2) & 7);
    } else if (flags & LTM_PRIME_2MSB_OFF) {
@@ -5855,7 +5859,7 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
    } 
 
    /* get the maskOR_lsb */
-   maskOR_lsb         = 0;
+   maskOR_lsb         = 1;
    if (flags & LTM_PRIME_BBS) {
       maskOR_lsb     |= 3;
    }
@@ -6080,7 +6084,7 @@ mp_rand (mp_int * a, int digits)
  */
 
 /* read a string [ASCII] in a given radix */
-int mp_read_radix (mp_int * a, char *str, int radix)
+int mp_read_radix (mp_int * a, const char *str, int radix)
 {
   int     y, res, neg;
   char    ch;
@@ -6263,8 +6267,7 @@ mp_read_unsigned_bin (mp_int * a, unsigned char *b, int c)
  * precomputed via mp_reduce_setup.
  * From HAC pp.604 Algorithm 14.42
  */
-int
+int mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
-mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
 {
   mp_int  q;
   int     res, um = m->used;
@@ -6361,8 +6364,7 @@ CLEANUP:
  */
 
 /* reduces a modulo n where n is of the form 2**p - d */
-int
+int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
-mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
 {
    mp_int q;
    int    p, res;
@@ -6404,6 +6406,68 @@ ERR:
 
 /* End: bn_mp_reduce_2k.c */
 
+/* Start: bn_mp_reduce_2k_l.c */
+#include <tommath.h>
+#ifdef BN_MP_REDUCE_2K_L_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
+ */
+
+/* reduces a modulo n where n is of the form 2**p - d 
+   This differs from reduce_2k since "d" can be larger
+   than a single digit.
+*/
+int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
+{
+   mp_int q;
+   int    p, res;
+   
+   if ((res = mp_init(&q)) != MP_OKAY) {
+      return res;
+   }
+   
+   p = mp_count_bits(n);    
+top:
+   /* q = a/2**p, a = a mod 2**p */
+   if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
+      goto ERR;
+   }
+   
+   /* q = q * d */
+   if ((res = mp_mul(&q, d, &q)) != MP_OKAY) { 
+      goto ERR;
+   }
+   
+   /* a = a + q */
+   if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
+      goto ERR;
+   }
+   
+   if (mp_cmp_mag(a, n) != MP_LT) {
+      s_mp_sub(a, n, a);
+      goto top;
+   }
+   
+ERR:
+   mp_clear(&q);
+   return res;
+}
+
+#endif
+
+/* End: bn_mp_reduce_2k_l.c */
+
 /* Start: bn_mp_reduce_2k_setup.c */
 #include <tommath.h>
 #ifdef BN_MP_REDUCE_2K_SETUP_C
@@ -6423,8 +6487,7 @@ ERR:
  */
 
 /* determines the setup value */
-int 
+int mp_reduce_2k_setup(mp_int *a, mp_digit *d)
-mp_reduce_2k_setup(mp_int *a, mp_digit *d)
 {
    int res, p;
    mp_int tmp;
@@ -6452,6 +6515,50 @@ mp_reduce_2k_setup(mp_int *a, mp_digit *d)
 
 /* End: bn_mp_reduce_2k_setup.c */
 
+/* Start: bn_mp_reduce_2k_setup_l.c */
+#include <tommath.h>
+#ifdef BN_MP_REDUCE_2K_SETUP_L_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 the setup value */
+int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
+{
+   int    res;
+   mp_int tmp;
+   
+   if ((res = mp_init(&tmp)) != MP_OKAY) {
+      return res;
+   }
+   
+   if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
+      goto ERR;
+   }
+   
+   if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
+      goto ERR;
+   }
+   
+ERR:
+   mp_clear(&tmp);
+   return res;
+}
+#endif
+
+/* End: bn_mp_reduce_2k_setup_l.c */
+
 /* Start: bn_mp_reduce_is_2k.c */
 #include <tommath.h>
 #ifdef BN_MP_REDUCE_IS_2K_C
@@ -6477,9 +6584,9 @@ int mp_reduce_is_2k(mp_int *a)
    mp_digit iz;
    
    if (a->used == 0) {
-      return 0;
+      return MP_NO;
    } else if (a->used == 1) {
-      return 1;
+      return MP_YES;
    } else if (a->used > 1) {
       iy = mp_count_bits(a);
       iz = 1;
@@ -6488,7 +6595,7 @@ int mp_reduce_is_2k(mp_int *a)
       /* Test every bit from the second digit up, must be 1 */
       for (ix = DIGIT_BIT; ix < iy; ix++) {
           if ((a->dp[iw] & iz) == 0) {
-             return 0;
+             return MP_NO;
           }
           iz <<= 1;
           if (iz > (mp_digit)MP_MASK) {
@@ -6497,13 +6604,57 @@ int mp_reduce_is_2k(mp_int *a)
           }
       }
    }
-   return 1;
+   return MP_YES;
 }
 
 #endif
 
 /* End: bn_mp_reduce_is_2k.c */
 
+/* Start: bn_mp_reduce_is_2k_l.c */
+#include <tommath.h>
+#ifdef BN_MP_REDUCE_IS_2K_L_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 reduce_2k_l can be used */
+int mp_reduce_is_2k_l(mp_int *a)
+{
+   int ix, iy;
+   
+   if (a->used == 0) {
+      return MP_NO;
+   } else if (a->used == 1) {
+      return MP_YES;
+   } else if (a->used > 1) {
+      /* if more than half of the digits are -1 we're sold */
+      for (iy = ix = 0; ix < a->used; ix++) {
+          if (a->dp[ix] == MP_MASK) {
+              ++iy;
+          }
+      }
+      return (iy >= (a->used/2)) ? MP_YES : MP_NO;
+      
+   }
+   return MP_NO;
+}
+
+#endif
+
+/* End: bn_mp_reduce_is_2k_l.c */
+
 /* Start: bn_mp_reduce_setup.c */
 #include <tommath.h>
 #ifdef BN_MP_REDUCE_SETUP_C
@@ -7138,8 +7289,7 @@ mp_submod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
  */
 
 /* store in signed [big endian] format */
-int
+int mp_to_signed_bin (mp_int * a, unsigned char *b)
-mp_to_signed_bin (mp_int * a, unsigned char *b)
 {
   int     res;
 
@@ -7153,6 +7303,37 @@ mp_to_signed_bin (mp_int * a, unsigned char *b)
 
 /* End: bn_mp_to_signed_bin.c */
 
+/* Start: bn_mp_to_signed_bin_n.c */
+#include <tommath.h>
+#ifdef BN_MP_TO_SIGNED_BIN_N_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
+ */
+
+/* store in signed [big endian] format */
+int mp_to_signed_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen)
+{
+   if (*outlen < (unsigned long)mp_signed_bin_size(a)) {
+      return MP_VAL;
+   }
+   *outlen = mp_signed_bin_size(a);
+   return mp_to_signed_bin(a, b);
+}
+#endif
+
+/* End: bn_mp_to_signed_bin_n.c */
+
 /* Start: bn_mp_to_unsigned_bin.c */
 #include <tommath.h>
 #ifdef BN_MP_TO_UNSIGNED_BIN_C
@@ -7172,8 +7353,7 @@ mp_to_signed_bin (mp_int * a, unsigned char *b)
  */
 
 /* store in unsigned [big endian] format */
-int
+int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
-mp_to_unsigned_bin (mp_int * a, unsigned char *b)
 {
   int     x, res;
   mp_int  t;
@@ -7202,6 +7382,37 @@ mp_to_unsigned_bin (mp_int * a, unsigned char *b)
 
 /* End: bn_mp_to_unsigned_bin.c */
 
+/* Start: bn_mp_to_unsigned_bin_n.c */
+#include <tommath.h>
+#ifdef BN_MP_TO_UNSIGNED_BIN_N_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
+ */
+
+/* store in unsigned [big endian] format */
+int mp_to_unsigned_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen)
+{
+   if (*outlen < (unsigned long)mp_unsigned_bin_size(a)) {
+      return MP_VAL;
+   }
+   *outlen = mp_unsigned_bin_size(a);
+   return mp_to_unsigned_bin(a, b);
+}
+#endif
+
+/* End: bn_mp_to_unsigned_bin_n.c */
+
 /* Start: bn_mp_toom_mul.c */
 #include <tommath.h>
 #ifdef BN_MP_TOOM_MUL_C
@@ -7894,8 +8105,7 @@ int mp_toradix_n(mp_int * a, char *str, int radix, int maxlen)
  */
 
 /* get the size for an unsigned equivalent */
-int
+int mp_unsigned_bin_size (mp_int * a)
-mp_unsigned_bin_size (mp_int * a)
 {
   int     size = mp_count_bits (a);
   return (size / 8 + ((size & 7) != 0 ? 1 : 0));
@@ -8218,11 +8428,12 @@ s_mp_add (mp_int * a, mp_int * b, mp_int * c)
    #define TAB_SIZE 256
 #endif
 
-int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
+int s_mp_exptmod (mp_int * G, mp_int * X, 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*,mp_int*,mp_int*);
 
   /* find window size */
   x = mp_count_bits (X);
@@ -8269,9 +8480,18 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
   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
    *
@@ -8293,11 +8513,14 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
   }
 
   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;
     }
-    if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
+
+    /* reduce modulo P */
+    if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
       goto LBL_MU;
     }
   }
@@ -8309,7 +8532,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
     if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
       goto LBL_MU;
     }
-    if ((err = mp_reduce (&M[x], P, &mu)) != MP_OKAY) {
+    if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
       goto LBL_MU;
     }
   }
@@ -8358,7 +8581,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
         goto LBL_RES;
       }
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
         goto LBL_RES;
       }
       continue;
@@ -8375,7 +8598,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
         if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
           goto LBL_RES;
         }
-        if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+        if ((err = redux (&res, P, &mu)) != MP_OKAY) {
           goto LBL_RES;
         }
       }
@@ -8384,7 +8607,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
       if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
         goto LBL_RES;
       }
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
         goto LBL_RES;
       }
 
@@ -8402,7 +8625,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
       if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
         goto LBL_RES;
       }
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
         goto LBL_RES;
       }
 
@@ -8412,7 +8635,7 @@ int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
         if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
           goto LBL_RES;
         }
-        if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+        if ((err = redux (&res, P, &mu)) != MP_OKAY) {
           goto LBL_RES;
         }
       }
@@ -8803,11 +9026,12 @@ s_mp_sub (mp_int * a, mp_int * b, mp_int * c)
  CPU                    /Compiler     /MUL CUTOFF/SQR CUTOFF
 -------------------------------------------------------------
  Intel P4 Northwood     /GCC v3.4.1   /        88/       128/LTM 0.32 ;-)
+ AMD Athlon64           /GCC v3.4.4   /        74/       124/LTM 0.34
  
 */
 
-int     KARATSUBA_MUL_CUTOFF = 88,      /* Min. number of digits before Karatsuba multiplication is used. */
+int     KARATSUBA_MUL_CUTOFF = 74,      /* Min. number of digits before Karatsuba multiplication is used. */
-        KARATSUBA_SQR_CUTOFF = 128,     /* Min. number of digits before Karatsuba squaring is used. */
+        KARATSUBA_SQR_CUTOFF = 124,     /* Min. number of digits before Karatsuba squaring is used. */
         
         TOOM_MUL_CUTOFF      = 350,      /* no optimal values of these are known yet so set em high */
         TOOM_SQR_CUTOFF      = 400; 

tommath.h

@@ -429,6 +429,15 @@ int mp_reduce_2k_setup(mp_int *a, mp_digit *d);
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
 int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d);
 
+/* returns true if a can be reduced with mp_reduce_2k_l */
+int mp_reduce_is_2k_l(mp_int *a);
+
+/* determines k value for 2k reduction */
+int mp_reduce_2k_setup_l(mp_int *a, mp_int *d);
+
+/* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
+int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d);
+
 /* d = a**b (mod c) */
 int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
 
@@ -511,12 +520,14 @@ int mp_count_bits(mp_int *a);
 int mp_unsigned_bin_size(mp_int *a);
 int mp_read_unsigned_bin(mp_int *a, unsigned char *b, int c);
 int mp_to_unsigned_bin(mp_int *a, unsigned char *b);
+int mp_to_unsigned_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen);
 
 int mp_signed_bin_size(mp_int *a);
 int mp_read_signed_bin(mp_int *a, unsigned char *b, int c);
 int mp_to_signed_bin(mp_int *a, unsigned char *b);
+int mp_to_signed_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen);
 
-int mp_read_radix(mp_int *a, char *str, int radix);
+int mp_read_radix(mp_int *a, const char *str, int radix);
 int mp_toradix(mp_int *a, char *str, int radix);
 int mp_toradix_n(mp_int * a, char *str, int radix, int maxlen);
 int mp_radix_size(mp_int *a, int radix, int *size);
@@ -554,7 +565,7 @@ int fast_mp_invmod(mp_int *a, mp_int *b, mp_int *c);
 int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c);
 int fast_mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
 int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int mode);
-int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y);
+int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int mode);
 void bn_reverse(unsigned char *s, int len);
 
 extern const char *mp_s_rmap;

tommath.src

@@ -66,7 +66,7 @@ QUALCOMM Australia \\
 }
 }
 \maketitle
-This text has been placed in the public domain.  This text corresponds to the v0.30 release of the 
+This text has been placed in the public domain.  This text corresponds to the v0.34 release of the 
 LibTomMath project.
 
 \begin{alltt}

tommath_class.h

@@ -90,8 +90,11 @@
 #define BN_MP_READ_UNSIGNED_BIN_C
 #define BN_MP_REDUCE_C
 #define BN_MP_REDUCE_2K_C
+#define BN_MP_REDUCE_2K_L_C
 #define BN_MP_REDUCE_2K_SETUP_C
+#define BN_MP_REDUCE_2K_SETUP_L_C
 #define BN_MP_REDUCE_IS_2K_C
+#define BN_MP_REDUCE_IS_2K_L_C
 #define BN_MP_REDUCE_SETUP_C
 #define BN_MP_RSHD_C
 #define BN_MP_SET_C
@@ -105,7 +108,9 @@
 #define BN_MP_SUB_D_C
 #define BN_MP_SUBMOD_C
 #define BN_MP_TO_SIGNED_BIN_C
+#define BN_MP_TO_SIGNED_BIN_N_C
 #define BN_MP_TO_UNSIGNED_BIN_C
+#define BN_MP_TO_UNSIGNED_BIN_N_C
 #define BN_MP_TOOM_MUL_C
 #define BN_MP_TOOM_SQR_C
 #define BN_MP_TORADIX_C
@@ -324,11 +329,12 @@
    #define BN_MP_CLEAR_C
    #define BN_MP_ABS_C
    #define BN_MP_CLEAR_MULTI_C
+   #define BN_MP_REDUCE_IS_2K_L_C
+   #define BN_S_MP_EXPTMOD_C
    #define BN_MP_DR_IS_MODULUS_C
    #define BN_MP_REDUCE_IS_2K_C
    #define BN_MP_ISODD_C
    #define BN_MP_EXPTMOD_FAST_C
-   #define BN_S_MP_EXPTMOD_C
 #endif
 
 #if defined(BN_MP_EXPTMOD_FAST_C)
@@ -725,6 +731,17 @@
    #define BN_MP_CLEAR_C
 #endif
 
+#if defined(BN_MP_REDUCE_2K_L_C)
+   #define BN_MP_INIT_C
+   #define BN_MP_COUNT_BITS_C
+   #define BN_MP_DIV_2D_C
+   #define BN_MP_MUL_C
+   #define BN_S_MP_ADD_C
+   #define BN_MP_CMP_MAG_C
+   #define BN_S_MP_SUB_C
+   #define BN_MP_CLEAR_C
+#endif
+
 #if defined(BN_MP_REDUCE_2K_SETUP_C)
    #define BN_MP_INIT_C
    #define BN_MP_COUNT_BITS_C
@@ -733,11 +750,22 @@
    #define BN_S_MP_SUB_C
 #endif
 
+#if defined(BN_MP_REDUCE_2K_SETUP_L_C)
+   #define BN_MP_INIT_C
+   #define BN_MP_2EXPT_C
+   #define BN_MP_COUNT_BITS_C
+   #define BN_S_MP_SUB_C
+   #define BN_MP_CLEAR_C
+#endif
+
 #if defined(BN_MP_REDUCE_IS_2K_C)
    #define BN_MP_REDUCE_2K_C
    #define BN_MP_COUNT_BITS_C
 #endif
 
+#if defined(BN_MP_REDUCE_IS_2K_L_C)
+#endif
+
 #if defined(BN_MP_REDUCE_SETUP_C)
    #define BN_MP_2EXPT_C
    #define BN_MP_DIV_C
@@ -815,6 +843,11 @@
    #define BN_MP_TO_UNSIGNED_BIN_C
 #endif
 
+#if defined(BN_MP_TO_SIGNED_BIN_N_C)
+   #define BN_MP_SIGNED_BIN_SIZE_C
+   #define BN_MP_TO_SIGNED_BIN_C
+#endif
+
 #if defined(BN_MP_TO_UNSIGNED_BIN_C)
    #define BN_MP_INIT_COPY_C
    #define BN_MP_ISZERO_C
@@ -822,6 +855,11 @@
    #define BN_MP_CLEAR_C
 #endif
 
+#if defined(BN_MP_TO_UNSIGNED_BIN_N_C)
+   #define BN_MP_UNSIGNED_BIN_SIZE_C
+   #define BN_MP_TO_UNSIGNED_BIN_C
+#endif
+
 #if defined(BN_MP_TOOM_MUL_C)
    #define BN_MP_INIT_MULTI_C
    #define BN_MP_MOD_2D_C
@@ -902,10 +940,12 @@
    #define BN_MP_INIT_C
    #define BN_MP_CLEAR_C
    #define BN_MP_REDUCE_SETUP_C
+   #define BN_MP_REDUCE_C
+   #define BN_MP_REDUCE_2K_SETUP_L_C
+   #define BN_MP_REDUCE_2K_L_C
    #define BN_MP_MOD_C
    #define BN_MP_COPY_C
    #define BN_MP_SQR_C
-   #define BN_MP_REDUCE_C
    #define BN_MP_MUL_C
    #define BN_MP_SET_C
    #define BN_MP_EXCH_C