Take over recent improvements from libtommath after version 1.1.0.

17 files changed, 87 insertions, 101 deletions

diff --git a/generic/tclTomMath.h b/generic/tclTomMath.h
index 3f23fd6..26eef26 100644
--- a/generic/tclTomMath.h
+++ b/generic/tclTomMath.h
@@ -7,8 +7,7 @@
  * 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.
+ * SPDX-License-Identifier: Unlicense
  */
 #ifndef BN_H_
 #define BN_H_
@@ -128,6 +127,7 @@ typedef unsigned long long   mp_word;
 #define MP_MEM        -2  /* out of mem */
 #define MP_VAL        -3  /* invalid input */
 #define MP_RANGE      MP_VAL
+#define MP_ITER       -4  /* Max. iterations reached */
 
 #define MP_YES        1   /* yes response */
 #define MP_NO         0   /* no response */
@@ -346,15 +346,20 @@ int mp_cnt_lsb(const mp_int *a);
 
 /* I Love Earth! */
 
-/* makes a pseudo-random int of a given size */
+/* makes a pseudo-random mp_int of a given size */
 /*
 int mp_rand(mp_int *a, int digits);
 */
+/* makes a pseudo-random small int of a given size */
+/*
+int mp_rand_digit(mp_digit *r);
+*/
 
 #ifdef MP_PRNG_ENABLE_LTM_RNG
-/* as last resort we will fall back to libtomcrypt's rng_get_bytes()
- * in case you don't use libtomcrypt or use it w/o rng_get_bytes()
- * you have to implement it somewhere else, as it's required */
+/* A last resort to provide random data on systems without any of the other
+ * implemented ways to gather entropy.
+ * It is compatible with `rng_get_bytes()` from libtomcrypt so you could
+ * provide that one and then set `ltm_rng = rng_get_bytes;` */
 extern unsigned long (*ltm_rng)(unsigned char *out, unsigned long outlen, void (*callback)(void));
 extern void (*ltm_rng_callback)(void);
 #endif
@@ -691,10 +696,17 @@ int mp_prime_miller_rabin(const mp_int *a, const mp_int *b, int *result);
 int mp_prime_rabin_miller_trials(int size);
 */
 
-/* performs t rounds of Miller-Rabin on "a" using the first
- * t prime bases.  Also performs an initial sieve of trial
+/* performs t random rounds of Miller-Rabin on "a" additional to
+ * bases 2 and 3.  Also performs an initial sieve of trial
  * division.  Determines if "a" is prime with probability
  * of error no more than (1/4)**t.
+ * Both a strong Lucas-Selfridge to complete the BPSW test
+ * and a separate Frobenius test are available at compile time.
+ * With t<0 a deterministic test is run for primes up to
+ * 318665857834031151167461. With t<13 (abs(t)-13) additional
+ * tests with sequential small primes are run starting at 43.
+ * Is Fips 186.4 compliant if called with t as computed by
+ * mp_prime_rabin_miller_trials();
  *
  * Sets result to 1 if probably prime, 0 otherwise
  */
diff --git a/generic/tclTomMathDecls.h b/generic/tclTomMathDecls.h
index 9b525d2..df66748 100644
--- a/generic/tclTomMathDecls.h
+++ b/generic/tclTomMathDecls.h
@@ -37,9 +37,10 @@
 /* MODULE_SCOPE void  TclBNFree( void* ); */
 #define TclBNFree(x) (ckfree((char*)(x)))
 
-#define XMALLOC(x) TclBNAlloc(x)
-#define XFREE(x) TclBNFree(x)
-#define XREALLOC(x,n) TclBNRealloc(x,n)
+#define XMALLOC(size)                   TclBNAlloc(size)
+#define XFREE(mem, size)                TclBNFree(mem)
+#define XREALLOC(mem, oldsize, newsize) TclBNRealloc(mem, newsize)
+
 
 /* Rename the global symbols in libtommath to avoid linkage conflicts */
 
diff --git a/libtommath/bn_mp_clear.c b/libtommath/bn_mp_clear.c
index 1f360b2..b8e724c 100644
--- a/libtommath/bn_mp_clear.c
+++ b/libtommath/bn_mp_clear.c
@@ -25,7 +25,7 @@ void mp_clear(mp_int *a)
       }
 
       /* free ram */
-      XFREE(a->dp);
+      XFREE(a->dp, sizeof (mp_digit) * (size_t)a->alloc);
 
       /* reset members to make debugging easier */
       a->dp    = NULL;
