added libtommath-0.19

bn.tex

@@ -1,7 +1,7 @@
 \documentclass[]{article}
 \begin{document}
 
-\title{LibTomMath v0.18 \\ A Free Multiple Precision Integer Library \\ http://math.libtomcrypt.org }
+\title{LibTomMath v0.19 \\ A Free Multiple Precision Integer Library \\ http://math.libtomcrypt.org }
 \author{Tom St Denis \\ tomstdenis@iahu.ca}
 \maketitle
 \newpage

bn_mp_exptmod.c

@@ -65,8 +65,8 @@ mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
      dr = mp_reduce_is_2k(P) << 1;
   }
     
-  /* if the modulus is odd use the fast method */
+  /* if the modulus is odd or dr != 0 use the fast method */
-  if ((mp_isodd (P) == 1 || dr !=  0) && P->used > 4) {
+  if (mp_isodd (P) == 1 || dr !=  0) {
     return mp_exptmod_fast (G, X, P, Y, dr);
   } else {
     return s_mp_exptmod (G, X, P, Y);

bn_mp_exptmod_fast.c

@@ -80,6 +80,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
      if (((P->used * 2 + 1) < MP_WARRAY) &&
           P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
         redux = fast_mp_montgomery_reduce;
+
      } else {
         /* use slower baselien method */
         redux = mp_montgomery_reduce;

bn_mp_karatsuba_mul.c

@@ -49,6 +49,7 @@ mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
   mp_int  x0, x1, y0, y1, t1, x0y0, x1y1;
   int     B, err;
 
+  /* default the return code to an error */
   err = MP_MEM;
 
   /* min # of digits */
@@ -149,6 +150,7 @@ mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
   if (mp_add (&t1, &x1y1, c) != MP_OKAY)
     goto X1Y1;          /* t1 = x0y0 + t1 + x1y1 */
 
+  /* Algorithm succeeded set the return code to MP_OKAY */
   err = MP_OKAY;
 
 X1Y1:mp_clear (&x1y1);

bn_mp_montgomery_reduce.c

@@ -76,6 +76,7 @@ mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
   }
 
   /* x = x/b**n.used */
+  mp_clamp(x);
   mp_rshd (x, n->used);
 
   /* if A >= m then A = A - m */

bn_mp_reduce_is_2k.c

@@ -27,7 +27,8 @@ mp_reduce_is_2k(mp_int *a)
    } else if (a->used > 1) {
       iy = mp_count_bits(a);
       for (ix = DIGIT_BIT; ix < iy; ix++) {
-          if ((a->dp[ix/DIGIT_BIT] & ((mp_digit)1 << (mp_digit)(ix % DIGIT_BIT))) == 0) {
+          if ((a->dp[ix/DIGIT_BIT] & 
+              ((mp_digit)1 << (mp_digit)(ix % DIGIT_BIT))) == 0) {
              return 0;
           }
       }

bn_mp_toom_mul.c

@@ -14,7 +14,7 @@
  */
 #include <tommath.h>
 
-/* multiplication using Toom-Cook 3-way algorithm */
+/* multiplication using the Toom-Cook 3-way algorithm */
 int 
 mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
 {
@@ -22,14 +22,16 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
     int res, B;
         
     /* init temps */
-    if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &b0, &b1, &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
+    if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4, 
+                             &a0, &a1, &a2, &b0, &b1, 
+                             &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
        return res;
     }
     
     /* B */
     B = MIN(a->used, b->used) / 3;
     
-    /* a = a2 * B^2 + a1 * B + a0 */
+    /* a = a2 * B**2 + a1 * B + a0 */
     if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
        goto ERR;
     }
@@ -45,7 +47,7 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
     }
     mp_rshd(&a2, B*2);
     
-    /* b = b2 * B^2 + b1 * B + b0 */
+    /* b = b2 * B**2 + b1 * B + b0 */
     if ((res = mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
        goto ERR;
     }
@@ -159,7 +161,8 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
        16 8  4  2  1
        1  0  0  0  0
        
-       using 12 subtractions, 4 shifts, 2 small divisions and 1 small multiplication 
+       using 12 subtractions, 4 shifts, 
+              2 small divisions and 1 small multiplication 
      */
      
      /* r1 - r4 */
@@ -262,7 +265,9 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
      }     
      
 ERR:
-     mp_clear_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &b0, &b1, &b2, &tmp1, &tmp2, NULL);
+     mp_clear_multi(&w0, &w1, &w2, &w3, &w4, 
+                    &a0, &a1, &a2, &b0, &b1, 
+                    &b2, &tmp1, &tmp2, NULL);
      return res;
 }     
      

bn_radix.c

@@ -73,6 +73,14 @@ mp_toradix (mp_int * a, char *str, int radix)
     return MP_VAL;
   }
   
+  /* quick out if its zero */
+  if (mp_iszero(a) == 1) {
+     *str++ = '0';
+     *str = '\0';
+     return MP_OKAY;
+  }
+  
+
   if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
     return res;
   }

bn_s_mp_exptmod.c

