import libtommath 1.0

2 files changed, 552 insertions, 286 deletions

diff --git a/libtommath/demo/demo.c b/libtommath/demo/demo.c
index 3e5663b..b46b7f8 100644
--- a/libtommath/demo/demo.c
+++ b/libtommath/demo/demo.c
@@ -1,3 +1,4 @@
+#include <string.h>
 #include <time.h>
 
 #ifdef IOWNANATHLON
@@ -7,6 +8,28 @@
 #define SLEEP
 #endif
 
+/*
+ * Configuration
+ */
+#ifndef LTM_DEMO_TEST_VS_MTEST
+#define LTM_DEMO_TEST_VS_MTEST 1
+#endif
+
+#ifndef LTM_DEMO_TEST_REDUCE_2K_L
+/* This test takes a moment so we disable it by default, but it can be:
+ * 0 to disable testing
+ * 1 to make the test with P = 2^1024 - 0x2A434 B9FDEC95 D8F9D550 FFFFFFFF FFFFFFFF
+ * 2 to make the test with P = 2^2048 - 0x1 00000000 00000000 00000000 00000000 4945DDBF 8EA2A91D 5776399B B83E188F
+ */
+#define LTM_DEMO_TEST_REDUCE_2K_L 0
+#endif
+
+#ifdef LTM_DEMO_REAL_RAND
+#define LTM_DEMO_RAND_SEED  time(NULL)
+#else
+#define LTM_DEMO_RAND_SEED  23
+#endif
+
 #include "tommath.h"
 
 void ndraw(mp_int * a, char *name)
@@ -16,12 +39,16 @@ void ndraw(mp_int * a, char *name)
    printf("%s: ", name);
    mp_toradix(a, buf, 10);
    printf("%s\n", buf);
+   mp_toradix(a, buf, 16);
+   printf("0x%s\n", buf);
 }
 
+#if LTM_DEMO_TEST_VS_MTEST
 static void draw(mp_int * a)
 {
    ndraw(a, "");
 }
+#endif
 
 
 unsigned long lfsr = 0xAAAAAAAAUL;
@@ -37,183 +64,392 @@ int lbit(void)
    }
 }
 
+#if defined(LTM_DEMO_REAL_RAND) && !defined(_WIN32)
+static FILE* fd_urandom;
+#endif
 int myrng(unsigned char *dst, int len, void *dat)
 {
    int x;
-
-   for (x = 0; x < len; x++)
-      dst[x] = rand() & 0xFF;
+   (void)dat;
+#if defined(LTM_DEMO_REAL_RAND)
+   if (!fd_urandom) {
+#if !defined(_WIN32)
+      fprintf(stderr, "\nno /dev/urandom\n");
+#endif
+   }
+   else {
+      return fread(dst, 1, len, fd_urandom);
+   }
+#endif
+   for (x = 0; x < len; ) {
+      unsigned int r = (unsigned int)rand();
+      do {
+         dst[x++] = r & 0xFF;
+         r >>= 8;
+      } while((r != 0) && (x < len));
+   }
    return len;
 }
 
+#if LTM_DEMO_TEST_VS_MTEST != 0
+static void _panic(int l)
+{
+  fprintf(stderr, "\n%d: fgets failed\n", l);
+  exit(EXIT_FAILURE);
+}
+#endif
+
+mp_int a, b, c, d, e, f;
+
+static void _cleanup(void)
+{
+  mp_clear_multi(&a, &b, &c, &d, &e, &f, NULL);
+  printf("\n");
 
+#ifdef LTM_DEMO_REAL_RAND
+  if(fd_urandom)
+     fclose(fd_urandom);
+#endif
+}
+struct mp_sqrtmod_prime_st {
+   unsigned long p;
+   unsigned long n;
+   mp_digit r;
+};
+struct mp_sqrtmod_prime_st sqrtmod_prime[] = {
+      { 5, 14, 3 },
+      { 7, 9, 4 },
+      { 113, 2, 62 }
+};
+struct mp_jacobi_st {
+   unsigned long n;
+   int c[16];
+};
+struct mp_jacobi_st jacobi[] = {
+      { 3, {  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1 } },
+      { 5, {  0,  1, -1, -1,  1,  0,  1, -1, -1,  1,  0,  1, -1, -1,  1,  0 } },
+      { 7, {  1, -1,  1, -1, -1,  0,  1,  1, -1,  1, -1, -1,  0,  1,  1, -1 } },
+      { 9, { -1,  1,  0,  1,  1,  0,  1,  1,  0,  1,  1,  0,  1,  1,  0,  1 } },
+};
 
 char cmd[4096], buf[4096];
 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, 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;
+   int cnt, ix;
+#if LTM_DEMO_TEST_VS_MTEST
+   unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n,
+      gcd_n, lcm_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n;
+   char* ret;
+#else
+   unsigned long s, t;
+   unsigned long long q, r;
    mp_digit mp;
+   int i, n, err, should;
+#endif
 
+   if (mp_init_multi(&a, &b, &c, &d, &e, &f, NULL)!= MP_OKAY)
+     return EXIT_FAILURE;
 