diff --git a/libtommath/bn_mp_fwrite.c b/libtommath/bn_mp_fwrite.c
index 9f0c3df..85a942f 100644
--- a/libtommath/bn_mp_fwrite.c
+++ b/libtommath/bn_mp_fwrite.c
@@ -22,24 +22,24 @@ int mp_fwrite(const mp_int *a, int radix, FILE *stream)
       return err;
    }
 
-   buf = OPT_CAST(char) XMALLOC((size_t)len);
+   buf = (char *) XMALLOC((size_t)len);
    if (buf == NULL) {
       return MP_MEM;
    }
 
    if ((err = mp_toradix(a, buf, radix)) != MP_OKAY) {
-      XFREE(buf);
+      XFREE(buf, len);
       return err;
    }
 
    for (x = 0; x < len; x++) {
       if (fputc((int)buf[x], stream) == EOF) {
-         XFREE(buf);
+         XFREE(buf, len);
          return MP_VAL;
       }
    }
 
-   XFREE(buf);
+   XFREE(buf, len);
    return MP_OKAY;
 }
 #endif
diff --git a/libtommath/bn_mp_get_double.c b/libtommath/bn_mp_get_double.c
index 3ed5a71..629eae3 100644
--- a/libtommath/bn_mp_get_double.c
+++ b/libtommath/bn_mp_get_double.c
@@ -19,10 +19,10 @@ double mp_get_double(const mp_int *a)
    for (i = 0; i < DIGIT_BIT; ++i) {
       fac *= 2.0;
    }
-   for (i = USED(a); i --> 0;) {
-      d = (d * fac) + (double)DIGIT(a, i);
+   for (i = a->used; i --> 0;) {
+      d = (d * fac) + (double)a->dp[i];
    }
-   return (mp_isneg(a) != MP_NO) ? -d : d;
+   return (a->sign == MP_NEG) ? -d : d;
 }
 #endif
 
diff --git a/libtommath/bn_mp_get_int.c b/libtommath/bn_mp_get_int.c
index 13eddbf..d9c7a11 100644
--- a/libtommath/bn_mp_get_int.c
+++ b/libtommath/bn_mp_get_int.c
@@ -15,25 +15,8 @@
 /* get the lower 32-bits of an mp_int */
 unsigned long mp_get_int(const mp_int *a)
 {
-   int i;
-   mp_min_u32 res;
-
-   if (a->used == 0) {
-      return 0;
-   }
-
-   /* get number of digits of the lsb we have to read */
-   i = MIN(a->used, ((((int)sizeof(unsigned long) * CHAR_BIT) + DIGIT_BIT - 1) / DIGIT_BIT)) - 1;
-
-   /* get most significant digit of result */
-   res = DIGIT(a, i);
-
-   while (--i >= 0) {
-      res = (res << DIGIT_BIT) | DIGIT(a, i);
-   }
-
    /* force result to 32-bits always so it is consistent on non 32-bit platforms */
-   return res & 0xFFFFFFFFUL;
+   return mp_get_long(a) & 0xFFFFFFFFUL;
 }
 #endif
 
diff --git a/libtommath/bn_mp_get_long.c b/libtommath/bn_mp_get_long.c
index a4d05d6..b95bb8a 100644
--- a/libtommath/bn_mp_get_long.c
+++ b/libtommath/bn_mp_get_long.c
@@ -18,19 +18,19 @@ unsigned long mp_get_long(const mp_int *a)
    int i;
    unsigned long res;
 
