added libtommath-0.21

bn.tex

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

bn_fast_mp_montgomery_reduce.c

@@ -124,7 +124,7 @@ fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
     _W = W + n->used;
 
     for (ix = 0; ix < n->used + 1; ix++) {
-      *tmpx++ = *_W++ & ((mp_word) MP_MASK);
+      *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
     }
 
     /* zero oldused digits, if the input a was larger than

bn_mp_div.c

@@ -14,22 +14,25 @@
  */
 #include <tommath.h>
 
-/* integer signed division. c*b + d == a [e.g. a/b, c=quotient, d=remainder]
+/* integer signed division. 
+ * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
  * HAC pp.598 Algorithm 14.20
  *
- * Note that the description in HAC is horribly incomplete.  For example,
- * it doesn't consider the case where digits are removed from 'x' in the inner
- * loop.  It also doesn't consider the case that y has fewer than three digits, etc..
+ * Note that the description in HAC is horribly 
+ * incomplete.  For example, it doesn't consider 
+ * the case where digits are removed from 'x' in 
+ * the inner loop.  It also doesn't consider the 
+ * case that y has fewer than three digits, etc..
  *
- * The overall algorithm is as described as 14.20 from HAC but fixed to treat these cases.
+ * The overall algorithm is as described as 
+ * 14.20 from HAC but fixed to treat these cases.
 */
 int
 mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
 {
   mp_int  q, x, y, t1, t2;
   int     res, n, t, i, norm, neg;
 
-
   /* is divisor zero ? */
   if (mp_iszero (b) == 1) {
     return MP_VAL;
@@ -73,7 +76,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
   neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
   x.sign = y.sign = MP_ZPOS;
 
-  /* normalize both x and y, ensure that y >= b/2, [b == 2^DIGIT_BIT] */
+  /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
   norm = mp_count_bits(&y) % DIGIT_BIT;
   if (norm < (int)(DIGIT_BIT-1)) {
      norm = (DIGIT_BIT-1) - norm;
@@ -91,8 +94,8 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
   n = x.used - 1;
   t = y.used - 1;
 
-  /* step 2. while (x >= y*b^n-t) do { q[n-t] += 1; x -= y*b^{n-t} } */
-  if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b^{n-t} */
+  /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
+  if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
     goto __Y;
   }
 
@@ -111,7 +114,8 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
     if (i > x.used)
       continue;
 
-    /* step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
+    /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 
+     * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
     if (x.dp[i] == y.dp[t]) {
       q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
     } else {
@@ -124,7 +128,11 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
       q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
     }
 
-    /* step 3.2 while (q{i-t-1} * (yt * b + y{t-1})) > xi * b^2 + xi-1 * b + xi-2 do q{i-t-1} -= 1; */
+    /* while (q{i-t-1} * (yt * b + y{t-1})) > 
+             xi * b**2 + xi-1 * b + xi-2 
+     
+       do q{i-t-1} -= 1; 
+    */
     q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
     do {
       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
@@ -145,7 +153,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
       t2.used = 3;
     } while (mp_cmp_mag(&t1, &t2) == MP_GT);
 
-    /* step 3.3 x = x - q{i-t-1} * y * b^{i-t-1} */
+    /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
     if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
       goto __Y;
     }
@@ -158,7 +166,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
       goto __Y;
     }
 
-    /* step 3.4 if x < 0 then { x = x + y*b^{i-t-1}; q{i-t-1} -= 1; } */
+    /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
     if (x.sign == MP_NEG) {
       if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
         goto __Y;
@@ -174,7 +182,10 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
     }
   }
 
-  /* now q is the quotient and x is the remainder [which we have to normalize] */
+  /* now q is the quotient and x is the remainder 
+   * [which we have to normalize] 
+   */
+
   /* get sign before writing to c */
   x.sign = a->sign;
 

bn_mp_div_3.c

@@ -46,11 +46,11 @@ mp_div_3 (mp_int * a, mp_int *c, mp_digit * d)
       } else {
         t = 0;
       }