-   mp_init(&a);
-   mp_init(&b);
-   mp_init(&c);
-   mp_init(&d);
-   mp_init(&e);
-   mp_init(&f);
-
-   srand(time(NULL));
-
-#if 0
-   // test montgomery 
-   printf("Testing montgomery...\n");
-   for (i = 1; i < 10; i++) {
-      printf("Testing digit size: %d\n", i);
-      for (n = 0; n < 1000; n++) {
-         mp_rand(&a, i);
-         a.dp[0] |= 1;
+   atexit(_cleanup);
 
-         // let's see if R is right
-         mp_montgomery_calc_normalization(&b, &a);
-         mp_montgomery_setup(&a, &mp);
+#if defined(LTM_DEMO_REAL_RAND)
+   if (!fd_urandom) {
+      fd_urandom = fopen("/dev/urandom", "r");
+      if (!fd_urandom) {
+#if !defined(_WIN32)
+         fprintf(stderr, "\ncould not open /dev/urandom\n");
+#endif
+      }
+   }
+#endif
+   srand(LTM_DEMO_RAND_SEED);
 
-         // now test a random reduction 
-         for (ix = 0; ix < 100; ix++) {
-             mp_rand(&c, 1 + abs(rand()) % (2*i));
-             mp_copy(&c, &d);
-             mp_copy(&c, &e);
+#ifdef MP_8BIT
+   printf("Digit size 8 Bit \n");
+#endif
+#ifdef MP_16BIT
+   printf("Digit size 16 Bit \n");
+#endif
+#ifdef MP_32BIT
+   printf("Digit size 32 Bit \n");
+#endif
+#ifdef MP_64BIT
+   printf("Digit size 64 Bit \n");
+#endif
+   printf("Size of mp_digit: %u\n", (unsigned int)sizeof(mp_digit));
+   printf("Size of mp_word: %u\n", (unsigned int)sizeof(mp_word));
+   printf("DIGIT_BIT: %d\n", DIGIT_BIT);
+   printf("MP_PREC: %d\n", MP_PREC);
+
+#if LTM_DEMO_TEST_VS_MTEST == 0
+   // trivial stuff
+   mp_set_int(&a, 5);
+   mp_neg(&a, &b);
+   if (mp_cmp(&a, &b) != MP_GT) {
+      return EXIT_FAILURE;
+   }
+   if (mp_cmp(&b, &a) != MP_LT) {
+      return EXIT_FAILURE;
+   }
+   mp_neg(&a, &a);
+   if (mp_cmp(&b, &a) != MP_EQ) {
+      return EXIT_FAILURE;
+   }
+   mp_abs(&a, &b);
+   if (mp_isneg(&b) != MP_NO) {
+      return EXIT_FAILURE;
+   }
+   mp_add_d(&a, 1, &b);
+   mp_add_d(&a, 6, &b);
 
-             mp_mod(&d, &a, &d);
-             mp_montgomery_reduce(&c, &a, mp);
-             mp_mulmod(&c, &b, &a, &c);
 
-             if (mp_cmp(&c, &d) != MP_EQ) { 
-printf("d = e mod a, c = e MOD a\n");
-mp_todecimal(&a, buf); printf("a = %s\n", buf);
-mp_todecimal(&e, buf); printf("e = %s\n", buf);
-mp_todecimal(&d, buf); printf("d = %s\n", buf);
-mp_todecimal(&c, buf); printf("c = %s\n", buf);
-printf("compare no compare!\n"); exit(EXIT_FAILURE); }
+   mp_set_int(&a, 0);
+   mp_set_int(&b, 1);
+   if ((err = mp_jacobi(&a, &b, &i)) != MP_OKAY) {
+      printf("Failed executing mp_jacobi(0 | 1) %s.\n", mp_error_to_string(err));
+      return EXIT_FAILURE;
+   }
+   if (i != 1) {
+      printf("Failed trivial mp_jacobi(0 | 1) %d != 1\n", i);
+      return EXIT_FAILURE;
+   }
+   for (cnt = 0; cnt < (int)(sizeof(jacobi)/sizeof(jacobi[0])); ++cnt) {
+      mp_set_int(&b, jacobi[cnt].n);
+      /* only test positive values of a */
+      for (n = -5; n <= 10; ++n) {
+         mp_set_int(&a, abs(n));
+         should = MP_OKAY;
+         if (n < 0) {
+            mp_neg(&a, &a);
+            /* Until #44 is fixed the negative a's must fail */
+            should = MP_VAL;
+         }
+         if ((err = mp_jacobi(&a, &b, &i)) != should) {
+            printf("Failed executing mp_jacobi(%d | %lu) %s.\n", n, jacobi[cnt].n, mp_error_to_string(err));
+            return EXIT_FAILURE;
+         }
+         if (err == MP_OKAY && i != jacobi[cnt].c[n + 5]) {
+            printf("Failed trivial mp_jacobi(%d | %lu) %d != %d\n", n, jacobi[cnt].n, i, jacobi[cnt].c[n + 5]);
+            return EXIT_FAILURE;
          }
       }
    }
-   printf("done\n");
 
    // test mp_get_int
-   printf("Testing: mp_get_int\n");
+   printf("\n\nTesting: mp_get_int");
    for (i = 0; i < 1000; ++i) {
-      t = ((unsigned long) rand() * rand() + 1) & 0xFFFFFFFF;
-      mp_set_int(&a, t);
-      if (t != mp_get_int(&a)) {
-	 printf("mp_get_int() bad result!\n");
-	 return 1;
+      t = ((unsigned long) rand () * rand () + 1) & 0xFFFFFFFF;
+      mp_set_int (&a, t);
+      if (t != mp_get_int (&a)) {
+         printf ("\nmp_get_int() bad result!");
+         return EXIT_FAILURE;
       }
    }
    mp_set_int(&a, 0);
    if (mp_get_int(&a) != 0) {
-      printf("mp_get_int() bad result!\n");
-      return 1;
+      printf("\nmp_get_int() bad result!");
+      return EXIT_FAILURE;
    }
    mp_set_int(&a, 0xffffffff);
    if (mp_get_int(&a) != 0xffffffff) {
-      printf("mp_get_int() bad result!\n");
-      return 1;
+      printf("\nmp_get_int() bad result!");
+      return EXIT_FAILURE;
    }
+
+   printf("\n\nTesting: mp_get_long\n");
+   for (i = 0; i < (int)(sizeof(unsigned long)*CHAR_BIT) - 1; ++i) {
+      t = (1ULL << (i+1)) - 1;
+      if (!t)
+         t = -1;
+      printf(" t = 0x%lx i = %d\r", t, i);
+      do {
+         if (mp_set_long(&a, t) != MP_OKAY) {
+            printf("\nmp_set_long() error!");
+            return EXIT_FAILURE;
+         }
+         s = mp_get_long(&a);
+         if (s != t) {
+            printf("\nmp_get_long() bad result! 0x%lx != 0x%lx", s, t);
+            return EXIT_FAILURE;
+         }
+         t <<= 1;
+      } while(t);
+   }
+
+   printf("\n\nTesting: mp_get_long_long\n");
+   for (i = 0; i < (int)(sizeof(unsigned long long)*CHAR_BIT) - 1; ++i) {
+      r = (1ULL << (i+1)) - 1;
+      if (!r)
+         r = -1;
+      printf(" r = 0x%llx i = %d\r", r, i);
+      do {
+         if (mp_set_long_long(&a, r) != MP_OKAY) {
+            printf("\nmp_set_long_long() error!");
+            return EXIT_FAILURE;
+         }
+         q = mp_get_long_long(&a);
+         if (q != r) {
+            printf("\nmp_get_long_long() bad result! 0x%llx != 0x%llx", q, r);
+            return EXIT_FAILURE;
+         }
+         r <<= 1;
+      } while(r);
+   }
+
    // test mp_sqrt
-   printf("Testing: mp_sqrt\n");
+   printf("\n\nTesting: mp_sqrt\n");
    for (i = 0; i < 1000; ++i) {
-      printf("%6d\r", i);
-      fflush(stdout);
-      n = (rand() & 15) + 1;
-      mp_rand(&a, n);
-      if (mp_sqrt(&a, &b) != MP_OKAY) {
-	 printf("mp_sqrt() error!\n");
-	 return 1;
+      printf ("%6d\r", i);
+      fflush (stdout);
+      n = (rand () & 15) + 1;
+      mp_rand (&a, n);
+      if (mp_sqrt (&a, &b) != MP_OKAY) {
+         printf ("\nmp_sqrt() error!");
+         return EXIT_FAILURE;
+      }
+      mp_n_root_ex (&a, 2, &c, 0);
+      mp_n_root_ex (&a, 2, &d, 1);
+      if (mp_cmp_mag (&c, &d) != MP_EQ) {
+         printf ("\nmp_n_root_ex() bad result!");
+         return EXIT_FAILURE;
       }
-      mp_n_root(&a, 2, &a);
-      if (mp_cmp_mag(&b, &a) != MP_EQ) {
-	 printf("mp_sqrt() bad result!\n");
-	 return 1;
+      if (mp_cmp_mag (&b, &c) != MP_EQ) {
+         printf ("mp_sqrt() bad result!\n");
+         return EXIT_FAILURE;
       }
    }
 
-   printf("\nTesting: mp_is_square\n");
+   printf("\n\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;
-      mp_rand(&a, n);
-      mp_sqr(&a, &a);
-      if (mp_is_square(&a, &n) != MP_OKAY) {
-	 printf("fn:mp_is_square() error!\n");
-	 return 1;
+      n = (rand () & 7) + 1;
+      mp_rand (&a, n);
+      mp_sqr (&a, &a);
+      if (mp_is_square (&a, &n) != MP_OKAY) {
+         printf ("\nfn:mp_is_square() error!");
+         return EXIT_FAILURE;
       }
       if (n == 0) {
-	 printf("fn:mp_is_square() bad result!\n");
-	 return 1;
+         printf ("\nfn:mp_is_square() bad result!");
+         return EXIT_FAILURE;
       }
 
       /* test for false positives */
-      mp_add_d(&a, 1, &a);
-      if (mp_is_square(&a, &n) != MP_OKAY) {
-	 printf("fp:mp_is_square() error!\n");
-	 return 1;
+      mp_add_d (&a, 1, &a);
+      if (mp_is_square (&a, &n) != MP_OKAY) {
+         printf ("\nfp:mp_is_square() error!");
+         return EXIT_FAILURE;
       }
       if (n == 1) {
-	 printf("fp:mp_is_square() bad result!\n");
-	 return 1;
+         printf ("\nfp:mp_is_square() bad result!");
+         return EXIT_FAILURE;
       }
 
    }
    printf("\n\n");
 
+   // r^2 = n (mod p)
+   for (i = 0; i < (int)(sizeof(sqrtmod_prime)/sizeof(sqrtmod_prime[0])); ++i) {
+      mp_set_int(&a, sqrtmod_prime[i].p);
+      mp_set_int(&b, sqrtmod_prime[i].n);
+      if (mp_sqrtmod_prime(&b, &a, &c) != MP_OKAY) {
+         printf("Failed executing %d. mp_sqrtmod_prime\n", (i+1));
+         return EXIT_FAILURE;
+      }
+      if (mp_cmp_d(&c, sqrtmod_prime[i].r) != MP_EQ) {
+         printf("Failed %d. trivial mp_sqrtmod_prime\n", (i+1));
+         ndraw(&c, "r");
+         return EXIT_FAILURE;
+      }
+   }
+
    /* test for size */
    for (ix = 10; ix < 128; ix++) {
-      printf("Testing (not safe-prime): %9d bits    \r", ix);
-      fflush(stdout);
-      err =
-	 mp_prime_random_ex(&a, 8, ix,
-			    (rand() & 1) ? LTM_PRIME_2MSB_OFF :
-			    LTM_PRIME_2MSB_ON, myrng, NULL);
+      printf ("Testing (not safe-prime): %9d bits    \r", ix);
+      fflush (stdout);
+      err = mp_prime_random_ex (&a, 8, ix,
+                                (rand () & 1) ? 0 : LTM_PRIME_2MSB_ON, myrng,
+                                NULL);
       if (err != MP_OKAY) {
-	 printf("failed with err code %d\n", err);
-	 return EXIT_FAILURE;
+         printf ("failed with err code %d\n", err);
+         return EXIT_FAILURE;
       }
-      if (mp_count_bits(&a) != ix) {
-	 printf("Prime is %d not %d bits!!!\n", mp_count_bits(&a), ix);
-	 return EXIT_FAILURE;
+      if (mp_count_bits (&a) != ix) {
+         printf ("Prime is %d not %d bits!!!\n", mp_count_bits (&a), ix);
+         return EXIT_FAILURE;
       }
    }
+   printf("\n");
 
    for (ix = 16; ix < 128; ix++) {
-      printf("Testing (   safe-prime): %9d bits    \r", ix);
-      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);
+      printf ("Testing (    safe-prime): %9d bits    \r", ix);
+      fflush (stdout);
+      err = mp_prime_random_ex (
+            &a, 8, ix, ((rand () & 1) ? 0 : LTM_PRIME_2MSB_ON) | LTM_PRIME_SAFE,
+            myrng, NULL);
       if (err != MP_OKAY) {
-	 printf("failed with err code %d\n", err);
-	 return EXIT_FAILURE;
+         printf ("failed with err code %d\n", err);
+         return EXIT_FAILURE;
       }
-      if (mp_count_bits(&a) != ix) {
-	 printf("Prime is %d not %d bits!!!\n", mp_count_bits(&a), ix);
-	 return EXIT_FAILURE;
+      if (mp_count_bits (&a) != ix) {
+         printf ("Prime is %d not %d bits!!!\n", mp_count_bits (&a), ix);
+         return EXIT_FAILURE;
       }
       /* let's see if it's really a safe prime */
-      mp_sub_d(&a, 1, &a);
-      mp_div_2(&a, &a);
-      mp_prime_is_prime(&a, 8, &cnt);
+      mp_sub_d (&a, 1, &a);
+      mp_div_2 (&a, &a);
+      mp_prime_is_prime (&a, 8, &cnt);
       if (cnt != MP_YES) {
-	 printf("sub is not prime!\n");
-	 return EXIT_FAILURE;
+         printf ("sub is not prime!\n");
+         return EXIT_FAILURE;
+      }
+   }
+
+   printf("\n\n");
+
+   // test montgomery
+   printf("Testing: montgomery...\n");
+   for (i = 1; i <= 10; i++) {
+      if (i == 10)
+         i = 1000;
+      printf(" digit size: %2d\r", i);
+      fflush(stdout);
+      for (n = 0; n < 1000; n++) {
+         mp_rand(&a, i);
+         a.dp[0] |= 1;
+
+         // let's see if R is right
+         mp_montgomery_calc_normalization(&b, &a);
+         mp_montgomery_setup(&a, &mp);
+
+         // now test a random reduction
+         for (ix = 0; ix < 100; ix++) {
+             mp_rand(&c, 1 + abs(rand()) % (2*i));
+             mp_copy(&c, &d);
+             mp_copy(&c, &e);
+
+             mp_mod(&d, &a, &d);
+             mp_montgomery_reduce(&c, &a, mp);
+             mp_mulmod(&c, &b, &a, &c);
+
+             if (mp_cmp(&c, &d) != MP_EQ) {
+printf("d = e mod a, c = e MOD a\n");
+mp_todecimal(&a, buf); printf("a = %s\n", buf);
+mp_todecimal(&e, buf); printf("e = %s\n", buf);
+mp_todecimal(&d, buf); printf("d = %s\n", buf);
+mp_todecimal(&c, buf); printf("c = %s\n", buf);
+printf("compare no compare!\n"); return EXIT_FAILURE; }
+             /* only one big montgomery reduction */
+             if (i > 10)
+             {
+                n = 1000;
+                ix = 100;
+             }
+         }
       }
    }
 
@@ -239,120 +475,123 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 #endif
 
    /* test mp_cnt_lsb */
-   printf("testing mp_cnt_lsb...\n");
+   printf("\n\nTesting: mp_cnt_lsb");
    mp_set(&a, 1);
    for (ix = 0; ix < 1024; ix++) {
-      if (mp_cnt_lsb(&a) != ix) {
-	 printf("Failed at %d, %d\n", ix, mp_cnt_lsb(&a));
-	 return 0;
+      if (mp_cnt_lsb (&a) != ix) {
+         printf ("Failed at %d, %d\n", ix, mp_cnt_lsb (&a));
+         return EXIT_FAILURE;
       }
-      mp_mul_2(&a, &a);
+      mp_mul_2 (&a, &a);
    }
 
 /* test mp_reduce_2k */
-   printf("Testing mp_reduce_2k...\n");
+   printf("\n\nTesting: 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 */
-
+      mp_2expt (&a, cnt);
+      mp_sub_d (&a, 2, &a); /* a = 2**cnt - 2 */
 
-      printf("\nTesting %4d bits", cnt);
-      printf("(%d)", mp_reduce_is_2k(&a));
-      mp_reduce_2k_setup(&a, &tmp);
-      printf("(%d)", tmp);
+      printf ("\r %4d bits", cnt);
+      printf ("(%d)", mp_reduce_is_2k (&a));
+      mp_reduce_2k_setup (&a, &tmp);
+      printf ("(%lu)", (unsigned long) tmp);
       for (ix = 0; ix < 1000; ix++) {
-	 if (!(ix & 127)) {
-	    printf(".");
-	    fflush(stdout);
-	 }
-	 mp_rand(&b, (cnt / DIGIT_BIT + 1) * 2);
-	 mp_copy(&c, &b);
-	 mp_mod(&c, &a, &c);
-	 mp_reduce_2k(&b, &a, 2);
-	 if (mp_cmp(&c, &b)) {
-	    printf("FAILED\n");
-	    exit(0);
-	 }
+         if (!(ix & 127)) {
+            printf (".");
+            fflush (stdout);
+         }
+         mp_rand (&b, (cnt / DIGIT_BIT + 1) * 2);
+         mp_copy (&c, &b);
+         mp_mod (&c, &a, &c);
+         mp_reduce_2k (&b, &a, 2);
+         if (mp_cmp (&c, &b)) {
+            printf ("FAILED\n");
+            return EXIT_FAILURE;
+         }
       }
    }
 
 /* test mp_div_3  */
-   printf("Testing mp_div_3...\n");
+   printf("\n\nTesting: mp_div_3...\n");
    mp_set(&d, 3);
    for (cnt = 0; cnt < 10000;) {
-      mp_digit r1, r2;
+      mp_digit r2;
 
       if (!(++cnt & 127))
-	 printf("%9d\r", cnt);
+      {
+        printf("%9d\r", cnt);
+        fflush(stdout);
+      }
       mp_rand(&a, abs(rand()) % 128 + 1);
       mp_div(&a, &d, &b, &e);
       mp_div_3(&a, &c, &r2);
 
       if (mp_cmp(&b, &c) || mp_cmp_d(&e, r2)) {
-	 printf("\n\nmp_div_3 => Failure\n");
+	 printf("\nmp_div_3 => Failure\n");
       }
    }
-   printf("\n\nPassed div_3 testing\n");
+   printf("\nPassed div_3 testing");
 
 /* test the DR reduction */
-   printf("testing mp_dr_reduce...\n");
+   printf("\n\nTesting: mp_dr_reduce...\n");
    for (cnt = 2; cnt < 32; cnt++) {
-      printf("%d digit modulus\n", cnt);
-      mp_grow(&a, cnt);
-      mp_zero(&a);
+      printf ("\r%d digit modulus", cnt);
+      mp_grow (&a, cnt);
+      mp_zero (&a);
       for (ix = 1; ix < cnt; ix++) {
-	 a.dp[ix] = MP_MASK;
+         a.dp[ix] = MP_MASK;
       }
       a.used = cnt;
       a.dp[0] = 3;
 
-      mp_rand(&b, cnt - 1);
-      mp_copy(&b, &c);
+      mp_rand (&b, cnt - 1);
+      mp_copy (&b, &c);
 
       rr = 0;
       do {
-	 if (!(rr & 127)) {
-	    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 (!(rr & 127)) {
+            printf (".");
+            fflush (stdout);
+         }
+         mp_sqr (&b, &b);
+         mp_add_d (&b, 1, &b);
+         mp_copy (&b, &c);
 
-	 if (mp_cmp(&b, &c) != MP_EQ) {
-	    printf("Failed on trial %lu\n", rr);
-	    exit(-1);
+         mp_mod (&b, &a, &b);
+         mp_dr_setup(&a, &mp),
+         mp_dr_reduce (&c, &a, mp);
 
-	 }
+         if (mp_cmp (&b, &c) != MP_EQ) {
+            printf ("Failed on trial %u\n", rr);
+            return EXIT_FAILURE;
+         }
       } while (++rr < 500);
-      printf("Passed DR test for %d digits\n", cnt);
+      printf (" passed");
+      fflush (stdout);
    }
 
-#endif
-
+#if LTM_DEMO_TEST_REDUCE_2K_L
 /* test the mp_reduce_2k_l code */
-#if 0
-#if 0
+#if LTM_DEMO_TEST_REDUCE_2K_L == 1
 /* 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
+#elif LTM_DEMO_TEST_REDUCE_2K_L == 2
 /*  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);
+#else
+#error oops
 #endif
 
    mp_todecimal(&a, buf);
-   printf("p==%s\n", buf);
+   printf("\n\np==%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");
@@ -363,9 +602,9 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
    mp_rand(&b, 64);
    mp_mod(&b, &a, &b);
    mp_copy(&b, &c);
-   printf("testing mp_reduce_2k_l...");
+   printf("Testing: mp_reduce_2k_l...");
    fflush(stdout);
-   for (cnt = 0; cnt < (1UL << 20); cnt++) {
+   for (cnt = 0; cnt < (int)(1UL << 20); cnt++) {
       mp_sqr(&b, &b);
       mp_add_d(&b, 1, &b);
       mp_reduce_2k_l(&b, &a, &d);
@@ -373,7 +612,7 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
       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);
+	 printf("mp_reduce_2k_l() failed at step %d\n", cnt);
 	 mp_tohex(&b, buf);
 	 printf("b == %s\n", buf);
 	 mp_tohex(&c, buf);
@@ -382,7 +621,9 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
       }
    }
    printf("...Passed\n");
-#endif
+#endif /* LTM_DEMO_TEST_REDUCE_2K_L */
+
+#else
 
    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 =
@@ -428,17 +669,17 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	 ("%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);
+      ret=fgets(cmd, 4095, stdin); if(!ret){_panic(__LINE__);}
       cmd[strlen(cmd) - 1] = 0;
-      printf("%s  ]\r", cmd);
+      printf("%-6s ]\r", cmd);
       fflush(stdout);
       if (!strcmp(cmd, "mul2d")) {
 	 ++mul2d_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 sscanf(buf, "%d", &rr);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 
 	 mp_mul_2d(&a, rr, &a);
@@ -447,15 +688,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    printf("mul2d failed, rr == %d\n", rr);
 	    draw(&a);
 	    draw(&b);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "div2d")) {
 	 ++div2d_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 sscanf(buf, "%d", &rr);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 
 	 mp_div_2d(&a, rr, &a, &e);
@@ -467,15 +708,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    printf("div2d failed, rr == %d\n", rr);
 	    draw(&a);
 	    draw(&b);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "add")) {
 	 ++add_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_add(&d, &b, &d);
@@ -485,7 +726,7 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
 
 	 /* test the sign/unsigned storage functions */
@@ -498,7 +739,7 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    printf("mp_signed_bin failure!\n");
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
 
 
@@ -510,16 +751,16 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    printf("mp_unsigned_bin failure!\n");
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
 
       } else if (!strcmp(cmd, "sub")) {
 	 ++sub_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_sub(&d, &b, &d);
@@ -529,15 +770,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "mul")) {
 	 ++mul_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_mul(&d, &b, &d);
@@ -547,17 +788,17 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "div")) {
 	 ++div_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&d, buf, 64);
 
 	 mp_div(&a, &b, &e, &f);
@@ -570,14 +811,14 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&d);
 	    draw(&e);
 	    draw(&f);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
 
       } else if (!strcmp(cmd, "sqr")) {
 	 ++sqr_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 	 mp_copy(&a, &c);
 	 mp_sqr(&c, &c);
@@ -586,15 +827,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&a);
 	    draw(&b);
 	    draw(&c);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "gcd")) {
 	 ++gcd_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_gcd(&d, &b, &d);
@@ -605,15 +846,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "lcm")) {
 	 ++lcm_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
 	 mp_copy(&a, &d);
 	 mp_lcm(&d, &b, &d);
@@ -624,17 +865,17 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    draw(&d);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "expt")) {
 	 ++expt_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&d, buf, 64);
 	 mp_copy(&a, &e);
 	 mp_exptmod(&e, &b, &c, &e);
@@ -645,15 +886,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&c);
 	    draw(&d);
 	    draw(&e);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "invmod")) {
 	 ++inv_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&c, buf, 64);
 	 mp_invmod(&a, &b, &d);
 	 mp_mulmod(&d, &a, &b, &e);
