Update some libtommath functions to the latest trunk versions. Small step forward in the upgrade to (upcoming) libtommath 1.2.

Advantage: simplify Tcl code accessing those functions.

7 files changed, 422 insertions, 382 deletions

diff --git a/libtommath/bn_mp_and.c b/libtommath/bn_mp_and.c
index 02bef18..54c0b4e 100644
--- a/libtommath/bn_mp_and.c
+++ b/libtommath/bn_mp_and.c
@@ -1,53 +1,56 @@
 #include <tommath.h>
 #ifdef BN_MP_AND_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
-
-/* AND two ints together */
-int
-mp_and (mp_int * a, mp_int * b, mp_int * c)
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* two complement and */
+mp_err mp_and(const mp_int *a, const mp_int *b, mp_int *c)
 {
-  int     res, ix, px;
-  mp_int  t, *x;
-
-  if (a->used > b->used) {
-    if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
-      return res;
-    }
-    px = b->used;
-    x = b;
-  } else {
-    if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
-      return res;
-    }
-    px = a->used;
-    x = a;
-  }
-
-  for (ix = 0; ix < px; ix++) {
-    t.dp[ix] &= x->dp[ix];
-  }
-
-  /* zero digits above the last from the smallest mp_int */
-  for (; ix < t.used; ix++) {
-    t.dp[ix] = 0;
-  }
-
-  mp_clamp (&t);
-  mp_exch (c, &t);
-  mp_clear (&t);
-  return MP_OKAY;
+   int used = MAX(a->used, b->used) + 1, i;
+   mp_err err;
+   mp_digit ac = 1, bc = 1, cc = 1;
+   mp_sign csign = ((a->sign == MP_NEG) && (b->sign == MP_NEG)) ? MP_NEG : MP_ZPOS;
+
+   if (c->alloc < used) {
+      if ((err = mp_grow(c, used)) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   for (i = 0; i < used; i++) {
+      mp_digit x, y;
+
+      /* convert to two complement if negative */
+      if (a->sign == MP_NEG) {
+         ac += (i >= a->used) ? MP_MASK : (~a->dp[i] & MP_MASK);
+         x = ac & MP_MASK;
+         ac >>= MP_DIGIT_BIT;
+      } else {
+         x = (i >= a->used) ? 0uL : a->dp[i];
+      }
+
+      /* convert to two complement if negative */
+      if (b->sign == MP_NEG) {
+         bc += (i >= b->used) ? MP_MASK : (~b->dp[i] & MP_MASK);
+         y = bc & MP_MASK;
+         bc >>= MP_DIGIT_BIT;
+      } else {
+         y = (i >= b->used) ? 0uL : b->dp[i];
+      }
+
+      c->dp[i] = x & y;
+
+      /* convert to to sign-magnitude if negative */
+      if (csign == MP_NEG) {
+         cc += ~c->dp[i] & MP_MASK;
+         c->dp[i] = cc & MP_MASK;
+         cc >>= MP_DIGIT_BIT;
+      }
+   }
+
+   c->used = used;
+   c->sign = csign;
+   mp_clamp(c);
+   return MP_OKAY;
 }
 #endif
diff --git a/libtommath/bn_mp_cmp.c b/libtommath/bn_mp_cmp.c
index b965d4b..c042b63 100644
--- a/libtommath/bn_mp_cmp.c
+++ b/libtommath/bn_mp_cmp.c
@@ -1,39 +1,26 @@
 #include <tommath.h>
 #ifdef BN_MP_CMP_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
 
 /* compare two ints (signed)*/
-int
-mp_cmp (mp_int * a, mp_int * b)
+mp_ord mp_cmp(const mp_int *a, const mp_int *b)
 {
-  /* compare based on sign */
-  if (a->sign != b->sign) {
-     if (a->sign == MP_NEG) {
-        return MP_LT;
-     } else {
-        return MP_GT;
-     }
-  }
-  
-  /* compare digits */
-  if (a->sign == MP_NEG) {
-     /* if negative compare opposite direction */
-     return mp_cmp_mag(b, a);
-  } else {
-     return mp_cmp_mag(a, b);
-  }
+   /* compare based on sign */
+   if (a->sign != b->sign) {
+      if (a->sign == MP_NEG) {
+         return MP_LT;
+      } else {
+         return MP_GT;
+      }
+   }
+
+   /* compare digits */
+   if (a->sign == MP_NEG) {
+      /* if negative compare opposite direction */
+      return mp_cmp_mag(b, a);
+   } else {
+      return mp_cmp_mag(a, b);
+   }
 }
 #endif
diff --git a/libtommath/bn_mp_cmp_d.c b/libtommath/bn_mp_cmp_d.c
index a446bb4..947c57a 100644
--- a/libtommath/bn_mp_cmp_d.c
+++ b/libtommath/bn_mp_cmp_d.c
@@ -1,40 +1,28 @@
 #include <tommath.h>
 #ifdef BN_MP_CMP_D_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
 
 /* compare a digit */
-int mp_cmp_d(mp_int * a, mp_digit b)
+mp_ord mp_cmp_d(const mp_int *a, mp_digit b)
 {
-  /* compare based on sign */
-  if (a->sign == MP_NEG) {
-    return MP_LT;
-  }
+   /* compare based on sign */
+   if (a->sign == MP_NEG) {
+      return MP_LT;
+   }
 
-  /* compare based on magnitude */
-  if (a->used > 1) {
-    return MP_GT;
-  }
+   /* compare based on magnitude */
+   if (a->used > 1) {
+      return MP_GT;
+   }
 
-  /* compare the only digit of a to b */
-  if (a->dp[0] > b) {
-    return MP_GT;
-  } else if (a->dp[0] < b) {
-    return MP_LT;
-  } else {
-    return MP_EQ;
-  }
+   /* compare the only digit of a to b */
+   if (a->dp[0] > b) {
+      return MP_GT;
+   } else if (a->dp[0] < b) {
+      return MP_LT;
+   } else {
+      return MP_EQ;
+   }
 }
 #endif
diff --git a/libtommath/bn_mp_cmp_mag.c b/libtommath/bn_mp_cmp_mag.c
index 3506d2b..850e083 100644
--- a/libtommath/bn_mp_cmp_mag.c
+++ b/libtommath/bn_mp_cmp_mag.c
@@ -1,51 +1,39 @@
 #include <tommath.h>
 #ifdef BN_MP_CMP_MAG_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
 
 /* compare maginitude of two ints (unsigned) */
-int mp_cmp_mag (mp_int * a, mp_int * b)
+mp_ord mp_cmp_mag(const mp_int *a, const mp_int *b)
 {
-  int     n;
-  mp_digit *tmpa, *tmpb;
+   int     n;
+   const mp_digit *tmpa, *tmpb;
 
-  /* compare based on # of non-zero digits */
-  if (a->used > b->used) {
-    return MP_GT;
-  }
-  
-  if (a->used < b->used) {
-    return MP_LT;
-  }
+   /* compare based on # of non-zero digits */
+   if (a->used > b->used) {
+      return MP_GT;
+   }
 
-  /* alias for a */
-  tmpa = a->dp + (a->used - 1);
+   if (a->used < b->used) {
+      return MP_LT;
+   }
 
-  /* alias for b */
-  tmpb = b->dp + (a->used - 1);
+   /* alias for a */
+   tmpa = a->dp + (a->used - 1);
 
-  /* compare based on digits  */
-  for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
-    if (*tmpa > *tmpb) {
-      return MP_GT;
-    }
+   /* alias for b */
+   tmpb = b->dp + (a->used - 1);
 
-    if (*tmpa < *tmpb) {
-      return MP_LT;
-    }
-  }
-  return MP_EQ;
+   /* compare based on digits  */
+   for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
+      if (*tmpa > *tmpb) {
+         return MP_GT;
+      }
+
+      if (*tmpa < *tmpb) {
+         return MP_LT;
+      }
+   }
+   return MP_EQ;
 }
 #endif