-      q.dp[ix] = t;
+      q.dp[ix] = (mp_digit)t;
   }
 
   if (d != NULL) {
-     *d = w;
+     *d = (mp_digit)w;
   }
 
   if (c != NULL) {

bn_mp_div_d.c

@@ -19,7 +19,8 @@ int
 mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
 {
   mp_int  q;
-  mp_word w, t;
+  mp_word w;
+  mp_digit t;
   int     res, ix;
 
   if (b == 0) {
@@ -41,16 +42,16 @@ mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
      w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);
 
      if (w >= b) {
-        t = w / b;
+        t = (mp_digit)(w / b);
         w = w % b;
       } else {
         t = 0;
       }
-      q.dp[ix] = t;
+      q.dp[ix] = (mp_digit)t;
   }
 
   if (d != NULL) {
-     *d = w;
+     *d = (mp_digit)w;
   }
 
   if (c != NULL) {

bn_mp_dr_reduce.c

@@ -60,8 +60,8 @@ mp_dr_reduce (mp_int * x, mp_int * n, mp_digit k)
   /* compute (x mod B**m) + mp * [x/B**m] inline and inplace */
   for (i = 0; i < m; i++) {
       r         = ((mp_word)*tmpx2++) * ((mp_word)k) + *tmpx1 + mu;
-      *tmpx1++  = r & MP_MASK;
-      mu        = r >> ((mp_word)DIGIT_BIT);
+      *tmpx1++  = (mp_digit)(r & MP_MASK);
+      mu        = (mp_digit)(r >> ((mp_word)DIGIT_BIT));
   }
 
   /* set final carry */

bn_mp_montgomery_reduce.c

@@ -61,10 +61,10 @@ mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
 
       /* Multiply and add in place */
       for (iy = 0; iy < n->used; iy++) {
-        r = ((mp_word) mu) * ((mp_word) * tmpn++) + 
-            ((mp_word) u) + ((mp_word) * tmpx);
-        u = (r >> ((mp_word) DIGIT_BIT));
-        *tmpx++ = (r & ((mp_word) MP_MASK));
+        r       = ((mp_word) mu) * ((mp_word) * tmpn++) + 
+                  ((mp_word) u) + ((mp_word) * tmpx);
+        u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
+        *tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
       }
       /* propagate carries */
       while (u) {

bn_mp_mul_d.c

@@ -33,6 +33,7 @@ mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
 
   /* set the new temporary used count */
   c->used = pa + 1;
+  c->sign = a->sign;
 
   {
     register mp_digit u, *tmpa, *tmpc;

bn_mp_n_root.c

@@ -16,11 +16,13 @@
 
 /* find the n'th root of an integer 
  *
- * Result found such that (c)^b <= a and (c+1)^b > a 
+ * Result found such that (c)**b <= a and (c+1)**b > a 
  *
- * This algorithm uses Newton's approximation x[i+1] = x[i] - f(x[i])/f'(x[i]) 
- * which will find the root in log(N) time where each step involves a fair bit.  This
- * is not meant to find huge roots [square and cube at most].
+ * This algorithm uses Newton's approximation 
+ * x[i+1] = x[i] - f(x[i])/f'(x[i]) 
+ * which will find the root in log(N) time where 
+ * each step involves a fair bit.  This is not meant to 
+ * find huge roots [square and cube, etc].
  */
 int
 mp_n_root (mp_int * a, mp_digit b, mp_int * c)
@@ -58,33 +60,39 @@ mp_n_root (mp_int * a, mp_digit b, mp_int * c)
       goto __T3;
     }
 
-    /* t2 = t1 - ((t1^b - a) / (b * t1^(b-1))) */
-    if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {	/* t3 = t1^(b-1) */
+    /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
+
+    /* t3 = t1**(b-1) */
+    if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {   
       goto __T3;
     }
 
     /* numerator */
-    if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) {	/* t2 = t1^b */
+    /* t2 = t1**b */
+    if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) {    
       goto __T3;
     }
 
-    if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) {	/* t2 = t1^b - a */
+    /* t2 = t1**b - a */
+    if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) {  
       goto __T3;
     }
 
-    if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) {	/* t3 = t1^(b-1) * b  */
+    /* denominator */
+    /* t3 = t1**(b-1) * b  */
+    if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) {    
       goto __T3;
     }
 
-    if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {	/* t3 = (t1^b - a)/(b * t1^(b-1)) */
+    /* t3 = (t1**b - a)/(b * t1**(b-1)) */
+    if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {  
       goto __T3;
     }
 
     if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
       goto __T3;
     }
-  }
-  while (mp_cmp (&t1, &t2) != MP_EQ);
+  }  while (mp_cmp (&t1, &t2) != MP_EQ);
 
   /* result can be off by a few so check */
   for (;;) {
@@ -94,7 +102,7 @@ mp_n_root (mp_int * a, mp_digit b, mp_int * c)
 
     if (mp_cmp (&t2, a) == MP_GT) {
       if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
-	goto __T3;
+    goto __T3;
       }
     } else {
       break;

bn_mp_reduce.c

@@ -32,8 +32,8 @@ mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
   /* q1 = x / b**(k-1)  */
   mp_rshd (&q, um - 1);         
 
-  /* according to HAC this is optimization is ok */
-  if (((unsigned long) m->used) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
+  /* according to HAC this optimization is ok */
+  if (((unsigned long) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
     if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) {
       goto CLEANUP;
     }
@@ -73,7 +73,7 @@ mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
   /* Back off if it's too big */
   while (mp_cmp (x, m) != MP_LT) {
     if ((res = s_mp_sub (x, m, x)) != MP_OKAY) {
-      break;
+      goto CLEANUP;
     }
   }
 

bn_radix.c

@@ -56,7 +56,9 @@ mp_read_radix (mp_int * a, char *str, int radix)
     }
     ++str;
   }
-  a->sign = neg;
+  if (mp_iszero(a) != 1) {
+     a->sign = neg;
+  }
   return MP_OKAY;
 }
 