@@ -663,16 +904,17 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    draw(&d);
+	    draw(&e);
 	    mp_gcd(&a, &b, &e);
 	    draw(&e);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
 
       } else if (!strcmp(cmd, "div2")) {
 	 ++div2_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 	 mp_div_2(&a, &c);
 	 if (mp_cmp(&c, &b) != MP_EQ) {
@@ -680,13 +922,13 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&a);
 	    draw(&b);
 	    draw(&c);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "mul2")) {
 	 ++mul2_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 	 mp_mul_2(&a, &c);
 	 if (mp_cmp(&c, &b) != MP_EQ) {
@@ -694,15 +936,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&a);
 	    draw(&b);
 	    draw(&c);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "add_d")) {
 	 ++add_d_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 sscanf(buf, "%d", &ix);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 	 mp_add_d(&a, ix, &c);
 	 if (mp_cmp(&b, &c) != MP_EQ) {
@@ -711,15 +953,15 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    printf("d == %d\n", ix);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
       } else if (!strcmp(cmd, "sub_d")) {
 	 ++sub_d_n;
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&a, buf, 64);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 sscanf(buf, "%d", &ix);
-	 fgets(buf, 4095, stdin);
+	 ret=fgets(buf, 4095, stdin); if(!ret){_panic(__LINE__);}
 	 mp_read_radix(&b, buf, 64);
 	 mp_sub_d(&a, ix, &c);
 	 if (mp_cmp(&b, &c) != MP_EQ) {
@@ -728,13 +970,17 @@ printf("compare no compare!\n"); exit(EXIT_FAILURE); }
 	    draw(&b);
 	    draw(&c);
 	    printf("d == %d\n", ix);
-	    return 0;
+	    return EXIT_FAILURE;
 	 }