diff --git a/libtommath/bn_mp_or.c b/libtommath/bn_mp_or.c
index aa5b1bd..afcdd9b 100644
--- a/libtommath/bn_mp_or.c
+++ b/libtommath/bn_mp_or.c
@@ -1,46 +1,56 @@
 #include <tommath.h>
 #ifdef BN_MP_OR_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
-
-/* OR two ints together */
-int mp_or (mp_int * a, mp_int * b, mp_int * c)
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* two complement or */
+mp_err mp_or(const mp_int *a, const mp_int *b, mp_int *c)
 {
-  int     res, ix, px;
-  mp_int  t, *x;
-
-  if (a->used > b->used) {
-    if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
-      return res;
-    }
-    px = b->used;
-    x = b;
-  } else {
-    if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
-      return res;
-    }
-    px = a->used;
-    x = a;
-  }
-
-  for (ix = 0; ix < px; ix++) {
-    t.dp[ix] |= x->dp[ix];
-  }
-  mp_clamp (&t);
-  mp_exch (c, &t);
-  mp_clear (&t);
-  return MP_OKAY;
+   int used = MAX(a->used, b->used) + 1, i;
+   mp_err err;
+   mp_digit ac = 1, bc = 1, cc = 1;
+   mp_sign csign = ((a->sign == MP_NEG) || (b->sign == MP_NEG)) ? MP_NEG : MP_ZPOS;
+
+   if (c->alloc < used) {
+      if ((err = mp_grow(c, used)) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   for (i = 0; i < used; i++) {
+      mp_digit x, y;
+
+      /* convert to two complement if negative */
+      if (a->sign == MP_NEG) {
+         ac += (i >= a->used) ? MP_MASK : (~a->dp[i] & MP_MASK);
+         x = ac & MP_MASK;
+         ac >>= MP_DIGIT_BIT;
+      } else {
+         x = (i >= a->used) ? 0uL : a->dp[i];
+      }
+
+      /* convert to two complement if negative */
+      if (b->sign == MP_NEG) {
+         bc += (i >= b->used) ? MP_MASK : (~b->dp[i] & MP_MASK);
+         y = bc & MP_MASK;
+         bc >>= MP_DIGIT_BIT;
+      } else {
+         y = (i >= b->used) ? 0uL : b->dp[i];
+      }
+
+      c->dp[i] = x | y;
+
+      /* convert to to sign-magnitude if negative */
+      if (csign == MP_NEG) {
+         cc += ~c->dp[i] & MP_MASK;
+         c->dp[i] = cc & MP_MASK;
+         cc >>= MP_DIGIT_BIT;
+      }
+   }
+
+   c->used = used;
+   c->sign = csign;
+   mp_clamp(c);
+   return MP_OKAY;
 }
 #endif
diff --git a/libtommath/bn_mp_xor.c b/libtommath/bn_mp_xor.c
index 432f42e..fba6617 100644
--- a/libtommath/bn_mp_xor.c
+++ b/libtommath/bn_mp_xor.c
@@ -1,47 +1,56 @@
 #include <tommath.h>
 #ifdef BN_MP_XOR_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
-
-/* XOR two ints together */
-int
-mp_xor (mp_int * a, mp_int * b, mp_int * c)
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* two complement xor */
+mp_err mp_xor(const mp_int *a, const mp_int *b, mp_int *c)
 {
-  int     res, ix, px;
-  mp_int  t, *x;
-
-  if (a->used > b->used) {
-    if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
-      return res;
-    }
-    px = b->used;
-    x = b;
-  } else {
-    if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
-      return res;
-    }
-    px = a->used;
-    x = a;
-  }
-
-  for (ix = 0; ix < px; ix++) {
-     t.dp[ix] ^= x->dp[ix];
-  }
-  mp_clamp (&t);
-  mp_exch (c, &t);
-  mp_clear (&t);
-  return MP_OKAY;
+   int used = MAX(a->used, b->used) + 1, i;
+   mp_err err;
+   mp_digit ac = 1, bc = 1, cc = 1;
+   mp_sign csign = (a->sign != b->sign) ? MP_NEG : MP_ZPOS;
+
+   if (c->alloc < used) {
+      if ((err = mp_grow(c, used)) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   for (i = 0; i < used; i++) {
+      mp_digit x, y;
+
+      /* convert to two complement if negative */
+      if (a->sign == MP_NEG) {
+         ac += (i >= a->used) ? MP_MASK : (~a->dp[i] & MP_MASK);
+         x = ac & MP_MASK;
+         ac >>= MP_DIGIT_BIT;
+      } else {
+         x = (i >= a->used) ? 0uL : a->dp[i];
+      }
+
+      /* convert to two complement if negative */
+      if (b->sign == MP_NEG) {
+         bc += (i >= b->used) ? MP_MASK : (~b->dp[i] & MP_MASK);
+         y = bc & MP_MASK;
+         bc >>= MP_DIGIT_BIT;
+      } else {
+         y = (i >= b->used) ? 0uL : b->dp[i];
+      }
+
+      c->dp[i] = x ^ y;
+
+      /* convert to to sign-magnitude if negative */
+      if (csign == MP_NEG) {
+         cc += ~c->dp[i] & MP_MASK;
+         c->dp[i] = cc & MP_MASK;
+         cc >>= MP_DIGIT_BIT;
+      }
+   }
+
+   c->used = used;
+   c->sign = csign;
+   mp_clamp(c);
+   return MP_OKAY;
 }
 #endif
diff --git a/libtommath/tommath.h b/libtommath/tommath.h
index 49b2c2a..df460f6 100644
--- a/libtommath/tommath.h
+++ b/libtommath/tommath.h
@@ -1,28 +1,14 @@
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
- *
- * LibTomMath is a library that provides multiple-precision
- * integer arithmetic as well as number theoretic functionality.
- *
- * The library was designed directly after the MPI library by
- * Michael Fromberger but has been written from scratch with
- * additional optimizations in place.
- *
- * The library is free for all purposes without any express
- * guarantee it works.
- *
- * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
- */
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
 #ifndef BN_H_
 #define BN_H_
 
-#include <stdio.h>
 #include <string.h>
 #include <stdlib.h>
 #include <ctype.h>
 #include <limits.h>
 
-#include <tommath_class.h>
-
 #ifndef MIN
    #define MIN(x,y) ((x)<(y)?(x):(y))
 #endif
@@ -31,6 +17,10 @@
    #define MAX(x,y) ((x)>(y)?(x):(y))
 #endif
 
+#ifndef MP_NO_FILE
+#  include <stdio.h>
+#endif
+
 #ifdef __cplusplus
 extern "C" {
 
@@ -131,45 +121,74 @@ extern "C" {
 #define MP_MASK          ((((mp_digit)1)<<((mp_digit)DIGIT_BIT))-((mp_digit)1))
 #define MP_DIGIT_MAX     MP_MASK
 
-/* equalities */
-#define MP_LT        -1   /* less than */
-#define MP_EQ         0   /* equal to */
-#define MP_GT         1   /* greater than */
+/* Primality generation flags */
+#define LTM_PRIME_BBS      0x0001 /* BBS style prime */
+#define LTM_PRIME_SAFE     0x0002 /* Safe prime (p-1)/2 == prime */
+#define LTM_PRIME_2MSB_ON  0x0008 /* force 2nd MSB to 1 */
 
+#ifdef MP_USE_ENUMS
+typedef enum {
+   MP_ZPOS = 0,
+   MP_NEG = 1
+} mp_sign;
+typedef enum {
+   MP_LT = -1,
+   MP_EQ = 0,
+   MP_GT = 1
+} mp_ord;
+typedef enum {
+   MP_NO = 0,
+   MP_YES = 1
+} mp_bool;
+typedef enum {
+   MP_OKAY  = 0,
+   MP_ERR   = -1,
+   MP_MEM   = -2,
+   MP_VAL   = -3,
+   MP_ITER  = -4
+} mp_err;
+#else
+typedef int mp_sign;
 #define MP_ZPOS       0   /* positive integer */
 #define MP_NEG        1   /* negative */
-
+typedef int mp_ord;
+#define MP_LT        -1   /* less than */
+#define MP_EQ         0   /* equal to */
+#define MP_GT         1   /* greater than */
+typedef int mp_bool;
+#define MP_YES        1   /* yes response */
+#define MP_NO         0   /* no response */
+typedef int mp_err;
 #define MP_OKAY       0   /* ok result */
+#define MP_ERR        -1  /* unknown error */
 #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 */
+#endif
 
-#define MP_YES        1   /* yes response */
-#define MP_NO         0   /* no response */
-
-/* Primality generation flags */
-#define LTM_PRIME_BBS      0x0001 /* BBS style prime */
-#define LTM_PRIME_SAFE     0x0002 /* Safe prime (p-1)/2 == prime */
-#define LTM_PRIME_2MSB_ON  0x0008 /* force 2nd MSB to 1 */
-
-typedef int           mp_err;
+/* tunable cutoffs */
 
-/* you'll have to tune these... */
-extern int KARATSUBA_MUL_CUTOFF,
-           KARATSUBA_SQR_CUTOFF,
-           TOOM_MUL_CUTOFF,
-           TOOM_SQR_CUTOFF;
+#ifndef MP_FIXED_CUTOFFS
+extern int
+KARATSUBA_MUL_CUTOFF,
+KARATSUBA_SQR_CUTOFF,
+TOOM_MUL_CUTOFF,
+TOOM_SQR_CUTOFF;
+#endif
 
 /* define this to use lower memory usage routines (exptmods mostly) */
 /* #define MP_LOW_MEM */
 
 /* default precision */
 #ifndef MP_PREC
-   #ifndef MP_LOW_MEM
-      #define MP_PREC                 32     /* default digits of precision */
-   #else
-      #define MP_PREC                 8      /* default digits of precision */
-   #endif   
+#   ifndef MP_LOW_MEM
+#      define MP_PREC 32        /* default digits of precision */
+#   elif defined(MP_8BIT)
+#      define MP_PREC 16        /* default digits of precision */
+#   else
+#      define MP_PREC 8         /* default digits of precision */
+#   endif
 #endif
 
 /* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */
@@ -177,8 +196,9 @@ extern int KARATSUBA_MUL_CUTOFF,
 
 /* the infamous mp_int structure */
 typedef struct  {
-    int used, alloc, sign;
-    mp_digit *dp;
+   int used, alloc;
+   mp_sign sign;
+   mp_digit *dp;
 } mp_int;
 
 /* callback for mp_prime_random, should fill dst with random bytes and return how many read [upto len] */
@@ -190,17 +210,17 @@ typedef int ltm_prime_callback(unsigned char *dst, int len, void *dat);
 #define SIGN(m)    ((m)->sign)
 
 /* error code to char* string */
-char *mp_error_to_string(int code);
+const char *mp_error_to_string(mp_err code);
 
 /* ---> init and deinit bignum functions <--- */
 /* init a bignum */
-int mp_init(mp_int *a);
+mp_err mp_init(mp_int *a);
 
 /* free a bignum */
 void mp_clear(mp_int *a);
 
 /* init a null terminated series of arguments */
-int mp_init_multi(mp_int *mp, ...);
+mp_err mp_init_multi(mp_int *mp, ...);
 
 /* clear a null terminated series of arguments */
 void mp_clear_multi(mp_int *mp, ...);
@@ -209,18 +229,19 @@ void mp_clear_multi(mp_int *mp, ...);
 void mp_exch(mp_int *a, mp_int *b);
 
 /* shrink ram required for a bignum */
-int mp_shrink(mp_int *a);
+mp_err mp_shrink(mp_int *a);
 
 /* grow an int to a given size */
-int mp_grow(mp_int *a, int size);
+mp_err mp_grow(mp_int *a, int size);
 
 /* init to a given number of digits */
-int mp_init_size(mp_int *a, int size);
+mp_err mp_init_size(mp_int *a, int size);
 
 /* ---> Basic Manipulations <--- */
 #define mp_iszero(a) (((a)->used == 0) ? MP_YES : MP_NO)
 #define mp_iseven(a) (((a)->used == 0 || (((a)->dp[0] & 1) == 0)) ? MP_YES : MP_NO)
 #define mp_isodd(a)  (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? MP_YES : MP_NO)
+#define mp_isneg(a)  (((a)->sign != MP_ZPOS) ? MP_YES : MP_NO)
 
 /* set to zero */
 void mp_zero(mp_int *a);
@@ -241,141 +262,153 @@ int mp_init_set (mp_int * a, mp_digit b);
 int mp_init_set_int (mp_int * a, unsigned long b);
 
 /* copy, b = a */
-int mp_copy(mp_int *a, mp_int *b);
+mp_err mp_copy(const mp_int *a, mp_int *b);
 
 /* inits and copies, a = b */
-int mp_init_copy(mp_int *a, mp_int *b);
+mp_err mp_init_copy(mp_int *a, const mp_int *b);
 
 /* trim unused digits */
 void mp_clamp(mp_int *a);
 
+/* import binary data */
+mp_err mp_import(mp_int *rop, size_t count, int order, size_t size, int endian, size_t nails, const void *op);
+
+/* export binary data */
+mp_err mp_export(void *rop, size_t *countp, int order, size_t size, int endian, size_t nails, const mp_int *op);
+
 /* ---> digit manipulation <--- */
 
 /* right shift by "b" digits */
 void mp_rshd(mp_int *a, int b);
 
 /* left shift by "b" digits */
-int mp_lshd(mp_int *a, int b);
+mp_err mp_lshd(mp_int *a, int b);
 
-/* c = a / 2**b */
-int mp_div_2d(mp_int *a, int b, mp_int *c, mp_int *d);
+/* c = a / 2**b, implemented as c = a >> b */
+mp_err mp_div_2d(const mp_int *a, int b, mp_int *c, mp_int *d);
 
 /* b = a/2 */
-int mp_div_2(mp_int *a, mp_int *b);
+mp_err mp_div_2(const mp_int *a, mp_int *b);
 
-/* c = a * 2**b */
-int mp_mul_2d(mp_int *a, int b, mp_int *c);
+/* c = a * 2**b, implemented as c = a << b */
+mp_err mp_mul_2d(const mp_int *a, int b, mp_int *c);
 
 /* b = a*2 */
-int mp_mul_2(mp_int *a, mp_int *b);
+mp_err mp_mul_2(const mp_int *a, mp_int *b);
 
-/* c = a mod 2**d */
-int mp_mod_2d(mp_int *a, int b, mp_int *c);
+/* c = a mod 2**b */
+mp_err mp_mod_2d(const mp_int *a, int b, mp_int *c);
 
 /* computes a = 2**b */
-int mp_2expt(mp_int *a, int b);
+mp_err mp_2expt(mp_int *a, int b);
 
 /* Counts the number of lsbs which are zero before the first zero bit */
-int mp_cnt_lsb(mp_int *a);
+int mp_cnt_lsb(const mp_int *a);
 
 /* I Love Earth! */
 
 /* makes a pseudo-random int of a given size */
-int mp_rand(mp_int *a, int digits);
+mp_err mp_rand(mp_int *a, int digits);
 
 /* ---> binary operations <--- */
-/* c = a XOR b  */
-int mp_xor(mp_int *a, mp_int *b, mp_int *c);
+/* c = a XOR b (two complement) */
+mp_err mp_xor(const mp_int *a, const mp_int *b, mp_int *c);
 
-/* c = a OR b */
-int mp_or(mp_int *a, mp_int *b, mp_int *c);
+/* c = a OR b (two complement) */
+mp_err mp_or(const mp_int *a, const mp_int *b, mp_int *c);
 
-/* c = a AND b */
-int mp_and(mp_int *a, mp_int *b, mp_int *c);
+/* c = a AND b (two complement) */
+mp_err mp_and(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* ---> Basic arithmetic <--- */
 
 /* b = -a */
-int mp_neg(mp_int *a, mp_int *b);
+mp_err mp_neg(const mp_int *a, mp_int *b);
 
 /* b = |a| */
-int mp_abs(mp_int *a, mp_int *b);
+mp_err mp_abs(const mp_int *a, mp_int *b);
 
 /* compare a to b */
-int mp_cmp(mp_int *a, mp_int *b);
+mp_ord mp_cmp(const mp_int *a, const mp_int *b);
 
 /* compare |a| to |b| */
-int mp_cmp_mag(mp_int *a, mp_int *b);
+mp_ord mp_cmp_mag(const mp_int *a, const mp_int *b);
 
 /* c = a + b */
-int mp_add(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_add(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = a - b */
-int mp_sub(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_sub(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = a * b */
-int mp_mul(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_mul(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* b = a*a  */
-int mp_sqr(mp_int *a, mp_int *b);
+mp_err mp_sqr(const mp_int *a, mp_int *b);
 
 /* a/b => cb + d == a */
-int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d);
 
 /* c = a mod b, 0 <= c < b  */
-int mp_mod(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_mod(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* ---> single digit functions <--- */
 
 /* compare against a single digit */
-int mp_cmp_d(mp_int *a, mp_digit b);
+mp_ord mp_cmp_d(const mp_int *a, mp_digit b);
 
 /* c = a + b */
-int mp_add_d(mp_int *a, mp_digit b, mp_int *c);
+mp_err mp_add_d(const mp_int *a, mp_digit b, mp_int *c);
+
+/* Increment "a" by one like "a++". Changes input! */
+mp_err mp_incr(mp_int *a);
 
 /* c = a - b */
-int mp_sub_d(mp_int *a, mp_digit b, mp_int *c);
+mp_err mp_sub_d(const mp_int *a, mp_digit b, mp_int *c);
+
+/* Decrement "a" by one like "a--". Changes input! */
+mp_err mp_decr(mp_int *a);
 
 /* c = a * b */
-int mp_mul_d(mp_int *a, mp_digit b, mp_int *c);
+mp_err mp_mul_d(const mp_int *a, mp_digit b, mp_int *c);
 
 /* a/b => cb + d == a */
-int mp_div_d(mp_int *a, mp_digit b, mp_int *c, mp_digit *d);
+mp_err mp_div_d(const mp_int *a, mp_digit b, mp_int *c, mp_digit *d);
 
 /* a/3 => 3c + d == a */
-int mp_div_3(mp_int *a, mp_int *c, mp_digit *d);
+mp_err mp_div_3(const mp_int *a, mp_int *c, mp_digit *d);
 
 /* c = a**b */
-int mp_expt_d(mp_int *a, mp_digit b, mp_int *c);
+mp_err mp_expt_d(const mp_int *a, mp_digit b, mp_int *c);
 
 /* c = a mod b, 0 <= c < b  */
-int mp_mod_d(mp_int *a, mp_digit b, mp_digit *c);
+mp_err mp_mod_d(const mp_int *a, mp_digit b, mp_digit *c);
 
 /* ---> number theory <--- */
 
 /* d = a + b (mod c) */
-int mp_addmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+mp_err mp_addmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* d = a - b (mod c) */
-int mp_submod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+mp_err mp_submod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* d = a * b (mod c) */
-int mp_mulmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+mp_err mp_mulmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* c = a * a (mod b) */
-int mp_sqrmod(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_sqrmod(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = 1/a (mod b) */
-int mp_invmod(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_invmod(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = (a, b) */
-int mp_gcd(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_gcd(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* produces value such that U1*a + U2*b = U3 */
-int mp_exteuclid(mp_int *a, mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3);
+mp_err mp_exteuclid(const mp_int *a, const mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3);
 
 /* c = [a, b] or (a*b)/(a, b) */
-int mp_lcm(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_lcm(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* finds one of the b'th root of a, such that |c|**b <= |a|
  *
@@ -384,64 +417,64 @@ int mp_lcm(mp_int *a, mp_int *b, mp_int *c);
 int mp_n_root(mp_int *a, mp_digit b, mp_int *c);
 
 /* special sqrt algo */
-int mp_sqrt(mp_int *arg, mp_int *ret);
+mp_err mp_sqrt(const mp_int *arg, mp_int *ret);
 
 /* is number a square? */
-int mp_is_square(mp_int *arg, int *ret);
+mp_err mp_is_square(const mp_int *arg, mp_bool *ret);
 
 /* computes the jacobi c = (a | n) (or Legendre if b is prime)  */
 int mp_jacobi(mp_int *a, mp_int *n, int *c);
 
 /* used to setup the Barrett reduction for a given modulus b */
-int mp_reduce_setup(mp_int *a, mp_int *b);
+mp_err mp_reduce_setup(mp_int *a, const mp_int *b);
 
 /* Barrett Reduction, computes a (mod b) with a precomputed value c
  *
- * Assumes that 0 < a <= b*b, note if 0 > a > -(b*b) then you can merely
- * compute the reduction as -1 * mp_reduce(mp_abs(a)) [pseudo code].
+ * Assumes that 0 < x <= m*m, note if 0 > x > -(m*m) then you can merely
+ * compute the reduction as -1 * mp_reduce(mp_abs(x)) [pseudo code].
  */
-int mp_reduce(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_reduce(mp_int *x, const mp_int *m, const mp_int *mu);
 
 /* setups the montgomery reduction */
-int mp_montgomery_setup(mp_int *a, mp_digit *mp);
+mp_err mp_montgomery_setup(const mp_int *n, mp_digit *rho);
 
 /* computes a = B**n mod b without division or multiplication useful for
  * normalizing numbers in a Montgomery system.
  */
-int mp_montgomery_calc_normalization(mp_int *a, mp_int *b);
+mp_err mp_montgomery_calc_normalization(mp_int *a, const mp_int *b);
 
 /* computes x/R == x (mod N) via Montgomery Reduction */
-int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
+mp_err mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho);
 
 /* returns 1 if a is a valid DR modulus */
-int mp_dr_is_modulus(mp_int *a);
+mp_bool mp_dr_is_modulus(const mp_int *a);
 
 /* sets the value of "d" required for mp_dr_reduce */
-void mp_dr_setup(mp_int *a, mp_digit *d);
+void mp_dr_setup(const mp_int *a, mp_digit *d);
 
-/* reduces a modulo b using the Diminished Radix method */
-int mp_dr_reduce(mp_int *a, mp_int *b, mp_digit mp);
+/* reduces a modulo n using the Diminished Radix method */
+mp_err mp_dr_reduce(mp_int *x, const mp_int *n, mp_digit k);
 
 /* returns true if a can be reduced with mp_reduce_2k */
-int mp_reduce_is_2k(mp_int *a);
+mp_bool mp_reduce_is_2k(const mp_int *a);
 
 /* determines k value for 2k reduction */
-int mp_reduce_2k_setup(mp_int *a, mp_digit *d);
+mp_err mp_reduce_2k_setup(const mp_int *a, mp_digit *d);
 
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
-int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d);
+mp_err mp_reduce_2k(mp_int *a, const mp_int *n, mp_digit d);
 
 /* returns true if a can be reduced with mp_reduce_2k_l */
-int mp_reduce_is_2k_l(mp_int *a);
+mp_bool mp_reduce_is_2k_l(const mp_int *a);
 
 /* determines k value for 2k reduction */
-int mp_reduce_2k_setup_l(mp_int *a, mp_int *d);
+mp_err mp_reduce_2k_setup_l(const mp_int *a, mp_int *d);
 
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
-int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d);
+mp_err mp_reduce_2k_l(mp_int *a, const mp_int *n, const mp_int *d);
 
-/* d = a**b (mod c) */
-int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+/* Y = G**X (mod P) */
+mp_err mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y);
 
 /* ---> Primes <--- */
 
@@ -456,41 +489,58 @@ int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
 extern const mp_digit ltm_prime_tab[];
 
 /* result=1 if a is divisible by one of the first PRIME_SIZE primes */
-int mp_prime_is_divisible(mp_int *a, int *result);
+mp_err mp_prime_is_divisible(mp_int *a, mp_bool *result);
 
 /* performs one Fermat test of "a" using base "b".
  * Sets result to 0 if composite or 1 if probable prime
  */
-int mp_prime_fermat(mp_int *a, mp_int *b, int *result);
+mp_err mp_prime_fermat(const mp_int *a, const mp_int *b, mp_bool *result);
 
 /* performs one Miller-Rabin test of "a" using base "b".
  * Sets result to 0 if composite or 1 if probable prime
  */
-int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result);
+mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, mp_bool *result);
 
 /* This gives [for a given bit size] the number of trials required
- * such that Miller-Rabin gives a prob of failure lower than 2^-96 
+ * such that Miller-Rabin gives a prob of failure lower than 2^-96
  */
 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 one strong Lucas-Selfridge test of "a".
+ * Sets result to 0 if composite or 1 if probable prime
+ */
+mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, mp_bool *result);
+
+/* performs one Frobenius test of "a" as described by Paul Underwood.
+ * Sets result to 0 if composite or 1 if probable prime
+ */
+mp_err mp_prime_frobenius_underwood(const mp_int *N, mp_bool *result);
+
+/* 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
  */
-int mp_prime_is_prime(mp_int *a, int t, int *result);
+mp_err mp_prime_is_prime(const mp_int *a, int t, mp_bool *result);
 
 /* finds the next prime after the number "a" using "t" trials
  * of Miller-Rabin.
  *
  * bbs_style = 1 means the prime must be congruent to 3 mod 4
  */
-int mp_prime_next_prime(mp_int *a, int t, int bbs_style);
+mp_err mp_prime_next_prime(mp_int *a, int t, int bbs_style);
 
 /* makes a truly random prime of a given size (bytes),
- * call with bbs = 1 if you want it to be congruent to 3 mod 4 
+ * call with bbs = 1 if you want it to be congruent to 3 mod 4
  *
  * You have to supply a callback which fills in a buffer with random bytes.  "dat" is a parameter you can
  * have passed to the callback (e.g. a state or something).  This function doesn't use "dat" itself
@@ -503,38 +553,43 @@ int mp_prime_next_prime(mp_int *a, int t, int bbs_style);
 /* makes a truly random prime of a given size (bits),
  *
  * Flags are as follows:
- * 
- *   LTM_PRIME_BBS      - make prime congruent to 3 mod 4
- *   LTM_PRIME_SAFE     - make sure (p-1)/2 is prime as well (implies LTM_PRIME_BBS)
- *   LTM_PRIME_2MSB_ON  - make the 2nd highest bit one
+ *
+ *   MP_PRIME_BBS      - make prime congruent to 3 mod 4
+ *   MP_PRIME_SAFE     - make sure (p-1)/2 is prime as well (implies MP_PRIME_BBS)
+ *   MP_PRIME_2MSB_ON  - make the 2nd highest bit one
  *
  * You have to supply a callback which fills in a buffer with random bytes.  "dat" is a parameter you can
  * have passed to the callback (e.g. a state or something).  This function doesn't use "dat" itself
  * so it can be NULL
  *
  */
-int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat);
+mp_err mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat);
 
-/* ---> radix conversion <--- */
-int mp_count_bits(mp_int *a);
-
-int mp_unsigned_bin_size(mp_int *a);
-int mp_read_unsigned_bin(mp_int *a, const unsigned char *b, int c);
-int mp_to_unsigned_bin(mp_int *a, unsigned char *b);
-int mp_to_unsigned_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen);
-
-int mp_signed_bin_size(mp_int *a);
-int mp_read_signed_bin(mp_int *a, const unsigned char *b, int c);
-int mp_to_signed_bin(mp_int *a,  unsigned char *b);
-int mp_to_signed_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen);
+/* Integer logarithm to integer base */
+mp_err mp_ilogb(const mp_int *a, mp_digit base, mp_int *c);
 
-int mp_read_radix(mp_int *a, const char *str, int radix);
-int mp_toradix(mp_int *a, char *str, int radix);
-int mp_toradix_n(mp_int * a, char *str, int radix, int maxlen);
-int mp_radix_size(mp_int *a, int radix, int *size);
-
-int mp_fread(mp_int *a, int radix, FILE *stream);
-int mp_fwrite(mp_int *a, int radix, FILE *stream);
+/* ---> radix conversion <--- */
+int mp_count_bits(const mp_int *a);
+
+int mp_unsigned_bin_size(const mp_int *a);
+mp_err mp_read_unsigned_bin(mp_int *a, const unsigned char *b, int c);
+mp_err mp_to_unsigned_bin(const mp_int *a, unsigned char *b);
+mp_err mp_to_unsigned_bin_n(const mp_int * a, unsigned char *b, unsigned long *outlen);
+
+int mp_signed_bin_size(const mp_int *a);
+mp_err mp_read_signed_bin(mp_int *a, const unsigned char *b, int c);
+mp_err mp_to_signed_bin(const mp_int *a,  unsigned char *b);
+mp_err mp_to_signed_bin_n(const mp_int * a, unsigned char *b, unsigned long *outlen);
+
+mp_err mp_read_radix(mp_int *a, const char *str, int radix);
+mp_err mp_toradix(const mp_int *a, char *str, int radix);
+mp_err mp_toradix_n(const mp_int * a, char *str, int radix, int maxlen);
+mp_err mp_radix_size(const mp_int *a, int radix, int *size);
+
+#ifndef MP_NO_FILE
+mp_err mp_fread(mp_int *a, int radix, FILE *stream);
+mp_err mp_fwrite(const mp_int *a, int radix, FILE *stream);
+#endif
 
 #define mp_read_raw(mp, str, len) mp_read_signed_bin((mp), (str), (len))
 #define mp_raw_size(mp)           mp_signed_bin_size(mp)