@@ -65,21 +65,26 @@ s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 
   /* create M table
    *
-   * The M table contains powers of the input base, e.g. M[x] = G**x mod P
+   * The M table contains powers of the base, 
+   * e.g. M[x] = G**x mod P
    *
-   * The first half of the table is not computed though accept for M[0] and M[1]
+   * The first half of the table is not 
+   * computed though accept for M[0] and M[1]
    */
   if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
     goto __MU;
   }
 
-  /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
+  /* compute the value at M[1<<(winsize-1)] by squaring 
+   * M[1] (winsize-1) times 
+   */
   if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
     goto __MU;
   }
 
   for (x = 0; x < (winsize - 1); x++) {
-    if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
+    if ((err = mp_sqr (&M[1 << (winsize - 1)], 
+                       &M[1 << (winsize - 1)])) != MP_OKAY) {
       goto __MU;
     }
     if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {

bn_s_mp_sqr.c

@@ -24,19 +24,19 @@ s_mp_sqr (mp_int * a, mp_int * b)
   mp_digit u, tmpx, *tmpt;
 
   pa = a->used;
-  if ((res = mp_init_size (&t, pa + pa + 1)) != MP_OKAY) {
+  if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
     return res;
   }
-  t.used = pa + pa + 1;
+  t.used = 2*pa + 1;
 
   for (ix = 0; ix < pa; ix++) {
     /* first calculate the digit at 2*ix */
     /* calculate double precision result */
-    r = ((mp_word) t.dp[ix + ix]) + 
+    r = ((mp_word) t.dp[2*ix]) + 
         ((mp_word) a->dp[ix]) * ((mp_word) a->dp[ix]);
 
     /* store lower part in result */
-    t.dp[ix + ix] = (mp_digit) (r & ((mp_word) MP_MASK));
+    t.dp[2*ix] = (mp_digit) (r & ((mp_word) MP_MASK));
 
     /* get the carry */
     u = (r >> ((mp_word) DIGIT_BIT));
@@ -45,14 +45,14 @@ s_mp_sqr (mp_int * a, mp_int * b)
     tmpx = a->dp[ix];
 
     /* alias for where to store the results */
-    tmpt = t.dp + (ix + ix + 1);
+    tmpt = t.dp + (2*ix + 1);
     
     for (iy = ix + 1; iy < pa; iy++) {
       /* first calculate the product */
       r = ((mp_word) tmpx) * ((mp_word) a->dp[iy]);
 
       /* now calculate the double precision result, note we use
-       * addition instead of *2 since its easier to optimize
+       * addition instead of *2 since it's easier to optimize
        */
       r = ((mp_word) * tmpt) + r + r + ((mp_word) u);
 

changes.txt

@@ -1,3 +1,10 @@
+June 6th, 2003
+v0.19  -- Fixed a bug in mp_montgomery_reduce() which was introduced when I tweaked mp_rshd() in the previous release.
+          Essentially the digits were not trimmed before the compare which cause a subtraction to occur all the time.
+       -- Fixed up etc/tune.c a bit to stop testing new cutoffs after 16 failures [to find more optimal points]. 
+          Brute force ho!
+          
+
 May 29th, 2003
 v0.18  -- Fixed a bug in s_mp_sqr which would handle carries properly just not very elegantly.
           (e.g. correct result, just bad looking code)

etc/drprimes.txt

@@ -0,0 +1,6 @@
+224-bit prime:
+p == 26959946667150639794667015087019630673637144422540572481103341844143
+
+532-bit prime:
+p == 14059105607947488696282932836518693308967803494693489478439861164411992439598399594747002144074658928593502845729752797260025831423419686528151609940203368691747
+

etc/mersenne.c

@@ -57,7 +57,7 @@ is_mersenne (long s, int *pp)
 
   /* if u == 0 then its prime */
   if (mp_iszero (&u) == 1) {
-    mp_prime_is_prime(&n, 3, pp);
+    mp_prime_is_prime(&n, 8, pp);
   if (*pp != 1) printf("FAILURE\n");
   }
 

etc/tune.c

@@ -5,6 +5,12 @@
 #include <tommath.h>
 #include <time.h>
 
+/* how many times todo each size mult.  Depends on your computer.  For slow computers
+ * this can be low like 5 or 10.  For fast [re: Athlon] should be 25 - 50 or so 
+ */
+#define TIMES 50
+
+
 #ifndef X86_TIMER
 
 /* generic ISO C timer */
@@ -31,7 +37,7 @@ time_mult (int max)
   for (x = 32; x <= max; x += 4) {
     mp_rand (&a, x);
     mp_rand (&b, x);
-    for (y = 0; y < 100; y++) {
+    for (y = 0; y < TIMES; y++) {
       mp_mul (&a, &b, &c);
     }
   }
@@ -53,7 +59,7 @@ time_sqr (int max)
   t_start();
   for (x = 32; x <= max; x += 4) {
     mp_rand (&a, x);
-    for (y = 0; y < 100; y++) {
+    for (y = 0; y < TIMES; y++) {
       mp_sqr (&a, &b);
     }
   }
@@ -65,7 +71,7 @@ time_sqr (int max)
 int
 main (void)
 {
-  int     best_kmult, best_tmult, best_ksquare, best_tsquare;
+  int     best_kmult, best_tmult, best_ksquare, best_tsquare, counter;
   ulong64 best, ti;
   FILE   *log;
 
@@ -77,39 +83,52 @@ main (void)
     
   log = fopen ("mult.log", "w");
   best = -1;
+  counter = 16;
   for (KARATSUBA_MUL_CUTOFF = 8; KARATSUBA_MUL_CUTOFF <= 200; KARATSUBA_MUL_CUTOFF++) {
     ti = time_mult (300);
     printf ("%4d : %9llu            \r", KARATSUBA_MUL_CUTOFF, ti);
     fprintf (log, "%d, %llu\n", KARATSUBA_MUL_CUTOFF, ti);
     fflush (stdout);
     if (ti < best) {
-      printf ("New best: %llu, %d         \n", ti, KARATSUBA_MUL_CUTOFF);
+      printf ("New best: %llu, %d         \r", ti, KARATSUBA_MUL_CUTOFF);
       best = ti;
       best_kmult = KARATSUBA_MUL_CUTOFF;
+      counter = 16;
+    } else if (--counter == 0) {
+       printf("No better found in 16 trials.\n");
+       break;
     }
   }
   fclose (log);
+  printf("Karatsuba Multiplier Cutoff (KARATSUBA_MUL_CUTOFF) == %d\n", best_kmult);
   
   /* tune squaring */
   log = fopen ("sqr.log", "w");
   best = -1;
+  counter = 16;
   for (KARATSUBA_SQR_CUTOFF = 8; KARATSUBA_SQR_CUTOFF <= 200; KARATSUBA_SQR_CUTOFF++) {
     ti = time_sqr (300);
     printf ("%4d : %9llu             \r", KARATSUBA_SQR_CUTOFF, ti);
     fprintf (log, "%d, %llu\n", KARATSUBA_SQR_CUTOFF, ti);
     fflush (stdout);
     if (ti < best) {
-      printf ("New best: %llu, %d         \n", ti, KARATSUBA_SQR_CUTOFF);
+      printf ("New best: %llu, %d         \r", ti, KARATSUBA_SQR_CUTOFF);
       best = ti;
       best_ksquare = KARATSUBA_SQR_CUTOFF;
+      counter = 16;
+    } else if (--counter == 0) {
+       printf("No better found in 16 trials.\n");
+       break;
     }
   }
   fclose (log);
+  printf("Karatsuba Squaring Cutoff (KARATSUBA_SQR_CUTOFF) == %d\n", best_ksquare);
   
   KARATSUBA_MUL_CUTOFF = best_kmult;
   KARATSUBA_SQR_CUTOFF = best_ksquare;
     
   /* tune TOOM mult */
+  counter = 16;
   log = fopen ("tmult.log", "w");
   best = -1;
   for (TOOM_MUL_CUTOFF = best_kmult*5; TOOM_MUL_CUTOFF <= 800; TOOM_MUL_CUTOFF++) {
@@ -118,32 +137,40 @@ main (void)
     fprintf (log, "%d, %llu\n", TOOM_MUL_CUTOFF, ti);
     fflush (stdout);
     if (ti < best) {
-      printf ("New best: %llu, %d         \n", ti, TOOM_MUL_CUTOFF);
+      printf ("New best: %llu, %d         \r", ti, TOOM_MUL_CUTOFF);
       best = ti;
       best_tmult = TOOM_MUL_CUTOFF;
+      counter = 16;
+    } else if (--counter == 0) {
+       printf("No better found in 16 trials.\n");
+       break;
     }
   }
   fclose (log);   
+  printf("Toom-Cook Multiplier Cutoff (TOOM_MUL_CUTOFF) == %d\n", best_tmult);
   
   /* tune TOOM sqr */
   log = fopen ("tsqr.log", "w");
   best = -1;
+  counter = 16;
   for (TOOM_SQR_CUTOFF = best_ksquare*3; TOOM_SQR_CUTOFF <= 800; TOOM_SQR_CUTOFF++) {
     ti = time_sqr (1200);
     printf ("%4d : %9llu           \r", TOOM_SQR_CUTOFF, ti);
     fprintf (log, "%d, %llu\n", TOOM_SQR_CUTOFF, ti);
     fflush (stdout);
     if (ti < best) {
-      printf ("New best: %llu, %d         \n", ti, TOOM_SQR_CUTOFF);
+      printf ("New best: %llu, %d         \r", ti, TOOM_SQR_CUTOFF);
       best = ti;
       best_tsquare = TOOM_SQR_CUTOFF;
+      counter = 16;
+    } else if (--counter == 0) {
+       printf("No better found in 16 trials.\n");
+       break;
     }
   }
   fclose (log);   
+  printf("Toom-Cook Squaring Cutoff (TOOM_SQR_CUTOFF) == %d\n", best_tsquare);
 
-  printf
-    ("\n\n\nKaratsuba Multiplier Cutoff: %d\nKaratsuba Squaring Cutoff: %d\nToom Multiplier Cutoff: %d\nToom Squaring Cutoff: %d\n",
-     best_kmult, best_ksquare, best_tmult, best_tsquare);
 
   return 0;
 }

makefile

@@ -1,6 +1,6 @@
 CFLAGS  +=  -I./ -Wall -W -Wshadow -O3 -fomit-frame-pointer -funroll-loops
 
-VERSION=0.18
+VERSION=0.19
 
 default: libtommath.a
 

makefile.bcc

@@ -32,6 +32,6 @@ TARGET = libtommath.lib
 
 $(TARGET): $(OBJECTS)
 
-.c.objbj:
+.c.obj:
 	$(CC) $(CFLAGS) $<
 	$(LIB) $(TARGET) -+$@

poster.tex

@@ -1,14 +1,15 @@
 \documentclass[landscape,11pt]{article}
 \usepackage{amsmath, amssymb}
+\usepackage{hyperref}
 \begin{document}
 
 \hspace*{-3in}
 \begin{tabular}{llllll}
-$c = a + b$  & {\tt mp\_add(\&a, \&b, \&c)} & $b = 2a$  & {\tt mp\_mul\_2(\&a, \&b)} & Greater Than & MP\_GT \\
+$c = a + b$  & {\tt mp\_add(\&a, \&b, \&c)} & $b = 2a$  & {\tt mp\_mul\_2(\&a, \&b)} & \\
-$c = a - b$  & {\tt mp\_sub(\&a, \&b, \&c)} & $b = a/2$ & {\tt mp\_div\_2(\&a, \&b)} & Equal To & MP\_EQ \\
+$c = a - b$  & {\tt mp\_sub(\&a, \&b, \&c)} & $b = a/2$ & {\tt mp\_div\_2(\&a, \&b)} & \\
-$c = ab $   & {\tt mp\_mul(\&a, \&b, \&c)}  & $c = 2^ba$  & {\tt mp\_mul\_2d(\&a, b, \&c)} & Less Than & MP\_LT \\
+$c = ab $   & {\tt mp\_mul(\&a, \&b, \&c)}  & $c = 2^ba$  & {\tt mp\_mul\_2d(\&a, b, \&c)}  \\
 $b = a^2 $  & {\tt mp\_sqr(\&a, \&b)}       & $c = a/2^b, d = a \mod 2^b$ & {\tt mp\_div\_2d(\&a, b, \&c, \&d)} \\
-$c = \lfloor a/b \rfloor, d = a \mod b$ & {\tt mp\_div(\&a, \&b, \&c, \&d)} & $c = a \mod 2^b $  & {\tt mp\_mod\_2d(\&a, b, \&c)}  & Bits per digit & DIGIT\_BIT \\
+$c = \lfloor a/b \rfloor, d = a \mod b$ & {\tt mp\_div(\&a, \&b, \&c, \&d)} & $c = a \mod 2^b $  & {\tt mp\_mod\_2d(\&a, b, \&c)}  \\
  && \\
 $a = b $  & {\tt mp\_set\_int(\&a, b)}  & $c = a \vee b$  & {\tt mp\_or(\&a, \&b, \&c)}  \\
 $b = a $  & {\tt mp\_copy(\&a, \&b)} & $c = a \wedge b$  & {\tt mp\_and(\&a, \&b, \&c)}  \\
@@ -28,5 +29,8 @@ $a \leftarrow buf[0..len-1]$          & {\tt mp\_read\_unsigned\_bin(\&a, buf, l
 &\\
 $b = \sqrt{a}$ & {\tt mp\_sqrt(\&a, \&b)}  & $c = \mbox{gcd}(a, b)$ & {\tt mp\_gcd(\&a, \&b, \&c)} \\
 $c = a^{1/b}$ & {\tt mp\_n\_root(\&a, b, \&c)} & $c = \mbox{lcm}(a, b)$ & {\tt mp\_lcm(\&a, \&b, \&c)} \\
+&\\
+Greater Than & MP\_GT & Equal To & MP\_EQ \\
+Less Than & MP\_LT & Bits per digit & DIGIT\_BIT \\
 \end{tabular}
 \end{document}

pre_gen/mpi.c

@@ -2060,8 +2060,8 @@ mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
      dr = mp_reduce_is_2k(P) << 1;
   }
     
-  /* if the modulus is odd use the fast method */
+  /* if the modulus is odd or dr != 0 use the fast method */
-  if ((mp_isodd (P) == 1 || dr !=  0) && P->used > 4) {
+  if (mp_isodd (P) == 1 || dr !=  0) {
     return mp_exptmod_fast (G, X, P, Y, dr);
   } else {
     return s_mp_exptmod (G, X, P, Y);
@@ -2155,6 +2155,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
      if (((P->used * 2 + 1) < MP_WARRAY) &&
           P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
         redux = fast_mp_montgomery_reduce;
+
      } else {
         /* use slower baselien method */
         redux = mp_montgomery_reduce;
@@ -2964,6 +2965,7 @@ mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
   mp_int  x0, x1, y0, y1, t1, x0y0, x1y1;
   int     B, err;
 
+  /* default the return code to an error */
   err = MP_MEM;
 
   /* min # of digits */
@@ -3064,6 +3066,7 @@ mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
   if (mp_add (&t1, &x1y1, c) != MP_OKAY)
     goto X1Y1;          /* t1 = x0y0 + t1 + x1y1 */
 
+  /* Algorithm succeeded set the return code to MP_OKAY */
   err = MP_OKAY;
 
 X1Y1:mp_clear (&x1y1);
@@ -3581,6 +3584,7 @@ mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
   }
 
   /* x = x/b**n.used */
+  mp_clamp(x);
   mp_rshd (x, n->used);
 
   /* if A >= m then A = A - m */
@@ -4969,7 +4973,8 @@ mp_reduce_is_2k(mp_int *a)
    } else if (a->used > 1) {
       iy = mp_count_bits(a);
       for (ix = DIGIT_BIT; ix < iy; ix++) {
-          if ((a->dp[ix/DIGIT_BIT] & ((mp_digit)1 << (mp_digit)(ix % DIGIT_BIT))) == 0) {
+          if ((a->dp[ix/DIGIT_BIT] & 
+              ((mp_digit)1 << (mp_digit)(ix % DIGIT_BIT))) == 0) {
              return 0;
           }
       }
@@ -5543,7 +5548,7 @@ mp_to_unsigned_bin (mp_int * a, unsigned char *b)
  */
 #include <tommath.h>
 
-/* multiplication using Toom-Cook 3-way algorithm */
+/* multiplication using the Toom-Cook 3-way algorithm */
 int 
 mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
 {
@@ -5551,14 +5556,16 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
     int res, B;
         
     /* init temps */
-    if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &b0, &b1, &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
+    if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4, 
+                             &a0, &a1, &a2, &b0, &b1, 
+                             &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
        return res;
     }
     
     /* B */
     B = MIN(a->used, b->used) / 3;
     
-    /* a = a2 * B^2 + a1 * B + a0 */
+    /* a = a2 * B**2 + a1 * B + a0 */
     if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
        goto ERR;
     }
@@ -5574,7 +5581,7 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
     }
     mp_rshd(&a2, B*2);
     
-    /* b = b2 * B^2 + b1 * B + b0 */
+    /* b = b2 * B**2 + b1 * B + b0 */
     if ((res = mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
        goto ERR;
     }
@@ -5688,7 +5695,8 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
        16 8  4  2  1
        1  0  0  0  0
        
-       using 12 subtractions, 4 shifts, 2 small divisions and 1 small multiplication 
+       using 12 subtractions, 4 shifts, 
+              2 small divisions and 1 small multiplication 
      */
      
      /* r1 - r4 */
@@ -5791,7 +5799,9 @@ mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
      }     
      
 ERR:
-     mp_clear_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &b0, &b1, &b2, &tmp1, &tmp2, NULL);
+     mp_clear_multi(&w0, &w1, &w2, &w3, &w4, 
+                    &a0, &a1, &a2, &b0, &b1, 
+                    &b2, &tmp1, &tmp2, NULL);
      return res;
 }     
      
@@ -6267,6 +6277,14 @@ mp_toradix (mp_int * a, char *str, int radix)
     return MP_VAL;
   }
   
+  /* quick out if its zero */
+  if (mp_iszero(a) == 1) {
+     *str++ = '0';
+     *str = '\0';
+     return MP_OKAY;
+  }
+  
+
   if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
     return res;
   }
@@ -6624,21 +6642,26 @@ s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 
   /* create M table
    *
-   * The M table contains powers of the input base, e.g. M[x] = G**x mod P
+   * The M table contains powers of the base, 
+   * e.g. M[x] = G**x mod P
    *
-   * The first half of the table is not computed though accept for M[0] and M[1]
+   * The first half of the table is not 
+   * computed though accept for M[0] and M[1]
    */
   if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
     goto __MU;
   }
 
-  /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
+  /* compute the value at M[1<<(winsize-1)] by squaring 
+   * M[1] (winsize-1) times 
+   */
   if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
     goto __MU;
   }
 
   for (x = 0; x < (winsize - 1); x++) {
-    if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
+    if ((err = mp_sqr (&M[1 << (winsize - 1)], 
+                       &M[1 << (winsize - 1)])) != MP_OKAY) {
       goto __MU;
     }
     if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
@@ -6966,19 +6989,19 @@ s_mp_sqr (mp_int * a, mp_int * b)
   mp_digit u, tmpx, *tmpt;
 
   pa = a->used;
-  if ((res = mp_init_size (&t, pa + pa + 1)) != MP_OKAY) {
+  if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
     return res;
   }
-  t.used = pa + pa + 1;
+  t.used = 2*pa + 1;
 
   for (ix = 0; ix < pa; ix++) {
     /* first calculate the digit at 2*ix */
     /* calculate double precision result */
-    r = ((mp_word) t.dp[ix + ix]) + 
+    r = ((mp_word) t.dp[2*ix]) + 
         ((mp_word) a->dp[ix]) * ((mp_word) a->dp[ix]);
 
     /* store lower part in result */
-    t.dp[ix + ix] = (mp_digit) (r & ((mp_word) MP_MASK));
+    t.dp[2*ix] = (mp_digit) (r & ((mp_word) MP_MASK));
 
     /* get the carry */
     u = (r >> ((mp_word) DIGIT_BIT));
@@ -6987,14 +7010,14 @@ s_mp_sqr (mp_int * a, mp_int * b)
     tmpx = a->dp[ix];
 
     /* alias for where to store the results */
-    tmpt = t.dp + (ix + ix + 1);
+    tmpt = t.dp + (2*ix + 1);
     
     for (iy = ix + 1; iy < pa; iy++) {
       /* first calculate the product */
       r = ((mp_word) tmpx) * ((mp_word) a->dp[iy]);
 
       /* now calculate the double precision result, note we use
-       * addition instead of *2 since its easier to optimize
+       * addition instead of *2 since it's easier to optimize
        */
       r = ((mp_word) * tmpt) + r + r + ((mp_word) u);
 

tommath.out

@@ -0,0 +1,143 @@
+\BOOKMARK [0][-]{chapter.1}{Introduction}{}
+\BOOKMARK [1][-]{section.1.1}{Multiple Precision Arithmetic}{chapter.1}
+\BOOKMARK [2][-]{subsection.1.1.1}{The Need for Multiple Precision Arithmetic}{section.1.1}
+\BOOKMARK [2][-]{subsection.1.1.2}{Multiple Precision Arithmetic}{section.1.1}
+\BOOKMARK [2][-]{subsection.1.1.3}{Benefits of Multiple Precision Arithmetic}{section.1.1}
+\BOOKMARK [2][-]{subsection.1.1.4}{Basis of Operations}{section.1.1}
+\BOOKMARK [1][-]{section.1.2}{Purpose of This Text}{chapter.1}
+\BOOKMARK [1][-]{section.1.3}{Discussion and Notation}{chapter.1}
+\BOOKMARK [2][-]{subsection.1.3.1}{Notation}{section.1.3}
+\BOOKMARK [2][-]{subsection.1.3.2}{Work Effort}{section.1.3}
+\BOOKMARK [1][-]{section.1.4}{Exercises}{chapter.1}
+\BOOKMARK [0][-]{chapter.2}{Introduction to LibTomMath}{}
+\BOOKMARK [1][-]{section.2.1}{What is LibTomMath?}{chapter.2}
+\BOOKMARK [1][-]{section.2.2}{Goals of LibTomMath}{chapter.2}
+\BOOKMARK [1][-]{section.2.3}{Choice of LibTomMath}{chapter.2}
+\BOOKMARK [2][-]{subsection.2.3.1}{Code Base}{section.2.3}
+\BOOKMARK [2][-]{subsection.2.3.2}{API Simplicity}{section.2.3}
+\BOOKMARK [2][-]{subsection.2.3.3}{Optimizations}{section.2.3}
+\BOOKMARK [2][-]{subsection.2.3.4}{Portability and Stability}{section.2.3}
+\BOOKMARK [2][-]{subsection.2.3.5}{Choice}{section.2.3}
+\BOOKMARK [0][-]{chapter.3}{Getting Started}{}
+\BOOKMARK [1][-]{section.3.1}{Library Basics}{chapter.3}
+\BOOKMARK [1][-]{section.3.2}{What is a Multiple Precision Integer?}{chapter.3}
+\BOOKMARK [2][-]{subsection.3.2.1}{The mp\137int structure}{section.3.2}
+\BOOKMARK [1][-]{section.3.3}{Argument Passing}{chapter.3}
+\BOOKMARK [1][-]{section.3.4}{Return Values}{chapter.3}
+\BOOKMARK [1][-]{section.3.5}{Initialization and Clearing}{chapter.3}
+\BOOKMARK [2][-]{subsection.3.5.1}{Initializing an mp\137int}{section.3.5}
+\BOOKMARK [2][-]{subsection.3.5.2}{Clearing an mp\137int}{section.3.5}
+\BOOKMARK [1][-]{section.3.6}{Other Initialization Routines}{chapter.3}
+\BOOKMARK [2][-]{subsection.3.6.1}{Initializing Variable Sized mp\137int Structures}{section.3.6}
+\BOOKMARK [2][-]{subsection.3.6.2}{Creating a Clone}{section.3.6}
+\BOOKMARK [2][-]{subsection.3.6.3}{Multiple Integer Initializations And Clearings}{section.3.6}
+\BOOKMARK [1][-]{section.3.7}{Maintenance}{chapter.3}
+\BOOKMARK [2][-]{subsection.3.7.1}{Augmenting Integer Precision}{section.3.7}
+\BOOKMARK [2][-]{subsection.3.7.2}{Clamping Excess Digits}{section.3.7}
+\BOOKMARK [0][-]{chapter.4}{Basic Operations}{}
+\BOOKMARK [1][-]{section.4.1}{Copying an Integer}{chapter.4}
+\BOOKMARK [1][-]{section.4.2}{Zeroing an Integer}{chapter.4}
+\BOOKMARK [1][-]{section.4.3}{Sign Manipulation}{chapter.4}
+\BOOKMARK [2][-]{subsection.4.3.1}{Absolute Value}{section.4.3}
+\BOOKMARK [2][-]{subsection.4.3.2}{Integer Negation}{section.4.3}
+\BOOKMARK [1][-]{section.4.4}{Small Constants}{chapter.4}
+\BOOKMARK [2][-]{subsection.4.4.1}{Setting Small Constants}{section.4.4}
+\BOOKMARK [2][-]{subsection.4.4.2}{Setting Large Constants}{section.4.4}
+\BOOKMARK [1][-]{section.4.5}{Comparisons}{chapter.4}
+\BOOKMARK [2][-]{subsection.4.5.1}{Unsigned Comparisions}{section.4.5}
+\BOOKMARK [2][-]{subsection.4.5.2}{Signed Comparisons}{section.4.5}
+\BOOKMARK [0][-]{chapter.5}{Basic Arithmetic}{}
+\BOOKMARK [1][-]{section.5.1}{Building Blocks}{chapter.5}
+\BOOKMARK [1][-]{section.5.2}{Addition and Subtraction}{chapter.5}
+\BOOKMARK [2][-]{subsection.5.2.1}{Low Level Addition}{section.5.2}
+\BOOKMARK [2][-]{subsection.5.2.2}{Low Level Subtraction}{section.5.2}
+\BOOKMARK [2][-]{subsection.5.2.3}{High Level Addition}{section.5.2}
+\BOOKMARK [2][-]{subsection.5.2.4}{High Level Subtraction}{section.5.2}
+\BOOKMARK [1][-]{section.5.3}{Bit and Digit Shifting}{chapter.5}
+\BOOKMARK [2][-]{subsection.5.3.1}{Multiplication by Two}{section.5.3}
+\BOOKMARK [2][-]{subsection.5.3.2}{Division by Two}{section.5.3}
+\BOOKMARK [1][-]{section.5.4}{Polynomial Basis Operations}{chapter.5}
+\BOOKMARK [2][-]{subsection.5.4.1}{Multiplication by x}{section.5.4}
+\BOOKMARK [2][-]{subsection.5.4.2}{Division by x}{section.5.4}
+\BOOKMARK [1][-]{section.5.5}{Powers of Two}{chapter.5}
+\BOOKMARK [2][-]{subsection.5.5.1}{Multiplication by Power of Two}{section.5.5}
+\BOOKMARK [2][-]{subsection.5.5.2}{Division by Power of Two}{section.5.5}
+\BOOKMARK [2][-]{subsection.5.5.3}{Remainder of Division by Power of Two}{section.5.5}
+\BOOKMARK [0][-]{chapter.6}{Multiplication and Squaring}{}
+\BOOKMARK [1][-]{section.6.1}{The Multipliers}{chapter.6}
+\BOOKMARK [1][-]{section.6.2}{Multiplication}{chapter.6}
+\BOOKMARK [2][-]{subsection.6.2.1}{The Baseline Multiplication}{section.6.2}
+\BOOKMARK [2][-]{subsection.6.2.2}{Faster Multiplication by the ``Comba'' Method}{section.6.2}
+\BOOKMARK [2][-]{subsection.6.2.3}{Polynomial Basis Multiplication}{section.6.2}
+\BOOKMARK [2][-]{subsection.6.2.4}{Karatsuba Multiplication}{section.6.2}
+\BOOKMARK [2][-]{subsection.6.2.5}{Toom-Cook 3-Way Multiplication}{section.6.2}
+\BOOKMARK [2][-]{subsection.6.2.6}{Signed Multiplication}{section.6.2}
+\BOOKMARK [1][-]{section.6.3}{Squaring}{chapter.6}
+\BOOKMARK [2][-]{subsection.6.3.1}{The Baseline Squaring Algorithm}{section.6.3}
+\BOOKMARK [2][-]{subsection.6.3.2}{Faster Squaring by the ``Comba'' Method}{section.6.3}
+\BOOKMARK [2][-]{subsection.6.3.3}{Polynomial Basis Squaring}{section.6.3}
+\BOOKMARK [2][-]{subsection.6.3.4}{Karatsuba Squaring}{section.6.3}
+\BOOKMARK [2][-]{subsection.6.3.5}{Toom-Cook Squaring}{section.6.3}
+\BOOKMARK [2][-]{subsection.6.3.6}{High Level Squaring}{section.6.3}
+\BOOKMARK [0][-]{chapter.7}{Modular Reduction}{}
+\BOOKMARK [1][-]{section.7.1}{Basics of Modular Reduction}{chapter.7}
+\BOOKMARK [1][-]{section.7.2}{The Barrett Reduction}{chapter.7}
+\BOOKMARK [2][-]{subsection.7.2.1}{Fixed Point Arithmetic}{section.7.2}
+\BOOKMARK [2][-]{subsection.7.2.2}{Choosing a Radix Point}{section.7.2}
+\BOOKMARK [2][-]{subsection.7.2.3}{Trimming the Quotient}{section.7.2}
+\BOOKMARK [2][-]{subsection.7.2.4}{Trimming the Residue}{section.7.2}
+\BOOKMARK [2][-]{subsection.7.2.5}{The Barrett Algorithm}{section.7.2}
+\BOOKMARK [2][-]{subsection.7.2.6}{The Barrett Setup Algorithm}{section.7.2}
+\BOOKMARK [1][-]{section.7.3}{The Montgomery Reduction}{chapter.7}
+\BOOKMARK [2][-]{subsection.7.3.1}{Digit Based Montgomery Reduction}{section.7.3}
+\BOOKMARK [2][-]{subsection.7.3.2}{Baseline Montgomery Reduction}{section.7.3}
+\BOOKMARK [2][-]{subsection.7.3.3}{Faster ``Comba'' Montgomery Reduction}{section.7.3}
+\BOOKMARK [2][-]{subsection.7.3.4}{Montgomery Setup}{section.7.3}
+\BOOKMARK [1][-]{section.7.4}{The Diminished Radix Algorithm}{chapter.7}
+\BOOKMARK [2][-]{subsection.7.4.1}{Choice of Moduli}{section.7.4}
+\BOOKMARK [2][-]{subsection.7.4.2}{Choice of k}{section.7.4}
+\BOOKMARK [2][-]{subsection.7.4.3}{Restricted Diminished Radix Reduction}{section.7.4}
+\BOOKMARK [2][-]{subsection.7.4.4}{Unrestricted Diminished Radix Reduction}{section.7.4}
+\BOOKMARK [1][-]{section.7.5}{Algorithm Comparison}{chapter.7}
+\BOOKMARK [0][-]{chapter.8}{Exponentiation}{}
+\BOOKMARK [1][-]{section.8.1}{Exponentiation Basics}{chapter.8}
+\BOOKMARK [2][-]{subsection.8.1.1}{Single Digit Exponentiation}{section.8.1}
+\BOOKMARK [1][-]{section.8.2}{k-ary Exponentiation}{chapter.8}
+\BOOKMARK [2][-]{subsection.8.2.1}{Optimal Values of k}{section.8.2}
+\BOOKMARK [2][-]{subsection.8.2.2}{Sliding-Window Exponentiation}{section.8.2}
+\BOOKMARK [1][-]{section.8.3}{Modular Exponentiation}{chapter.8}
+\BOOKMARK [2][-]{subsection.8.3.1}{Barrett Modular Exponentiation}{section.8.3}
+\BOOKMARK [1][-]{section.8.4}{Quick Power of Two}{chapter.8}
+\BOOKMARK [0][-]{chapter.9}{Higher Level Algorithms}{}
+\BOOKMARK [1][-]{section.9.1}{Integer Division with Remainder}{chapter.9}
+\BOOKMARK [1][-]{section.9.2}{Single Digit Helpers}{chapter.9}
+\BOOKMARK [2][-]{subsection.9.2.1}{Single Digit Addition}{section.9.2}
+\BOOKMARK [2][-]{subsection.9.2.2}{Single Digit Subtraction}{section.9.2}
+\BOOKMARK [2][-]{subsection.9.2.3}{Single Digit Multiplication}{section.9.2}
+\BOOKMARK [2][-]{subsection.9.2.4}{Single Digit Division}{section.9.2}
+\BOOKMARK [2][-]{subsection.9.2.5}{Single Digit Modulo}{section.9.2}
+\BOOKMARK [2][-]{subsection.9.2.6}{Single Digit Root Extraction}{section.9.2}
+\BOOKMARK [1][-]{section.9.3}{Random Number Generation}{chapter.9}
+\BOOKMARK [1][-]{section.9.4}{Formatted Output}{chapter.9}
+\BOOKMARK [2][-]{subsection.9.4.1}{Getting The Output Size}{section.9.4}
+\BOOKMARK [2][-]{subsection.9.4.2}{Generating Radix-n Output}{section.9.4}
+\BOOKMARK [2][-]{subsection.9.4.3}{Reading Radix-n Input}{section.9.4}
+\BOOKMARK [1][-]{section.9.5}{Unformatted Output}{chapter.9}
+\BOOKMARK [2][-]{subsection.9.5.1}{Getting The Output Size}{section.9.5}
+\BOOKMARK [2][-]{subsection.9.5.2}{Generating Output}{section.9.5}
+\BOOKMARK [2][-]{subsection.9.5.3}{Reading Input}{section.9.5}
+\BOOKMARK [0][-]{chapter.10}{Number Theoretic Algorithms}{}
+\BOOKMARK [1][-]{section.10.1}{Greatest Common Divisor}{chapter.10}
+\BOOKMARK [1][-]{section.10.2}{Least Common Multiple}{chapter.10}
+\BOOKMARK [1][-]{section.10.3}{Jacobi Symbol Computation}{chapter.10}
+\BOOKMARK [1][-]{section.10.4}{Modular Inverse}{chapter.10}
+\BOOKMARK [2][-]{subsection.10.4.1}{General Case}{section.10.4}
+\BOOKMARK [2][-]{subsection.10.4.2}{Odd Moduli}{section.10.4}
+\BOOKMARK [1][-]{section.10.5}{Primality Tests}{chapter.10}
+\BOOKMARK [2][-]{subsection.10.5.1}{Trial Division}{section.10.5}
+\BOOKMARK [2][-]{subsection.10.5.2}{The Fermat Test}{section.10.5}
+\BOOKMARK [2][-]{subsection.10.5.3}{The Miller-Rabin Test}{section.10.5}
+\BOOKMARK [2][-]{subsection.10.5.4}{Primality Test in a Bottle}{section.10.5}
+\BOOKMARK [2][-]{subsection.10.5.5}{The Next Prime}{section.10.5}
+\BOOKMARK [1][-]{section.10.6}{Root Extraction}{chapter.10}
+\BOOKMARK [0][-]{appendix*.16}{Appendix}{}