-   if (a->used == 0) {
+   if (IS_ZERO(a)) {
       return 0;
    }
 
    /* get number of digits of the lsb we have to read */
-   i = MIN(a->used, ((((int)sizeof(unsigned long) * CHAR_BIT) + DIGIT_BIT - 1) / DIGIT_BIT)) - 1;
+   i = MIN(a->used, (((CHAR_BIT * (int)sizeof(unsigned long)) + DIGIT_BIT - 1) / DIGIT_BIT)) - 1;
 
    /* get most significant digit of result */
-   res = DIGIT(a, i);
+   res = (unsigned long)a->dp[i];
 
-#if (ULONG_MAX != 0xffffffffuL) || (DIGIT_BIT < 32)
+#if (ULONG_MAX != 0xFFFFFFFFUL) || (DIGIT_BIT < 32)
    while (--i >= 0) {
-      res = (res << DIGIT_BIT) | DIGIT(a, i);
+      res = (res << DIGIT_BIT) | (unsigned long)a->dp[i];
    }
 #endif
    return res;
diff --git a/libtommath/bn_mp_get_long_long.c b/libtommath/bn_mp_get_long_long.c
index 4201b4d..cafd9a4 100644
--- a/libtommath/bn_mp_get_long_long.c
+++ b/libtommath/bn_mp_get_long_long.c
@@ -18,19 +18,19 @@ unsigned long long mp_get_long_long(const mp_int *a)
    int i;
    unsigned long long res;
 
-   if (a->used == 0) {
+   if (IS_ZERO(a)) {
       return 0;
    }
 
    /* get number of digits of the lsb we have to read */
-   i = MIN(a->used, ((((int)sizeof(unsigned long long) * CHAR_BIT) + DIGIT_BIT - 1) / DIGIT_BIT)) - 1;
+   i = MIN(a->used, (((CHAR_BIT * (int)sizeof(unsigned long long)) + DIGIT_BIT - 1) / DIGIT_BIT)) - 1;
 
    /* get most significant digit of result */
-   res = DIGIT(a, i);
+   res = (unsigned long long)a->dp[i];
 
 #if DIGIT_BIT < 64
    while (--i >= 0) {
-      res = (res << DIGIT_BIT) | DIGIT(a, i);
+      res = (res << DIGIT_BIT) | (unsigned long long)a->dp[i];
    }
 #endif
    return res;
diff --git a/libtommath/bn_mp_grow.c b/libtommath/bn_mp_grow.c
index 1d92b29..b120194 100644
--- a/libtommath/bn_mp_grow.c
+++ b/libtommath/bn_mp_grow.c
@@ -29,7 +29,9 @@ int mp_grow(mp_int *a, int size)
        * in case the operation failed we don't want
        * to overwrite the dp member of a.
        */
-      tmp = OPT_CAST(mp_digit) XREALLOC(a->dp, sizeof(mp_digit) * (size_t)size);
+      tmp = (mp_digit *) XREALLOC(a->dp,
+                                  (size_t)a->alloc * sizeof (mp_digit),
+                                  (size_t)size * sizeof(mp_digit));
       if (tmp == NULL) {
          /* reallocation failed but "a" is still valid [can be freed] */
          return MP_MEM;
diff --git a/libtommath/bn_mp_init.c b/libtommath/bn_mp_init.c
index 7520089..3c0c489 100644
--- a/libtommath/bn_mp_init.c
+++ b/libtommath/bn_mp_init.c
@@ -18,7 +18,7 @@ int mp_init(mp_int *a)
    int i;
 
    /* allocate memory required and clear it */
-   a->dp = OPT_CAST(mp_digit) XMALLOC(sizeof(mp_digit) * (size_t)MP_PREC);
+   a->dp = (mp_digit *) XMALLOC(MP_PREC * sizeof(mp_digit));
    if (a->dp == NULL) {
       return MP_MEM;
    }
diff --git a/libtommath/bn_mp_init_size.c b/libtommath/bn_mp_init_size.c
index 9b933fb..1becb23 100644
--- a/libtommath/bn_mp_init_size.c
+++ b/libtommath/bn_mp_init_size.c
@@ -21,7 +21,7 @@ int mp_init_size(mp_int *a, int size)
    size += (MP_PREC * 2) - (size % MP_PREC);
 
    /* alloc mem */
-   a->dp = OPT_CAST(mp_digit) XMALLOC(sizeof(mp_digit) * (size_t)size);
+   a->dp = (mp_digit *) XMALLOC((size_t)size * sizeof(mp_digit));
    if (a->dp == NULL) {
       return MP_MEM;
    }
diff --git a/libtommath/bn_mp_is_square.c b/libtommath/bn_mp_is_square.c
index 5363a47..1dd1d6c 100644
--- a/libtommath/bn_mp_is_square.c
+++ b/libtommath/bn_mp_is_square.c
@@ -49,13 +49,12 @@ int mp_is_square(const mp_int *arg, int *ret)
       return MP_VAL;
    }
 
-   /* digits used?  (TSD) */
-   if (arg->used == 0) {
+   if (IS_ZERO(arg)) {
       return MP_OKAY;
    }
 
    /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
-   if (rem_128[127u & DIGIT(arg, 0)] == (char)1) {
+   if (rem_128[127u & arg->dp[0]] == (char)1) {
       return MP_OKAY;
    }
 
diff --git a/libtommath/bn_mp_prime_random_ex.c b/libtommath/bn_mp_prime_random_ex.c
index b0b4632..0ca29ec 100644
--- a/libtommath/bn_mp_prime_random_ex.c
+++ b/libtommath/bn_mp_prime_random_ex.c
@@ -46,19 +46,19 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
    bsize = (size>>3) + ((size&7)?1:0);
 
    /* we need a buffer of bsize bytes */
-   tmp = OPT_CAST(unsigned char) XMALLOC((size_t)bsize);
+   tmp = (unsigned char *) XMALLOC((size_t)bsize);
    if (tmp == NULL) {
       return MP_MEM;
    }
 
    /* calc the maskAND value for the MSbyte*/
-   maskAND = ((size&7) == 0) ? 0xFF : (0xFF >> (8 - (size & 7)));
+   maskAND = ((size&7) == 0) ? 0xFF : (unsigned char)(0xFF >> (8 - (size & 7)));
 
    /* calc the maskOR_msb */
    maskOR_msb        = 0;
    maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
    if ((flags & LTM_PRIME_2MSB_ON) != 0) {
-      maskOR_msb       |= 0x80 >> ((9 - size) & 7);
+      maskOR_msb       |= (unsigned char)(0x80 >> ((9 - size) & 7));
    }
 
    /* get the maskOR_lsb */
@@ -76,7 +76,7 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
 
       /* work over the MSbyte */
       tmp[0]    &= maskAND;
-      tmp[0]    |= 1 << ((size - 1) & 7);
+      tmp[0]    |= (unsigned char)(1 << ((size - 1) & 7));
 
       /* mix in the maskORs */
       tmp[maskOR_msb_offset]   |= maskOR_msb;
@@ -123,7 +123,7 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
 
    err = MP_OKAY;
 error:
-   XFREE(tmp);
+   XFREE(tmp, bsize);
    return err;
 }
 
diff --git a/libtommath/bn_mp_read_radix.c b/libtommath/bn_mp_read_radix.c
index 200601e..a8723b7 100644
--- a/libtommath/bn_mp_read_radix.c
+++ b/libtommath/bn_mp_read_radix.c
@@ -12,6 +12,8 @@
  * SPDX-License-Identifier: Unlicense
  */
 
+#define MP_TOUPPER(c) ((((c) >= 'a') && ((c) <= 'z')) ? (((c) + 'A') - 'a') : (c))
+
 /* read a string [ASCII] in a given radix */
 int mp_read_radix(mp_int *a, const char *str, int radix)
 {
@@ -46,7 +48,7 @@ int mp_read_radix(mp_int *a, const char *str, int radix)
        * this allows numbers like 1AB and 1ab to represent the same  value
        * [e.g. in hex]
        */
-      ch = (radix <= 36) ? (char)toupper((int)*str) : *str;
+      ch = (radix <= 36) ? (char)MP_TOUPPER((int)*str) : *str;
       pos = (unsigned)(ch - '(');
       if (mp_s_rmap_reverse_sz < pos) {
          break;
diff --git a/libtommath/bn_mp_set_double.c b/libtommath/bn_mp_set_double.c
index 76f6293..c96a3b3 100644
--- a/libtommath/bn_mp_set_double.c
+++ b/libtommath/bn_mp_set_double.c
@@ -41,8 +41,8 @@ int mp_set_double(mp_int *a, double b)
       return res;
    }
 
-   if (((cast.bits >> 63) != 0ULL) && (mp_iszero(a) == MP_NO)) {
-      SIGN(a) = MP_NEG;
+   if (((cast.bits >> 63) != 0ULL) && !IS_ZERO(a)) {
+      a->sign = MP_NEG;
    }
 
    return MP_OKAY;
diff --git a/libtommath/bn_mp_shrink.c b/libtommath/bn_mp_shrink.c
index ff7905f..fa30184 100644
--- a/libtommath/bn_mp_shrink.c
+++ b/libtommath/bn_mp_shrink.c
@@ -23,7 +23,9 @@ int mp_shrink(mp_int *a)
    }
 
    if (a->alloc != used) {
-      if ((tmp = OPT_CAST(mp_digit) XREALLOC(a->dp, sizeof(mp_digit) * (size_t)used)) == NULL) {
+      if ((tmp = (mp_digit *) XREALLOC(a->dp,
+                                       (size_t)a->alloc * sizeof (mp_digit),
+                                       (size_t)used * sizeof(mp_digit))) == NULL) {
          return MP_MEM;
       }
       a->dp    = tmp;
diff --git a/libtommath/tommath_private.h b/libtommath/tommath_private.h
index 2096f77..b3600a0 100644
--- a/libtommath/tommath_private.h
+++ b/libtommath/tommath_private.h
@@ -13,7 +13,6 @@
 #define TOMMATH_PRIV_H_
 
 #include <tommath.h>
-#include <ctype.h>
 
 #ifndef MIN
 #define MIN(x, y) (((x) < (y)) ? (x) : (y))
@@ -25,28 +24,19 @@
 
 #ifdef __cplusplus
 extern "C" {
-
-/* C++ compilers don't like assigning void * to mp_digit * */
-#define OPT_CAST(x) (x *)
-
-#else
-
-/* C on the other hand doesn't care */
-#define OPT_CAST(x)
-
 #endif
 
 /* define heap macros */
 #ifndef XMALLOC
 /* default to libc stuff */
-#   define XMALLOC   malloc
-#   define XFREE     free
-#   define XREALLOC  realloc
+#   define XMALLOC(size)                   malloc(size)
+#   define XFREE(mem, size)                free(mem)
+#   define XREALLOC(mem, oldsize, newsize) realloc(mem, newsize)
 #elif 0
 /* prototypes for our heap functions */
-extern void *XMALLOC(size_t n);
-extern void *XREALLOC(void *p, size_t n);
-extern void XFREE(void *p);
+extern void *XMALLOC(size_t size);
+extern void *XREALLOC(void *mem, size_t oldsize, size_t newsize);
+extern void XFREE(void *mem, size_t size);
 #endif
 
 /* you'll have to tune these... */
@@ -58,6 +48,11 @@ extern void XFREE(void *p);
 /* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */
 #define MP_WARRAY               (1u << (((sizeof(mp_word) * CHAR_BIT) - (2 * DIGIT_BIT)) + 1))
 
+/* ---> Basic Manipulations <--- */
+#define IS_ZERO(a) ((a)->used == 0)
+#define IS_EVEN(a) (((a)->used == 0) || (((a)->dp[0] & 1u) == 0u))
+#define IS_ODD(a)  (((a)->used > 0) && (((a)->dp[0] & 1u) == 1u))
+
 /* lowlevel functions, do not call! */
 int s_mp_add(const mp_int *a, const mp_int *b, mp_int *c);
 int s_mp_sub(const mp_int *a, const mp_int *b, mp_int *c);
@@ -85,36 +80,26 @@ extern const size_t mp_s_rmap_reverse_sz;
 
 /* Fancy macro to set an MPI from another type.
  * There are several things assumed:
- *  x is the counter and unsigned
+ *  x is the counter
  *  a is the pointer to the MPI
  *  b is the original value that should be set in the MPI.
  */
 #define MP_SET_XLONG(func_name, type)                    \
 int func_name (mp_int * a, type b)                       \
 {                                                        \
-  unsigned int  x;                                       \
-  int           res;                                     \
-                                                         \
-  mp_zero (a);                                           \
-                                                         \
-  /* set four bits at a time */                          \
-  for (x = 0; x < (sizeof(type) * 2u); x++) {            \
-    /* shift the number up four bits */                  \
-    if ((res = mp_mul_2d (a, 4, a)) != MP_OKAY) {        \
-      return res;                                        \
-    }                                                    \
-                                                         \
-    /* OR in the top four bits of the source */          \
-    a->dp[0] |= (mp_digit)(b >> ((sizeof(type) * 8u) - 4u)) & 15uL;\
-                                                         \
-    /* shift the source up to the next four bits */      \
-    b <<= 4;                                             \
-                                                         \
-    /* ensure that digits are not clamped off */         \
-    a->used += 1;                                        \
-  }                                                      \
-  mp_clamp (a);                                          \
-  return MP_OKAY;                                        \
+   int x = 0;                                            \
+   int new_size = (((CHAR_BIT * sizeof(type)) + DIGIT_BIT) - 1) / DIGIT_BIT; \
+   int res = mp_grow(a, new_size);                       \
+   if (res == MP_OKAY) {                                 \
+     mp_zero(a);                                         \
+     while (b != 0u) {                                   \
+        a->dp[x++] = ((mp_digit)b & MP_MASK);            \
+        if ((CHAR_BIT * sizeof (b)) <= DIGIT_BIT) { break; } \
+        b >>= ((CHAR_BIT * sizeof (b)) <= DIGIT_BIT ? 0 : DIGIT_BIT); \
+     }                                                   \
+     a->used = x;                                        \
+   }                                                     \
+   return res;                                           \
 }
 
 #ifdef __cplusplus