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.

16 files changed, 603 insertions, 684 deletions

diff --git a/generic/tclBasic.c b/generic/tclBasic.c
index b148333..3b9fca9 100644
--- a/generic/tclBasic.c
+++ b/generic/tclBasic.c
@@ -1920,7 +1920,7 @@ Tcl_CreateCommand(
 
 	/* An existing command conflicts. Try to delete it.. */
 	cmdPtr = Tcl_GetHashValue(hPtr);
-	
+
 	/*
 	 * Be careful to preserve
 	 * any existing import links so we can restore them down below. That
@@ -2107,7 +2107,7 @@ Tcl_CreateObjCommand(
 	cmdPtr = Tcl_GetHashValue(hPtr);
 
 	/*
-	 * [***] This is wrong.  See Tcl Bug a16752c252.  
+	 * [***] This is wrong.  See Tcl Bug a16752c252.
 	 * However, this buggy behavior is kept under particular
 	 * circumstances to accommodate deployed binaries of the
 	 * "tclcompiler" program. http://sourceforge.net/projects/tclpro/
@@ -5188,7 +5188,7 @@ TclEvalObjEx(
 		TclStackAlloc(interp, sizeof(CmdFrame));
 
 	eoFramePtr->type = TCL_LOCATION_EVAL_LIST;
-	eoFramePtr->level = (iPtr->cmdFramePtr == NULL?  1 
+	eoFramePtr->level = (iPtr->cmdFramePtr == NULL?  1
 		: iPtr->cmdFramePtr->level + 1);
 	eoFramePtr->framePtr = iPtr->framePtr;
 	eoFramePtr->nextPtr = iPtr->cmdFramePtr;
@@ -6346,7 +6346,7 @@ ExprIsqrtFunc(
 	if (Tcl_GetBignumFromObj(interp, objv[1], &big) != TCL_OK) {
 	    return TCL_ERROR;
 	}
-	if (SIGN(&big) == MP_NEG) {
+	if (mp_isneg(&big)) {
 	    mp_clear(&big);
 	    goto negarg;
 	}
@@ -6617,8 +6617,7 @@ ExprAbsFunc(
 #endif
 
     if (type == TCL_NUMBER_BIG) {
-	/* TODO: const correctness ? */
-	if (mp_cmp_d((mp_int *) ptr, 0) == MP_LT) {
+	if (mp_cmp_d(ptr, 0) == MP_LT) {
 	    Tcl_GetBignumFromObj(NULL, objv[1], &big);
 	tooLarge:
 	    mp_neg(&big, &big);
@@ -6768,7 +6767,7 @@ ExprIntFunc(
 	mp_int big;
 
 	Tcl_GetBignumFromObj(NULL, objPtr, &big);
-	mp_mod_2d(&big, (int) CHAR_BIT * sizeof(long), &big);
+	mp_mod_2d(&big, CHAR_BIT * sizeof(long), &big);
 	objPtr = Tcl_NewBignumObj(&big);
 	Tcl_IncrRefCount(objPtr);
 	TclGetLongFromObj(NULL, objPtr, &iResult);
@@ -6800,7 +6799,7 @@ ExprWideFunc(
 	mp_int big;
 
 	Tcl_GetBignumFromObj(NULL, objPtr, &big);
-	mp_mod_2d(&big, (int) CHAR_BIT * sizeof(Tcl_WideInt), &big);
+	mp_mod_2d(&big, CHAR_BIT * sizeof(Tcl_WideInt), &big);
 	objPtr = Tcl_NewBignumObj(&big);
 	Tcl_IncrRefCount(objPtr);
 	Tcl_GetWideIntFromObj(NULL, objPtr, &wResult);
@@ -7007,7 +7006,7 @@ ExprSrandFunc(
 	    return TCL_ERROR;
 	}
 
-	mp_mod_2d(&big, (int) CHAR_BIT * sizeof(long), &big);
+	mp_mod_2d(&big, CHAR_BIT * sizeof(long), &big);
 	objPtr = Tcl_NewBignumObj(&big);
 	Tcl_IncrRefCount(objPtr);
 	TclGetLongFromObj(NULL, objPtr, &i);
diff --git a/generic/tclExecute.c b/generic/tclExecute.c
index 61d0ddc..265b82f 100644
--- a/generic/tclExecute.c
+++ b/generic/tclExecute.c
@@ -1361,7 +1361,7 @@ FreeExprCodeInternalRep(
 
 ByteCode *
 TclCompileObj(
-    Tcl_Interp *interp, 
+    Tcl_Interp *interp,
     Tcl_Obj *objPtr,
     const CmdFrame *invoker,
     int word)
@@ -4942,29 +4942,28 @@ TclExecuteByteCode(
 
 		    Tcl_TakeBignumFromObj(NULL, value2Ptr, &big2);
 
-		    /* TODO: internals intrusion */
-		    if ((l1 > 0) ^ (big2.sign == MP_ZPOS)) {
+		    if ((l1 > 0) ^ mp_isneg(&big2)) {
 			/*
-			 * Arguments are opposite sign; remainder is sum.
+			 * Arguments are same sign; remainder is first operand.
 			 */
 
-			mp_int big1;
-
-			TclBNInitBignumFromLong(&big1, l1);
-			mp_add(&big2, &big1, &big2);
-			mp_clear(&big1);
-			objResultPtr = Tcl_NewBignumObj(&big2);
-			TRACE(("%s\n", O2S(objResultPtr)));
-			NEXT_INST_F(1, 2, 1);
+			mp_clear(&big2);
+			TRACE(("%s\n", O2S(valuePtr)));
+			NEXT_INST_F(1, 1, 0);
 		    }
 
 		    /*
-		     * Arguments are same sign; remainder is first operand.
+		     * Arguments are opposite sign; remainder is sum.
 		     */
 
-		    mp_clear(&big2);
-		    TRACE(("%s\n", O2S(valuePtr)));
-		    NEXT_INST_F(1, 1, 0);
+		    mp_int big1;
+
+		    TclBNInitBignumFromLong(&big1, l1);
+		    mp_add(&big2, &big1, &big2);
+		    mp_clear(&big1);
+		    objResultPtr = Tcl_NewBignumObj(&big2);
+		    TRACE(("%s\n", O2S(objResultPtr)));
+		    NEXT_INST_F(1, 2, 1);
 		}
 	    }
 #ifndef NO_WIDE_TYPE
@@ -4999,28 +4998,28 @@ TclExecuteByteCode(
 		    Tcl_TakeBignumFromObj(NULL, value2Ptr, &big2);
 
 		    /* TODO: internals intrusion */
-		    if ((w1 > ((Tcl_WideInt) 0)) ^ (big2.sign == MP_ZPOS)) {
+		    if ((w1 > ((Tcl_WideInt) 0)) ^ mp_isneg(&big2)) {
 			/*
-			 * Arguments are opposite sign; remainder is sum.
+			 * Arguments are same sign; remainder is first operand.
 			 */
 
-			mp_int big1;
-
-			TclBNInitBignumFromWideInt(&big1, w1);
-			mp_add(&big2, &big1, &big2);
-			mp_clear(&big1);
-			objResultPtr = Tcl_NewBignumObj(&big2);
-			TRACE(("%s\n", O2S(objResultPtr)));
-			NEXT_INST_F(1, 2, 1);
+			mp_clear(&big2);
+			TRACE(("%s\n", O2S(valuePtr)));
+			NEXT_INST_F(1, 1, 0);
 		    }
-
 		    /*
-		     * Arguments are same sign; remainder is first operand.
+		     * Arguments are opposite sign; remainder is sum.
 		     */
 
-		    mp_clear(&big2);
-		    TRACE(("%s\n", O2S(valuePtr)));
-		    NEXT_INST_F(1, 1, 0);
+		    mp_int big1;
+
+		    TclBNInitBignumFromWideInt(&big1, w1);
+		    mp_add(&big2, &big1, &big2);
+		    mp_clear(&big1);
+		    objResultPtr = Tcl_NewBignumObj(&big2);
+		    TRACE(("%s\n", O2S(objResultPtr)));
+		    NEXT_INST_F(1, 2, 1);
+
 		}
 	    }
 #endif
@@ -5033,7 +5032,7 @@ TclExecuteByteCode(
 		mp_init(&bigRemainder);
 		mp_div(&big1, &big2, &bigResult, &bigRemainder);
 		if (!mp_iszero(&bigRemainder)
-			&& (bigRemainder.sign != big2.sign)) {
+			&& (mp_isneg(&bigRemainder) != mp_isneg(&big2))) {
 		    /*
 		     * Convert to Tcl's integer division rules.
 		     */
@@ -5316,139 +5315,24 @@ TclExecuteByteCode(
 	}
 
 	if ((type1 == TCL_NUMBER_BIG) || (type2 == TCL_NUMBER_BIG)) {
-	    mp_int big1, big2, bigResult, *First, *Second;
-	    int numPos;
+	    mp_int big1, big2, bigResult;
 
 	    Tcl_TakeBignumFromObj(NULL, valuePtr, &big1);
 	    Tcl_TakeBignumFromObj(NULL, value2Ptr, &big2);
 
-	    /*
-	     * Count how many positive arguments we have. If only one of the
-	     * arguments is negative, store it in 'Second'.
-	     */
-
-	    if (mp_cmp_d(&big1, 0) != MP_LT) {
-		numPos = 1 + (mp_cmp_d(&big2, 0) != MP_LT);
-		First = &big1;
-		Second = &big2;
-	    } else {
-		First = &big2;
-		Second = &big1;
-		numPos = (mp_cmp_d(First, 0) != MP_LT);
-	    }
 	    mp_init(&bigResult);
 
 	    switch (*pc) {
 	    case INST_BITAND:
-		switch (numPos) {
-		case 2:
-		    /*
-		     * Both arguments positive, base case.
-		     */
-
-		    mp_and(First, Second, &bigResult);
-		    break;
-		case 1:
-		    /*
-		     * First is positive; second negative:
-		     * P & N = P & ~~N = P&~(-N-1) = P & (P ^ (-N-1))
-		     */
-
-		    mp_neg(Second, Second);
-		    mp_sub_d(Second, 1, Second);
-		    mp_xor(First, Second, &bigResult);
-		    mp_and(First, &bigResult, &bigResult);
-		    break;
-		case 0:
-		    /*
-		     * Both arguments negative:
-		     * a & b = ~ (~a | ~b) = -(-a-1|-b-1)-1
-		     */
-
-		    mp_neg(First, First);
-		    mp_sub_d(First, 1, First);
-		    mp_neg(Second, Second);
-		    mp_sub_d(Second, 1, Second);
-		    mp_or(First, Second, &bigResult);
-		    mp_neg(&bigResult, &bigResult);
-		    mp_sub_d(&bigResult, 1, &bigResult);
-		    break;
-		}
+		mp_and(&big1, &big2, &bigResult);
 		break;
 
 	    case INST_BITOR:
-		switch (numPos) {
-		case 2:
-		    /*
-		     * Both arguments positive, base case.
-		     */
-
-		    mp_or(First, Second, &bigResult);
-		    break;
-		case 1:
-		    /*
-		     * First is positive; second negative:
-		     * N|P = ~(~N&~P) = ~((-N-1)&~P) = -((-N-1)&((-N-1)^P))-1
-		     */
-
-		    mp_neg(Second, Second);
-		    mp_sub_d(Second, 1, Second);
-		    mp_xor(First, Second, &bigResult);
-		    mp_and(Second, &bigResult, &bigResult);
-		    mp_neg(&bigResult, &bigResult);
-		    mp_sub_d(&bigResult, 1, &bigResult);
-		    break;
-		case 0:
-		    /*
-		     * Both arguments negative:
-		     * a | b = ~ (~a & ~b) = -(-a-1&-b-1)-1
-		     */
-
-		    mp_neg(First, First);
-		    mp_sub_d(First, 1, First);
-		    mp_neg(Second, Second);
-		    mp_sub_d(Second, 1, Second);
-		    mp_and(First, Second, &bigResult);
-		    mp_neg(&bigResult, &bigResult);
-		    mp_sub_d(&bigResult, 1, &bigResult);
-		    break;
-		}
+		mp_or(&big1, &big2, &bigResult);
 		break;
 
 	    case INST_BITXOR:
-		switch (numPos) {
-		case 2:
-		    /*
-		     * Both arguments positive, base case.
-		     */
-
-		    mp_xor(First, Second, &bigResult);
-		    break;
-		case 1:
-		    /*
-		     * First is positive; second negative:
-		     * P^N = ~(P^~N) = -(P^(-N-1))-1
-		     */
-
-		    mp_neg(Second, Second);
-		    mp_sub_d(Second, 1, Second);
-		    mp_xor(First, Second, &bigResult);
-		    mp_neg(&bigResult, &bigResult);
-		    mp_sub_d(&bigResult, 1, &bigResult);
-		    break;
-		case 0:
-		    /*
-		     * Both arguments negative:
-		     * a ^ b = (~a ^ ~b) = (-a-1^-b-1)
-		     */
-
-		    mp_neg(First, First);
-		    mp_sub_d(First, 1, First);
-		    mp_neg(Second, Second);
-		    mp_sub_d(Second, 1, Second);
-		    mp_xor(First, Second, &bigResult);
-		    break;
-		}
+		mp_xor(&big1, &big2, &bigResult);
 		break;
 	    }
 
@@ -6256,9 +6140,8 @@ TclExecuteByteCode(
 		}
 		mp_init(&bigRemainder);
 		mp_div(&big1, &big2, &bigResult, &bigRemainder);
-		/* TODO: internals intrusion */
 		if (!mp_iszero(&bigRemainder)
-			&& (bigRemainder.sign != big2.sign)) {
+			&& (mp_isneg(&bigRemainder) != mp_isneg(&big2))) {
 		    /*
 		     * Convert to Tcl's integer division rules.
 		     */
diff --git a/generic/tclObj.c b/generic/tclObj.c
index 1738985..6bff71c 100644
--- a/generic/tclObj.c
+++ b/generic/tclObj.c
@@ -189,7 +189,7 @@ static Tcl_ThreadDataKey pendingObjDataKey;
 	    mp_shrink(&(bignum)); \
 	} \
 	(objPtr)->internalRep.ptrAndLongRep.ptr = (void*) (bignum).dp; \
-	(objPtr)->internalRep.ptrAndLongRep.value = ( ((bignum).sign << 30) \
+	(objPtr)->internalRep.ptrAndLongRep.value = ( (mp_isneg(&bignum) << 30) \
 		| ((bignum).alloc << 15) | ((bignum).used)); \
     }
 
@@ -2787,7 +2787,7 @@ Tcl_GetLongFromObj(
 		    while (numBytes-- > 0) {
 			value = (value << CHAR_BIT) | *bytes++;
 		    }
-		    if (big.sign) {
+		    if (mp_isneg(&big)) {
 			*longPtr = - (long) value;
 		    } else {
 			*longPtr = (long) value;
@@ -3089,7 +3089,7 @@ Tcl_GetWideIntFromObj(
 		    while (numBytes-- > 0) {
 			value = (value << CHAR_BIT) | *bytes++;
 		    }
-		    if (big.sign) {
+		    if (mp_isneg(&big)) {
 			*wideIntPtr = - (Tcl_WideInt) value;
 		    } else {
 			*wideIntPtr = (Tcl_WideInt) value;
@@ -3508,10 +3508,10 @@ Tcl_SetBignumObj(
 	while (numBytes-- > 0) {
 	    value = (value << CHAR_BIT) | *bytes++;
 	}
-	if (value > (((~(unsigned long)0) >> 1) + bignumValue->sign)) {
+	if (value > (((~(unsigned long)0) >> 1) + mp_isneg(bignumValue))) {
 	    goto tooLargeForLong;
 	}
-	if (bignumValue->sign) {
+	if (mp_isneg(bignumValue)) {
 	    TclSetLongObj(objPtr, -(long)value);
 	} else {
 	    TclSetLongObj(objPtr, (long)value);
@@ -3533,10 +3533,10 @@ Tcl_SetBignumObj(
 	while (numBytes-- > 0) {
 	    value = (value << CHAR_BIT) | *bytes++;
 	}
-	if (value > (((~(Tcl_WideUInt)0) >> 1) + bignumValue->sign)) {
+	if (value > (((~(Tcl_WideUInt)0) >> 1) + mp_isneg(bignumValue))) {
 	    goto tooLargeForWide;
 	}
-	if (bignumValue->sign) {
+	if (mp_isneg(bignumValue)) {
 	    TclSetWideIntObj(objPtr, -(Tcl_WideInt)value);
 	} else {
 	    TclSetWideIntObj(objPtr, (Tcl_WideInt)value);
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index 3ed4349..17d630b 100644
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -129,7 +129,7 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
 				/* Highest power of two that is greater than
 				 * DBL_MAX_10_EXP, divided by 16 */
 #define DIGIT_GROUP		8
-				/* floor(DIGIT_BIT*log(2)/log(10)) */
+				/* floor(MP_DIGIT_BIT*log(2)/log(10)) */
 
 /* Union used to dismantle floating point numbers. */
 
@@ -1447,9 +1447,9 @@ AccumulateDecimalDigit(
 	 * More than single digit multiplication. Multiply by the appropriate
 	 * small powers of 5, and then shift. Large strings of zeroes are
 	 * eaten 256 at a time; this is less efficient than it could be, but
-	 * seems implausible. We presume that DIGIT_BIT is at least 27. The
+	 * seems implausible. We presume that MP_DIGIT_BIT is at least 27. The
 	 * first multiplication, by up to 10**7, is done with a one-DIGIT
-	 * multiply (this presumes that DIGIT_BIT >= 24).
+	 * multiply (this presumes that MP_DIGIT_BIT >= 24).
 	 */
 
 	n = numZeros + 1;
@@ -3100,7 +3100,7 @@ StrictInt64Conversion(Double* dPtr,
  *
  *	Test whether bankers' rounding should round a digit up. Assumption
  *	is made that the denominator of the fraction being tested is
- *	a power of 2**DIGIT_BIT.
+ *	a power of 2**MP_DIGIT_BIT.
  *
  * Results:
  *	Returns 1 iff the fraction is more than 1/2, or if the fraction
@@ -3112,7 +3112,7 @@ StrictInt64Conversion(Double* dPtr,
 inline static int
 ShouldBankerRoundUpPowD(mp_int* b,
 				/* Numerator of the fraction */
-			int sd,	/* Denominator is 2**(sd*DIGIT_BIT) */
+			int sd,	/* Denominator is 2**(sd*MP_DIGIT_BIT) */
 			int isodd)
 				/* 1 if the digit is odd, 0 if even */
 {
@@ -3153,7 +3153,7 @@ ShouldBankerRoundUpToNextPowD(mp_int* b,
 			      mp_int* m,
 				/* Numerator of the rounding tolerance */
 			      int sd,
-				/* Common denominator is 2**(sd*DIGIT_BIT) */
+				/* Common denominator is 2**(sd*MP_DIGIT_BIT) */
 			      int convType,
 				/* Conversion type: STEELE defeats
 				 * round-to-even (Not sure why one wants to
@@ -3168,7 +3168,7 @@ ShouldBankerRoundUpToNextPowD(mp_int* b,
     /*
      * Compare B and S-m -- which is the same as comparing B+m and S --
      * which we do by computing b+m and doing a bitwhack compare against
-     * 2**(DIGIT_BIT*sd)
+     * 2**(MP_DIGIT_BIT*sd)
      */
     mp_add(b, m, temp);
     if (temp->used <= sd) {	/* too few digits to be > S */
@@ -3200,7 +3200,7 @@ ShouldBankerRoundUpToNextPowD(mp_int* b,
  *	digits that reconverts exactly to the given number, or to
  *	'ilim' digits if that will yield a shorter result. The denominator
  *	in David Gay's conversion algorithm is known to be a power of
- *	2**DIGIT_BIT, and hence the division in the main loop may be replaced
+ *	2**MP_DIGIT_BIT, and hence the division in the main loop may be replaced
  *	by a digit shift and mask.
  *
  * Results:
@@ -3289,7 +3289,7 @@ ShorteningBignumConversionPowD(Double* dPtr,
     }
     mp_init(&temp);
 
-    /* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
+    /* Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
      * by mp_digit extraction */
 
     i = 0;
@@ -3396,7 +3396,7 @@ ShorteningBignumConversionPowD(Double* dPtr,
  *	Converts a double-precision number to a fixed-lengt string of
  *	'ilim' digits (or 'ilim1' if log10(d) has been overestimated.)
  *	The denominator in David Gay's conversion algorithm is known to
- *	be a power of 2**DIGIT_BIT, and hence the division in the main
+ *	be a power of 2**MP_DIGIT_BIT, and hence the division in the main
  *	loop may be replaced by a digit shift and mask.
  *
  * Results:
@@ -3465,7 +3465,7 @@ StrictBignumConversionPowD(Double* dPtr,
     mp_init(&temp);
 
     /*
-     * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
+     * Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
      * by mp_digit extraction
      */
 
@@ -4234,7 +4234,7 @@ TclDoubleDigits(double dv,	/* Number to convert */
 	    /*
 	     * The denominator is a power of 2, so we can replace division
 	     * by digit shifts. First we round up s2 to a multiple of
-	     * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
+	     * MP_DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
 	     * into a version of the comparison that's specialized for
 	     * the 'power of mp_digit in the denominator' case.
 	     */
@@ -4294,7 +4294,7 @@ TclDoubleDigits(double dv,	/* Number to convert */
 	    /*
 	     * The denominator is a power of 2, so we can replace division
 	     * by digit shifts. First we round up s2 to a multiple of
-	     * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
+	     * MP_DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
 	     * into a version of the comparison that's specialized for
 	     * the 'power of mp_digit in the denominator' case.
 	     */
@@ -4573,10 +4573,10 @@ TclBignumToDouble(
     bits = mp_count_bits(a);
     if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
 	errno = ERANGE;
-	if (a->sign == MP_ZPOS) {
-	    return HUGE_VAL;
-	} else {
+	if (mp_isneg(a)) {
 	    return -HUGE_VAL;
+	} else {
+	    return HUGE_VAL;
 	}
     }
     shift = mantBits - bits;
@@ -4606,10 +4606,10 @@ TclBignumToDouble(
 
 	    mp_div_2d(a, -shift, &b, NULL);
 	    if (mp_isodd(&b)) {
-		if (b.sign == MP_ZPOS) {
-		    mp_add_d(&b, 1, &b);
-		} else {
+		if (mp_isneg(&b)) {
 		    mp_sub_d(&b, 1, &b);
+		} else {
+		    mp_add_d(&b, 1, &b);
 		}
 	    }
 	} else {
@@ -4619,10 +4619,10 @@ TclBignumToDouble(
 	     */
 
 	    mp_div_2d(a, -1-shift, &b, NULL);
-	    if (b.sign == MP_ZPOS) {
-		mp_add_d(&b, 1, &b);
-	    } else {
+	    if (mp_isneg(&b)) {
 		mp_sub_d(&b, 1, &b);
+	    } else {
+		mp_add_d(&b, 1, &b);
 	    }
 	    mp_div_2d(&b, 1, &b, NULL);
 	}
@@ -4648,10 +4648,10 @@ TclBignumToDouble(
      * Return the result with the appropriate sign.
      */
 
-    if (a->sign == MP_ZPOS) {
-	return r;
-    } else {
+    if (mp_isneg(a)) {
 	return -r;
+    } else {
+	return r;
     }
 }
 
@@ -4824,7 +4824,7 @@ BignumToBiasedFrExp(
      */
 
     *machexp = bits - mantBits + 2;
-    return ((a->sign == MP_ZPOS) ? r : -r);
+    return (mp_isneg(a) ? -r : r);
 }
 
 /*
diff --git a/generic/tclStringObj.c b/generic/tclStringObj.c
index 699dc5a..aeb4285 100644
--- a/generic/tclStringObj.c
+++ b/generic/tclStringObj.c
@@ -2093,7 +2093,7 @@ Tcl_AppendFormatToObj(
 		    if (Tcl_GetBignumFromObj(interp,segment,&big) != TCL_OK) {
 			goto error;
 		    }
-		    mp_mod_2d(&big, (int) CHAR_BIT*sizeof(Tcl_WideInt), &big);
+		    mp_mod_2d(&big, CHAR_BIT*sizeof(Tcl_WideInt), &big);
 		    objPtr = Tcl_NewBignumObj(&big);
 		    Tcl_IncrRefCount(objPtr);
 		    Tcl_GetWideIntFromObj(NULL, objPtr, &w);
@@ -2107,7 +2107,7 @@ Tcl_AppendFormatToObj(
 		    if (Tcl_GetBignumFromObj(interp,segment,&big) != TCL_OK) {
 			goto error;
 		    }
-		    mp_mod_2d(&big, (int) CHAR_BIT * sizeof(long), &big);
+		    mp_mod_2d(&big, CHAR_BIT * sizeof(long), &big);
 		    objPtr = Tcl_NewBignumObj(&big);
 		    Tcl_IncrRefCount(objPtr);
 		    TclGetLongFromObj(NULL, objPtr, &l);
diff --git a/generic/tclTestObj.c b/generic/tclTestObj.c
index 4226d51..b5d0d6b 100644
--- a/generic/tclTestObj.c
+++ b/generic/tclTestObj.c
@@ -269,9 +269,9 @@ TestbignumobjCmd(
 	    return TCL_ERROR;
 	}
 	if (!Tcl_IsShared(varPtr[varIndex])) {
-	    Tcl_SetIntObj(varPtr[varIndex], mp_iseven(&bignumValue));
+	    Tcl_SetIntObj(varPtr[varIndex], !mp_isodd(&bignumValue));
 	} else {
-	    SetVarToObj(varIndex, Tcl_NewIntObj(mp_iseven(&bignumValue)));
+	    SetVarToObj(varIndex, Tcl_NewIntObj(!mp_isodd(&bignumValue)));
 	}
 	mp_clear(&bignumValue);
 	break;
diff --git a/generic/tclTomMath.decls b/generic/tclTomMath.decls
index ab39e83..db37e41 100644
--- a/generic/tclTomMath.decls
+++ b/generic/tclTomMath.decls
@@ -1,9 +1,8 @@
 # tclTomMath.decls --
 #
-#	This file contains the declarations for the functions in
-#	'libtommath' that are contained within the Tcl library.
-#	This file is used to generate the 'tclTomMathDecls.h' and
-#	'tclStubInit.c' files.
+#	This file contains the declarations for the functions in 'libtommath'
+#	that are contained within the Tcl library.  This file is used to
+#	generate the 'tclTomMathDecls.h' and 'tclStubInit.c' files.
 #
 # If you edit this file, advance the revision number (and the epoch
 # if the new stubs are not backward compatible) in tclTomMathDecls.h
@@ -19,196 +18,197 @@ library tcl
 
 interface tclTomMath
 # hooks {tclTomMathInt}
+scspec EXTERN
 
 # Declare each of the functions in the Tcl tommath interface
 
-declare 0 generic {
+declare 0 {
     int TclBN_epoch(void)
 }
-declare 1 generic {
+declare 1 {
     int TclBN_revision(void)
 }
 
-declare 2 generic {
+declare 2 {
     int TclBN_mp_add(mp_int *a, mp_int *b, mp_int *c)
 }
-declare 3 generic {
+declare 3 {
     int TclBN_mp_add_d(mp_int *a, mp_digit b, mp_int *c)
 }
-declare 4 generic {
-    int TclBN_mp_and(mp_int *a, mp_int *b, mp_int *c)
+declare 4 {
+    int TclBN_mp_and(const mp_int *a, const mp_int *b, mp_int *c)
 }
-declare 5 generic {
+declare 5 {
     void TclBN_mp_clamp(mp_int *a)
 }
-declare 6 generic {
+declare 6 {
     void TclBN_mp_clear(mp_int *a)
 }
-declare 7 generic {
+declare 7 {
     void TclBN_mp_clear_multi(mp_int *a, ...)
 }
-declare 8 generic {
-    int TclBN_mp_cmp(mp_int *a, mp_int *b)
+declare 8 {
+    int TclBN_mp_cmp(const mp_int *a, const mp_int *b)
 }
-declare 9 generic {
-    int TclBN_mp_cmp_d(mp_int *a, mp_digit b)
+declare 9 {
+    int TclBN_mp_cmp_d(const mp_int *a, mp_digit b)
 }
-declare 10 generic {
-    int TclBN_mp_cmp_mag(mp_int *a, mp_int *b)
+declare 10 {
+    int TclBN_mp_cmp_mag(const mp_int *a, const mp_int *b)
 }
-declare 11 generic {
+declare 11 {
     int TclBN_mp_copy(mp_int *a, mp_int *b)
 }
-declare 12 generic {
+declare 12 {
     int TclBN_mp_count_bits(mp_int *a)
 }
-declare 13 generic {
+declare 13 {
     int TclBN_mp_div(mp_int *a, mp_int *b, mp_int *q, mp_int *r)
 }
-declare 14 generic {
+declare 14 {
     int TclBN_mp_div_d(mp_int *a, mp_digit b, mp_int *q, mp_digit *r)
 }
-declare 15 generic {
+declare 15 {
     int TclBN_mp_div_2(mp_int *a, mp_int *q)
 }
-declare 16 generic {
+declare 16 {
     int TclBN_mp_div_2d(mp_int *a, int b, mp_int *q, mp_int *r)
 }
-declare 17 generic {
+declare 17 {
     int TclBN_mp_div_3(mp_int *a, mp_int *q, mp_digit *r)
 }
-declare 18 generic {
+declare 18 {
     void TclBN_mp_exch(mp_int *a, mp_int *b)
 }
-declare 19 generic {
+declare 19 {
     int TclBN_mp_expt_d(mp_int *a, mp_digit b, mp_int *c)
 }
-declare 20 generic {
+declare 20 {
     int TclBN_mp_grow(mp_int *a, int size)
 }
-declare 21 generic {
+declare 21 {
     int TclBN_mp_init(mp_int *a)
 }
-declare 22 generic {
+declare 22 {
     int TclBN_mp_init_copy(mp_int *a, mp_int *b)
 }
-declare 23 generic {
+declare 23 {
     int TclBN_mp_init_multi(mp_int *a, ...)
 }
-declare 24 generic {
+declare 24 {
     int TclBN_mp_init_set(mp_int *a, mp_digit b)
 }
-declare 25 generic {
+declare 25 {
     int TclBN_mp_init_size(mp_int *a, int size)
 }
-declare 26 generic {
+declare 26 {
     int TclBN_mp_lshd(mp_int *a, int shift)
 }
-declare 27 generic {
+declare 27 {
     int TclBN_mp_mod(mp_int *a, mp_int *b, mp_int *r)
 }
-declare 28 generic {
+declare 28 {
     int TclBN_mp_mod_2d(mp_int *a, int b, mp_int *r)
 }
-declare 29 generic {
+declare 29 {
     int TclBN_mp_mul(mp_int *a, mp_int *b, mp_int *p)
 }
-declare 30 generic {
+declare 30 {
     int TclBN_mp_mul_d(mp_int *a, mp_digit b, mp_int *p)
 }
-declare 31 generic {
+declare 31 {
     int TclBN_mp_mul_2(mp_int *a, mp_int *p)
 }
-declare 32 generic {
+declare 32 {
     int TclBN_mp_mul_2d(mp_int *a, int d, mp_int *p)
 }
-declare 33 generic {
+declare 33 {
     int TclBN_mp_neg(mp_int *a, mp_int *b)
 }
-declare 34 generic {
-    int TclBN_mp_or(mp_int *a, mp_int *b, mp_int *c)
+declare 34 {
+    int TclBN_mp_or(const mp_int *a, const mp_int *b, mp_int *c)
 }
-declare 35 generic {
+declare 35 {
     int TclBN_mp_radix_size(mp_int *a, int radix, int *size)
 }
-declare 36 generic {
+declare 36 {
     int TclBN_mp_read_radix(mp_int *a, const char *str, int radix)
 }
-declare 37 generic {
+declare 37 {
     void TclBN_mp_rshd(mp_int *a, int shift)
 }
-declare 38 generic {
+declare 38 {
     int TclBN_mp_shrink(mp_int *a)
 }
-declare 39 generic {
+declare 39 {
     void TclBN_mp_set(mp_int *a, mp_digit b)
 }
-declare 40 generic {
+declare 40 {
     int TclBN_mp_sqr(mp_int *a, mp_int *b)
 }
-declare 41 generic {
+declare 41 {
     int TclBN_mp_sqrt(mp_int *a, mp_int *b)
 }
-declare 42 generic {
+declare 42 {
     int TclBN_mp_sub(mp_int *a, mp_int *b, mp_int *c)
 }
-declare 43 generic {
+declare 43 {
     int TclBN_mp_sub_d(mp_int *a, mp_digit b, mp_int *c)
 }
-declare 44 generic {
+declare 44 {
     int TclBN_mp_to_unsigned_bin(mp_int *a, unsigned char *b)
 }
-declare 45 generic {
+declare 45 {
     int TclBN_mp_to_unsigned_bin_n(mp_int *a, unsigned char *b,
 	    unsigned long *outlen)
 }
-declare 46 generic {
+declare 46 {
     int TclBN_mp_toradix_n(mp_int *a, char *str, int radix, int maxlen)
 }
-declare 47 generic {
+declare 47 {
     int TclBN_mp_unsigned_bin_size(mp_int *a)
 }
-declare 48 generic {
-    int TclBN_mp_xor(mp_int *a, mp_int *b, mp_int *c)
+declare 48 {
+    int TclBN_mp_xor(const mp_int *a, const mp_int *b, mp_int *c)
 }
-declare 49 generic {
+declare 49 {
     void TclBN_mp_zero(mp_int *a)
 }
 
 # internal routines to libtommath - should not be called but must be
 # exported to accommodate the "tommath" extension
 
-declare 50 generic {
+declare 50 {
     void TclBN_reverse(unsigned char *s, int len)
 }
-declare 51 generic {
+declare 51 {
     int TclBN_fast_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
 }
-declare 52 generic {
+declare 52 {
     int TclBN_fast_s_mp_sqr(mp_int *a, mp_int *b)
 }
-declare 53 generic {
+declare 53 {
     int TclBN_mp_karatsuba_mul(mp_int *a, mp_int *b, mp_int *c)
 }
-declare 54 generic {
+declare 54 {
     int TclBN_mp_karatsuba_sqr(mp_int *a, mp_int *b)
 }
-declare 55 generic {
+declare 55 {
     int TclBN_mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
 }
-declare 56 generic {
+declare 56 {
     int TclBN_mp_toom_sqr(mp_int *a, mp_int *b)
 }
-declare 57 generic {
+declare 57 {
     int TclBN_s_mp_add(mp_int *a, mp_int *b, mp_int *c)
 }
-declare 58 generic {
+declare 58 {
     int TclBN_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
 }
-declare 59 generic {
+declare 59 {
     int TclBN_s_mp_sqr(mp_int *a, mp_int *b)
 }
-declare 60 generic {
+declare 60 {
     int TclBN_s_mp_sub(mp_int *a, mp_int *b, mp_int *c)
 }
 declare 61 {
diff --git a/generic/tclTomMath.h b/generic/tclTomMath.h
index b435d57..b219405 100644
--- a/generic/tclTomMath.h
+++ b/generic/tclTomMath.h
@@ -241,13 +241,13 @@ typedef int ltm_prime_callback(unsigned char *dst, int len, void *dat);
 
 /* 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 */
@@ -257,7 +257,7 @@ 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 */
@@ -272,23 +272,24 @@ 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 */
 /*
@@ -676,14 +677,14 @@ int mp_prime_is_divisible(mp_int *a, int *result);
  * 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
@@ -701,7 +702,7 @@ int mp_prime_rabin_miller_trials(int size);
  * 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
@@ -710,7 +711,7 @@ int mp_prime_is_prime(mp_int *a, int t, int *result);
  * 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),
@@ -728,9 +729,9 @@ int mp_prime_next_prime(mp_int *a, int t, int bbs_style);
  *
  * 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
diff --git a/generic/tclTomMathDecls.h b/generic/tclTomMathDecls.h
index 7113f69..bd801a3 100644
--- a/generic/tclTomMathDecls.h
+++ b/generic/tclTomMathDecls.h
@@ -44,13 +44,6 @@
 #define XREALLOC(x,n) TclBNRealloc(x,n)
 #define XCALLOC(n,x) TclBNCalloc(n,x)
 
-/* Rename the global symbols in libtommath to avoid linkage conflicts */
-
-#define KARATSUBA_MUL_CUTOFF TclBNKaratsubaMulCutoff
-#define KARATSUBA_SQR_CUTOFF TclBNKaratsubaSqrCutoff
-#define TOOM_MUL_CUTOFF TclBNToomMulCutoff
-#define TOOM_SQR_CUTOFF TclBNToomSqrCutoff
-
 #define bn_reverse TclBN_reverse
 #define s_mp_reverse TclBN_reverse
 #define fast_s_mp_mul_digs TclBN_fast_s_mp_mul_digs
@@ -172,7 +165,8 @@ EXTERN int		TclBN_mp_add_d(mp_int *a, mp_digit b, mp_int *c);
 #ifndef TclBN_mp_and_TCL_DECLARED
 #define TclBN_mp_and_TCL_DECLARED
 /* 4 */
-EXTERN int		TclBN_mp_and(mp_int *a, mp_int *b, mp_int *c);
+EXTERN int		TclBN_mp_and(CONST mp_int *a, CONST mp_int *b,
+				mp_int *c);
 #endif
 #ifndef TclBN_mp_clamp_TCL_DECLARED
 #define TclBN_mp_clamp_TCL_DECLARED
@@ -192,17 +186,17 @@ EXTERN void		TclBN_mp_clear_multi(mp_int *a, ...);
 #ifndef TclBN_mp_cmp_TCL_DECLARED
 #define TclBN_mp_cmp_TCL_DECLARED
 /* 8 */
-EXTERN int		TclBN_mp_cmp(mp_int *a, mp_int *b);
+EXTERN int		TclBN_mp_cmp(CONST mp_int *a, CONST mp_int *b);
 #endif
 #ifndef TclBN_mp_cmp_d_TCL_DECLARED
 #define TclBN_mp_cmp_d_TCL_DECLARED
 /* 9 */
-EXTERN int		TclBN_mp_cmp_d(mp_int *a, mp_digit b);
+EXTERN int		TclBN_mp_cmp_d(CONST mp_int *a, mp_digit b);
 #endif
 #ifndef TclBN_mp_cmp_mag_TCL_DECLARED
 #define TclBN_mp_cmp_mag_TCL_DECLARED
 /* 10 */
-EXTERN int		TclBN_mp_cmp_mag(mp_int *a, mp_int *b);
+EXTERN int		TclBN_mp_cmp_mag(CONST mp_int *a, CONST mp_int *b);
 #endif
 #ifndef TclBN_mp_copy_TCL_DECLARED
 #define TclBN_mp_copy_TCL_DECLARED
@@ -325,7 +319,8 @@ EXTERN int		TclBN_mp_neg(mp_int *a, mp_int *b);
 #ifndef TclBN_mp_or_TCL_DECLARED
 #define TclBN_mp_or_TCL_DECLARED
 /* 34 */
-EXTERN int		TclBN_mp_or(mp_int *a, mp_int *b, mp_int *c);
+EXTERN int		TclBN_mp_or(CONST mp_int *a, CONST mp_int *b,
+				mp_int *c);
 #endif
 #ifndef TclBN_mp_radix_size_TCL_DECLARED
 #define TclBN_mp_radix_size_TCL_DECLARED
@@ -398,7 +393,8 @@ EXTERN int		TclBN_mp_unsigned_bin_size(mp_int *a);
 #ifndef TclBN_mp_xor_TCL_DECLARED
 #define TclBN_mp_xor_TCL_DECLARED
 /* 48 */
-EXTERN int		TclBN_mp_xor(mp_int *a, mp_int *b, mp_int *c);
+EXTERN int		TclBN_mp_xor(CONST mp_int *a, CONST mp_int *b,
+				mp_int *c);
 #endif
 #ifndef TclBN_mp_zero_TCL_DECLARED
 #define TclBN_mp_zero_TCL_DECLARED
@@ -487,13 +483,13 @@ typedef struct TclTomMathStubs {
     int (*tclBN_revision) (void); /* 1 */
     int (*tclBN_mp_add) (mp_int *a, mp_int *b, mp_int *c); /* 2 */
     int (*tclBN_mp_add_d) (mp_int *a, mp_digit b, mp_int *c); /* 3 */
-    int (*tclBN_mp_and) (mp_int *a, mp_int *b, mp_int *c); /* 4 */
+    int (*tclBN_mp_and) (CONST mp_int *a, CONST mp_int *b, mp_int *c); /* 4 */
     void (*tclBN_mp_clamp) (mp_int *a); /* 5 */
     void (*tclBN_mp_clear) (mp_int *a); /* 6 */
     void (*tclBN_mp_clear_multi) (mp_int *a, ...); /* 7 */
-    int (*tclBN_mp_cmp) (mp_int *a, mp_int *b); /* 8 */
-    int (*tclBN_mp_cmp_d) (mp_int *a, mp_digit b); /* 9 */
-    int (*tclBN_mp_cmp_mag) (mp_int *a, mp_int *b); /* 10 */
+    int (*tclBN_mp_cmp) (CONST mp_int *a, CONST mp_int *b); /* 8 */
+    int (*tclBN_mp_cmp_d) (CONST mp_int *a, mp_digit b); /* 9 */
+    int (*tclBN_mp_cmp_mag) (CONST mp_int *a, CONST mp_int *b); /* 10 */
     int (*tclBN_mp_copy) (mp_int *a, mp_int *b); /* 11 */
     int (*tclBN_mp_count_bits) (mp_int *a); /* 12 */
     int (*tclBN_mp_div) (mp_int *a, mp_int *b, mp_int *q, mp_int *r); /* 13 */
@@ -517,7 +513,7 @@ typedef struct TclTomMathStubs {
     int (*tclBN_mp_mul_2) (mp_int *a, mp_int *p); /* 31 */
     int (*tclBN_mp_mul_2d) (mp_int *a, int d, mp_int *p); /* 32 */
     int (*tclBN_mp_neg) (mp_int *a, mp_int *b); /* 33 */
-    int (*tclBN_mp_or) (mp_int *a, mp_int *b, mp_int *c); /* 34 */
+    int (*tclBN_mp_or) (CONST mp_int *a, CONST mp_int *b, mp_int *c); /* 34 */
     int (*tclBN_mp_radix_size) (mp_int *a, int radix, int *size); /* 35 */
     int (*tclBN_mp_read_radix) (mp_int *a, CONST char *str, int radix); /* 36 */
     void (*tclBN_mp_rshd) (mp_int *a, int shift); /* 37 */
@@ -531,7 +527,7 @@ typedef struct TclTomMathStubs {
     int (*tclBN_mp_to_unsigned_bin_n) (mp_int *a, unsigned char *b, unsigned long *outlen); /* 45 */
     int (*tclBN_mp_toradix_n) (mp_int *a, char *str, int radix, int maxlen); /* 46 */
     int (*tclBN_mp_unsigned_bin_size) (mp_int *a); /* 47 */
-    int (*tclBN_mp_xor) (mp_int *a, mp_int *b, mp_int *c); /* 48 */
+    int (*tclBN_mp_xor) (CONST mp_int *a, CONST mp_int *b, mp_int *c); /* 48 */
     void (*tclBN_mp_zero) (mp_int *a); /* 49 */
     void (*tclBN_reverse) (unsigned char *s, int len); /* 50 */
     int (*tclBN_fast_s_mp_mul_digs) (mp_int *a, mp_int *b, mp_int *c, int digs); /* 51 */
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)