bn_s_mp_sqr.c

@@ -39,7 +39,7 @@ s_mp_sqr (mp_int * a, mp_int * b)
     t.dp[2*ix] = (mp_digit) (r & ((mp_word) MP_MASK));
 
     /* get the carry */
-    u = (r >> ((mp_word) DIGIT_BIT));
+    u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
 
     /* left hand side of A[ix] * A[iy] */
     tmpx = a->dp[ix];
@@ -60,13 +60,13 @@ s_mp_sqr (mp_int * a, mp_int * b)
       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
 
       /* get carry */
-      u = (r >> ((mp_word) DIGIT_BIT));
+      u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
     }
     /* propagate upwards */
     while (u != ((mp_digit) 0)) {
       r = ((mp_word) * tmpt) + ((mp_word) u);
       *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
-      u = (r >> ((mp_word) DIGIT_BIT));
+      u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
     }
   }
 

changes.txt

@@ -1,3 +1,7 @@
+June 19th, 2003
+v0.21  -- Fixed bug in mp_mul_d which would not handle sign correctly [would not always forward it]
+       -- Removed the #line lines from gen.pl [was in violation of ISO C]       
+
 June 8th, 2003
 v0.20  -- Removed the book from the package.  Added the TDCAL license document.  
        -- This release is officially pure-bred TDCAL again [last officially TDCAL based release was v0.16]

demo/demo.c

@@ -162,6 +162,8 @@ int main(void)
          fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
       }
       fclose(log);
+
+      return 0;
 
       log = fopen("logs/sub.log", "w");
       for (cnt = 8; cnt <= 128; cnt += 8) {

etc/2kprime.1

@@ -1 +1,2 @@
-259-bits (k = 17745) = 926336713898529563388567880069503262826159877325124512315660672063305037101743
+256-bits (k = 36113) = 115792089237316195423570985008687907853269984665640564039457584007913129603823
+512-bits (k = 38117) = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006045979

etc/makefile.msvc

@@ -2,7 +2,7 @@
 #
 #Tom St Denis
 
-CFLAGS = /I../ /Ogityb2 /Gs /DWIN32 /W3
+CFLAGS = /I../ /Ox /DWIN32 /W3
 
 pprime: pprime.obj
 	cl pprime.obj ../tommath.lib 

gen.pl

@@ -9,7 +9,6 @@
 foreach my $filename (glob "bn*.c") {
    open( SRC, "<$filename" ) or die "Couldn't open $filename for reading: $!";
    print OUT "/* Start: $filename */\n";
-   print OUT qq[#line 0 "$filename"\n];
    print OUT while <SRC>;
    print OUT "\n/* End: $filename */\n\n";
    close SRC or die "Error closing $filename after reading: $!";

logs/add.log

@@ -1,16 +0,0 @@
-224  11069160
-448   9156136
-672   8089755
-896   7399424
-1120   6389352
-1344   5818648
-1568   5257112
-1792   4982160
-2016   4527856
-2240   4325312
-2464   4051760
-2688   3767640
-2912   3612520
-3136   3415208
-3360   3258656
-3584   3113360