+      } else if (!strcmp(cmd, "exit")) {
+         printf("\nokay, exiting now\n");
+         break;
       }
    }
+#endif
    return 0;
 }
 
 /* $Source$ */
-/* $Revision: 0.36 $ */
-/* $Date: 2005-08-01 16:37:28 +0000 $ */
+/* $Revision$ */
+/* $Date$ */
diff --git a/libtommath/demo/timing.c b/libtommath/demo/timing.c
index 57bb6d4..1bd8489 100644
--- a/libtommath/demo/timing.c
+++ b/libtommath/demo/timing.c
@@ -1,5 +1,6 @@
 #include <tommath.h>
 #include <time.h>
+#include <unistd.h>
 
 ulong64 _tt;
 
@@ -10,6 +11,12 @@ ulong64 _tt;
 #define SLEEP
 #endif
 
+#ifdef LTM_TIMING_REAL_RAND
+#define LTM_TIMING_RAND_SEED  time(NULL)
+#else
+#define LTM_TIMING_RAND_SEED  23
+#endif
+
 
 void ndraw(mp_int * a, char *name)
 {
@@ -44,10 +51,12 @@ 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;
+  /* version from http://www.mcs.anl.gov/~kazutomo/rdtsc.html
+   * the old code always got a warning issued by gcc, clang did not complain...
+   */
+  unsigned hi, lo;
+  __asm__ __volatile__ ("rdtsc" : "=a"(lo), "=d"(hi));
+  return ((ulong64)lo)|( ((ulong64)hi)<<32);
 #else /* gcc-IA64 version */
    unsigned long result;
    __asm__ __volatile__("mov %0=ar.itc":"=r"(result)::"memory");
@@ -78,12 +87,24 @@ static ulong64 TIMFUNC(void)
 //#define DO8(x) DO4(x); DO4(x);
 //#define DO(x)  DO8(x); DO8(x);
 
+#ifdef TIMING_NO_LOGS
+#define FOPEN(a, b)     NULL
+#define FPRINTF(a,b,c,d)
+#define FFLUSH(a)
+#define FCLOSE(a)       (void)(a)
+#else
+#define FOPEN(a,b)       fopen(a,b)
+#define FPRINTF(a,b,c,d) fprintf(a,b,c,d)
+#define FFLUSH(a)        fflush(a)
+#define FCLOSE(a)        fclose(a)
+#endif
+
 int main(void)
 {
    ulong64 tt, gg, CLK_PER_SEC;
    FILE *log, *logb, *logc, *logd;
    mp_int a, b, c, d, e, f;
-   int n, cnt, ix, old_kara_m, old_kara_s;
+   int n, cnt, ix, old_kara_m, old_kara_s, old_toom_m, old_toom_s;
    unsigned rr;
 
    mp_init(&a);
@@ -93,19 +114,15 @@ int main(void)
    mp_init(&e);
    mp_init(&f);
 
-   srand(time(NULL));
-
+   srand(LTM_TIMING_RAND_SEED);
 
-   /* temp. turn off TOOM */
-   TOOM_MUL_CUTOFF = TOOM_SQR_CUTOFF = 100000;
 
    CLK_PER_SEC = TIMFUNC();
    sleep(1);
    CLK_PER_SEC = TIMFUNC() - CLK_PER_SEC;
 
    printf("CLK_PER_SEC == %llu\n", CLK_PER_SEC);
-   goto exptmod;
-   log = fopen("logs/add.log", "w");
+   log = FOPEN("logs/add.log", "w");
    for (cnt = 8; cnt <= 128; cnt += 8) {
       SLEEP;
       mp_rand(&a, cnt);
@@ -121,12 +138,12 @@ int main(void)
       } while (++rr < 100000);
       printf("Adding\t\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);
-      fflush(log);
+      FPRINTF(log, "%d %9llu\n", cnt * DIGIT_BIT, tt);
+      FFLUSH(log);
    }
-   fclose(log);
+   FCLOSE(log);
 
-   log = fopen("logs/sub.log", "w");
+   log = FOPEN("logs/sub.log", "w");
    for (cnt = 8; cnt <= 128; cnt += 8) {
       SLEEP;
       mp_rand(&a, cnt);
@@ -143,22 +160,26 @@ int main(void)
 
       printf("Subtracting\t\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);
-      fflush(log);
+      FPRINTF(log, "%d %9llu\n", cnt * DIGIT_BIT, tt);
+      FFLUSH(log);
    }
-   fclose(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 < 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");
+   /* currently toom-cook cut-off is too high to kick in, so we just use the karatsuba values */
+   old_toom_m = old_kara_m;
+   old_toom_s = old_kara_m;
+   for (ix = 0; ix < 3; ix++) {
+      printf("With%s Karatsuba, With%s Toom\n", (ix == 0) ? "out" : "", (ix == 1) ? "out" : "");
+
+      KARATSUBA_MUL_CUTOFF = (ix == 1) ? old_kara_m : 9999;
+      KARATSUBA_SQR_CUTOFF = (ix == 1) ? old_kara_s : 9999;
+      TOOM_MUL_CUTOFF = (ix == 2) ? old_toom_m : 9999;
+      TOOM_SQR_CUTOFF = (ix == 2) ? old_toom_s : 9999;
+
+      log = FOPEN((ix == 0) ? "logs/mult.log" : (ix == 1) ? "logs/mult_kara.log" : "logs/mult_toom.log", "w");
       for (cnt = 4; cnt <= 10240 / DIGIT_BIT; cnt += 2) {
 	 SLEEP;
 	 mp_rand(&a, cnt);
@@ -174,12 +195,12 @@ int main(void)
 	 } while (++rr < 100);
 	 printf("Multiplying\t%4d-bit => %9llu/sec, %9llu cycles\n",
 		mp_count_bits(&a), CLK_PER_SEC / tt, tt);
-	 fprintf(log, "%d %9llu\n", mp_count_bits(&a), tt);
-	 fflush(log);
+	 FPRINTF(log, "%d %9llu\n", mp_count_bits(&a), tt);
+	 FFLUSH(log);
       }
-      fclose(log);
+      FCLOSE(log);
 
-      log = fopen((ix == 0) ? "logs/sqr.log" : "logs/sqr_kara.log", "w");
+      log = FOPEN((ix == 0) ? "logs/sqr.log" : (ix == 1) ? "logs/sqr_kara.log" : "logs/sqr_toom.log", "w");
       for (cnt = 4; cnt <= 10240 / DIGIT_BIT; cnt += 2) {
 	 SLEEP;
 	 mp_rand(&a, cnt);
@@ -194,13 +215,12 @@ int main(void)
 	 } while (++rr < 100);
 	 printf("Squaring\t%4d-bit => %9llu/sec, %9llu cycles\n",
 		mp_count_bits(&a), CLK_PER_SEC / tt, tt);
-	 fprintf(log, "%d %9llu\n", mp_count_bits(&a), tt);
-	 fflush(log);
+	 FPRINTF(log, "%d %9llu\n", mp_count_bits(&a), tt);
+	 FFLUSH(log);
       }
-      fclose(log);
+      FCLOSE(log);
 
    }
- exptmod:
 
    {
       char *primes[] = {
@@ -235,10 +255,10 @@ int main(void)
 	 "1214855636816562637502584060163403830270705000634713483015101384881871978446801224798536155406895823305035467591632531067547890948695117172076954220727075688048751022421198712032848890056357845974246560748347918630050853933697792254955890439720297560693579400297062396904306270145886830719309296352765295712183040773146419022875165382778007040109957609739589875590885701126197906063620133954893216612678838507540777138437797705602453719559017633986486649523611975865005712371194067612263330335590526176087004421363598470302731349138773205901447704682181517904064735636518462452242791676541725292378925568296858010151852326316777511935037531017413910506921922450666933202278489024521263798482237150056835746454842662048692127173834433089016107854491097456725016327709663199738238442164843147132789153725513257167915555162094970853584447993125488607696008169807374736711297007473812256272245489405898470297178738029484459690836250560495461579533254473316340608217876781986188705928270735695752830825527963838355419762516246028680280988020401914551825487349990306976304093109384451438813251211051597392127491464898797406789175453067960072008590614886532333015881171367104445044718144312416815712216611576221546455968770801413440778423979",
 	 NULL
       };
-      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");
+      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);
@@ -271,17 +291,17 @@ int main(void)
 	 }
 	 printf("Exponentiating\t%4d-bit => %9llu/sec, %9llu cycles\n",
 		mp_count_bits(&a), CLK_PER_SEC / tt, tt);
-	 fprintf(n < 4 ? logd : (n < 9) ? logc : (n < 16) ? logb : log,
+	 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);
+   FCLOSE(log);
+   FCLOSE(logb);
+   FCLOSE(logc);
+   FCLOSE(logd);
 
-   log = fopen("logs/invmod.log", "w");
-   for (cnt = 4; cnt <= 128; cnt += 4) {
+   log = FOPEN("logs/invmod.log", "w");
+   for (cnt = 4; cnt <= 32; cnt += 4) {
       SLEEP;
       mp_rand(&a, cnt);
       mp_rand(&b, cnt);
@@ -307,13 +327,13 @@ int main(void)
       }
       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);
+      FPRINTF(log, "%d %9llu\n", cnt * DIGIT_BIT, tt);
    }
-   fclose(log);
+   FCLOSE(log);
 
    return 0;
 }
 
 /* $Source$ */
-/* $Revision: 0.36 $ */
-/* $Date: 2005-08-01 16:37:28 +0000 $ */
+/* $Revision$ */
+/* $Date$ */