Change the compiler's big integer library to use libTomMath

This now replaces Bill's crappy big int implementation

Makefile

@@ -41,13 +41,16 @@ demo:
 	./odin run examples/demo/demo.odin
 
 debug:
-	$(CC) src/main.cpp $(DISABLED_WARNINGS) $(CFLAGS) -g $(LDFLAGS) -o odin
+	$(CC) src/main.cpp src/libtommath.c $(DISABLED_WARNINGS) $(CFLAGS) -g $(LDFLAGS) -o odin
 
 release:
-	$(CC) src/main.cpp $(DISABLED_WARNINGS) $(CFLAGS) -O3 -march=native $(LDFLAGS) -o odin
+	$(CC) src/main.cpp src/libtommath.c $(DISABLED_WARNINGS) $(CFLAGS) -O3 $(LDFLAGS) -o odin
+
+release-native:
+    $(CC) src/main.cpp src/libtommath.c $(DISABLED_WARNINGS) $(CFLAGS) -O3 -march=native $(LDFLAGS) -o odin
 
 nightly:
-	$(CC) src/main.cpp $(DISABLED_WARNINGS) $(CFLAGS) -DNIGHTLY -O3 $(LDFLAGS) -o odin
+	$(CC) src/main.cpp src/libtommath.c $(DISABLED_WARNINGS) $(CFLAGS) -DNIGHTLY -O3 $(LDFLAGS) -o odin
 
 
 

build.bat

@@ -46,7 +46,7 @@ if %release_mode% EQU 0 ( rem Debug
 
 set compiler_warnings= ^
 	-W4 -WX ^
-	-wd4100 -wd4101 -wd4127 -wd4189 ^
+	-wd4100 -wd4101 -wd4127 -wd4146 -wd4189 ^
 	-wd4201 -wd4204 ^
 	-wd4456 -wd4457 -wd4480 ^
 	-wd4512
@@ -70,10 +70,19 @@ set linker_settings=%libs% %linker_flags%
 del *.pdb > NUL 2> NUL
 del *.ilk > NUL 2> NUL
 
-cl %compiler_settings% "src\main.cpp" /link %linker_settings% -OUT:%exe_name%
+cl %compiler_settings% "src\main.cpp" "src\libtommath.c" /link %linker_settings% -OUT:%exe_name%
+rem cl %compiler_settings% "src\main.cpp" /link %linker_settings% -OUT:%exe_name%
 
 if %errorlevel% neq 0 goto end_of_build
-if %release_mode% EQU 0 odin run examples/demo/demo.odin
+
+odin run examples/demo
+rem odin run examples/bug -show-more-timings -threaded-checker
+rem odin check examples/demo -threaded-checker
+rem odin build examples/demo
+rem odin run examples/demo -threaded-checker
+rem odin run examples/bug -show-more-timings
+rem if %release_mode% EQU 0 odin check examples/bug -show-more-timings
+rem if %release_mode% EQU 0 odin check examples/bug
 
 del *.obj > NUL 2> NUL
 

src/check_builtin.cpp

@@ -745,7 +745,7 @@ bool check_builtin_procedure(CheckerContext *c, Operand *operand, Ast *call, i32
 				return false;
 			}
 
-			if (op.value.value_integer.neg) {
+			if (op.value.value_integer.sign) {
 				error(op.expr, "Negative 'swizzle' index");
 				return false;
 			}
@@ -795,10 +795,12 @@ bool check_builtin_procedure(CheckerContext *c, Operand *operand, Ast *call, i32
 		convert_to_typed(c, &y, x.type); if (y.mode == Addressing_Invalid) return false;
 		if (x.mode == Addressing_Constant &&
 		    y.mode == Addressing_Constant) {
-			if (is_type_numeric(x.type) && exact_value_imag(x.value).value_float == 0) {
+			x.value = exact_value_to_float(x.value);
+		    	y.value = exact_value_to_float(y.value);
+			if (is_type_numeric(x.type) && x.value.kind == ExactValue_Float) {
 				x.type = t_untyped_float;
 			}
-			if (is_type_numeric(y.type) && exact_value_imag(y.value).value_float == 0) {
+			if (is_type_numeric(y.type) && y.value.kind == ExactValue_Float) {
 				y.type = t_untyped_float;
 			}
 		}
@@ -882,16 +884,20 @@ bool check_builtin_procedure(CheckerContext *c, Operand *operand, Ast *call, i32
 		    y.mode == Addressing_Constant &&
 		    z.mode == Addressing_Constant &&
 		    w.mode == Addressing_Constant) {
-			if (is_type_numeric(x.type) && exact_value_imag(x.value).value_float == 0) {
+		    	x.value = exact_value_to_float(x.value);
+		    	y.value = exact_value_to_float(y.value);
+		    	z.value = exact_value_to_float(z.value);
+		    	w.value = exact_value_to_float(w.value);
+			if (is_type_numeric(x.type) && x.value.kind == ExactValue_Float) {
 				x.type = t_untyped_float;
 			}
-			if (is_type_numeric(y.type) && exact_value_imag(y.value).value_float == 0) {
+			if (is_type_numeric(y.type) && y.value.kind == ExactValue_Float) {
 				y.type = t_untyped_float;
 			}
-			if (is_type_numeric(z.type) && exact_value_imag(z.value).value_float == 0) {
+			if (is_type_numeric(z.type) && z.value.kind == ExactValue_Float) {
 				z.type = t_untyped_float;
 			}
-			if (is_type_numeric(w.type) && exact_value_imag(w.value).value_float == 0) {
+			if (is_type_numeric(w.type) && w.value.kind == ExactValue_Float) {
 				w.type = t_untyped_float;
 			}
 		}
@@ -1484,7 +1490,7 @@ bool check_builtin_procedure(CheckerContext *c, Operand *operand, Ast *call, i32
 		if (operand->mode == Addressing_Constant) {
 			switch (operand->value.kind) {
 			case ExactValue_Integer:
-				operand->value.value_integer.neg = false;
+				mp_abs(&operand->value.value_integer, &operand->value.value_integer);
 				break;
 			case ExactValue_Float:
 				operand->value.value_float = gb_abs(operand->value.value_float);
@@ -1837,7 +1843,7 @@ bool check_builtin_procedure(CheckerContext *c, Operand *operand, Ast *call, i32
 			operand->type = t_invalid;
 			return false;
 		}
-		if (x.value.value_integer.neg) {
+		if (x.value.value_integer.sign) {
 			error(call, "Negative vector element length");
 			operand->mode = Addressing_Type;
 			operand->type = t_invalid;
@@ -1877,7 +1883,7 @@ bool check_builtin_procedure(CheckerContext *c, Operand *operand, Ast *call, i32
 			operand->type = t_invalid;
 			return false;
 		}
-		if (x.value.value_integer.neg) {
+		if (x.value.value_integer.sign) {
 			error(call, "Negative array element length");
 			operand->mode = Addressing_Type;
 			operand->type = t_invalid;

src/check_expr.cpp

@@ -1503,13 +1503,13 @@ bool check_representable_as_constant(CheckerContext *c, ExactValue in_value, Typ
 			big_int_from_i64(&bi127, 127);
 
 			big_int_shl_eq(&umax, &bi128);
-			big_int_sub_eq(&umax, &BIG_INT_ONE);
+			mp_decr(&umax);
 
 			big_int_shl_eq(&imin, &bi127);
 			big_int_neg(&imin, &imin);
 
 			big_int_shl_eq(&imax, &bi127);
-			big_int_sub_eq(&imax, &BIG_INT_ONE);
+			mp_decr(&imax);
 		}
 
 		switch (type->Basic.kind) {
@@ -1555,7 +1555,7 @@ bool check_representable_as_constant(CheckerContext *c, ExactValue in_value, Typ
 			{
 				// return 0ull <= i && i <= umax;
 				int b = big_int_cmp(&i, &umax);
-				return !i.neg && (b <= 0);
+				return !i.sign && (b <= 0);
 			}
 
 		case Basic_UntypedInteger:
@@ -1758,12 +1758,6 @@ void check_is_expressible(CheckerContext *ctx, Operand *o, Type *type) {
 			if (!is_type_integer(o->type) && is_type_integer(type)) {
 				error(o->expr, "'%s' truncated to '%s'", a, b);
 			} else {
-				#if 0
-				gb_printf_err("AddressingMode, %d\n", o->mode);
-				gb_printf_err("ExactValueKind, %d\n", o->value.kind);
-				bool ok = check_representable_as_constant(ctx, o->value, type, &out_value);
-				gb_printf_err("ok, %d\n", ok);
-				#endif
 				error(o->expr, "Cannot convert numeric value '%s' to '%s' from '%s", a, b, c);
 				check_assignment_error_suggestion(ctx, o, type);
 			}
@@ -2206,7 +2200,7 @@ void check_shift(CheckerContext *c, Operand *x, Operand *y, Ast *node, Type *typ
 			}
 
 			BigInt max_shift = {};
-			big_int_from_u64(&max_shift, 128);
+			big_int_from_u64(&max_shift, MAX_BIG_INT_SHIFT);
 
 			if (big_int_cmp(&y_val.value_integer, &max_shift) > 0) {
 				gbString err_str = expr_to_string(y->expr);
@@ -2248,7 +2242,7 @@ void check_shift(CheckerContext *c, Operand *x, Operand *y, Ast *node, Type *typ
 		}
 	}
 
-	if (y->mode == Addressing_Constant && y->value.value_integer.neg) {
+	if (y->mode == Addressing_Constant && y->value.value_integer.sign) {
 		gbString err_str = expr_to_string(y->expr);
 		error(node, "Shift amount cannot be negative: '%s'", err_str);
 		gb_string_free(err_str);
@@ -3320,7 +3314,7 @@ bool check_index_value(CheckerContext *c, bool open_range, Ast *index_value, i64
 	if (operand.mode == Addressing_Constant &&
 	    (c->state_flags & StateFlag_no_bounds_check) == 0) {
 		BigInt i = exact_value_to_integer(operand.value).value_integer;
-		if (i.neg && !is_type_enum(index_type)) {
+		if (i.sign && !is_type_enum(index_type)) {
 			gbString expr_str = expr_to_string(operand.expr);
 			error(operand.expr, "Index '%s' cannot be a negative value", expr_str);
 			gb_string_free(expr_str);
@@ -3366,7 +3360,7 @@ bool check_index_value(CheckerContext *c, bool open_range, Ast *index_value, i64
 
 			} else { // NOTE(bill): Do array bound checking
 				i64 v = -1;
-				if (i.len <= 1) {
+				if (i.used <= 1) {
 					v = big_int_to_i64(&i);
 				}
 				if (value) *value = v;

src/check_type.cpp

@@ -207,7 +207,7 @@ bool check_custom_align(CheckerContext *ctx, Ast *node, i64 *align_) {
 	if (is_type_untyped(type) || is_type_integer(type)) {
 		if (o.value.kind == ExactValue_Integer) {
 			BigInt v = o.value.value_integer;
-			if (v.len > 1) {
+			if (v.used > 1) {
 				gbAllocator a = heap_allocator();
 				String str = big_int_to_string(a, &v);
 				error(node, "#align too large, %.*s", LIT(str));
@@ -1998,16 +1998,16 @@ i64 check_array_count(CheckerContext *ctx, Operand *o, Ast *e) {
 	if (is_type_untyped(type) || is_type_integer(type)) {
 		if (o->value.kind == ExactValue_Integer) {
 			BigInt count = o->value.value_integer;
-			if (o->value.value_integer.neg) {
+			if (o->value.value_integer.sign) {
 				gbAllocator a = heap_allocator();
 				String str = big_int_to_string(a, &count);
 				error(e, "Invalid negative array count, %.*s", LIT(str));
 				gb_free(a, str.text);
 				return 0;
 			}
-			switch (count.len) {
+			switch (count.used) {
 			case 0: return 0;
-			case 1: return count.d.word;
+			case 1: return big_int_to_u64(&count);
 			}
 			gbAllocator a = heap_allocator();
 			String str = big_int_to_string(a, &count);

src/exact_value.cpp

@@ -75,8 +75,8 @@ HashKey hash_exact_value(ExactValue v) {
 		}
 	case ExactValue_Integer:
 		{
-			HashKey key = hashing_proc(big_int_ptr(&v.value_integer), v.value_integer.len * gb_size_of(u64));
-			u8 last = (u8)v.value_integer.neg;
+			HashKey key = hashing_proc(v.value_integer.dp, gb_size_of(*v.value_integer.dp) * v.value_integer.used);
+			u8 last = (u8)v.value_integer.sign;
 			key.key = (key.key ^ last) * 0x100000001b3ll;
 			return key;
 		}
@@ -719,7 +719,6 @@ ExactValue exact_binary_operator_value(TokenKind op, ExactValue x, ExactValue y)
 		case Token_Shr:    big_int_shr(&c, a, b);     break;
 		default: goto error;
 		}
-		big_int_normalize(&c);
 		ExactValue res = {ExactValue_Integer};
 		res.value_integer = c;
 		return res;

src/libtommath/LICENSE

@@ -0,0 +1,26 @@
+                          The LibTom license
+
+This is free and unencumbered software released into the public domain.
+
+Anyone is free to copy, modify, publish, use, compile, sell, or
+distribute this software, either in source code form or as a compiled
+binary, for any purpose, commercial or non-commercial, and by any
+means.
+
+In jurisdictions that recognize copyright laws, the author or authors
+of this software dedicate any and all copyright interest in the
+software to the public domain. We make this dedication for the benefit
+of the public at large and to the detriment of our heirs and
+successors. We intend this dedication to be an overt act of
+relinquishment in perpetuity of all present and future rights to this
+software under copyright law.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
+OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
+ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+OTHER DEALINGS IN THE SOFTWARE.
+
+For more information, please refer to <http://unlicense.org/>

src/libtommath/README.md

@@ -0,0 +1,44 @@
+# libtommath
+
+This is the git repository for [LibTomMath](http://www.libtom.net/LibTomMath/), a free open source portable number theoretic multiple-precision integer (MPI) library written entirely in C.
+
+## Build Status
+
+### Travis CI
+
+master: [![Build Status](https://api.travis-ci.org/libtom/libtommath.png?branch=master)](https://travis-ci.org/libtom/libtommath)
+
+develop: [![Build Status](https://api.travis-ci.org/libtom/libtommath.png?branch=develop)](https://travis-ci.org/libtom/libtommath)
+
+### AppVeyor
+
+master: [![Build status](https://ci.appveyor.com/api/projects/status/b80lpolw3i8m6hsh/branch/master?svg=true)](https://ci.appveyor.com/project/libtom/libtommath/branch/master)
+
+develop: [![Build status](https://ci.appveyor.com/api/projects/status/b80lpolw3i8m6hsh/branch/develop?svg=true)](https://ci.appveyor.com/project/libtom/libtommath/branch/develop)
+
+### ABI Laboratory
+
+API/ABI changes: [check here](https://abi-laboratory.pro/tracker/timeline/libtommath/)
+
+## Summary
+
+The `develop` branch contains the in-development version. Stable releases are tagged.
+
+Documentation is built from the LaTeX file `bn.tex`. There is also limited documentation in `tommath.h`.
+There is also a document, `tommath.pdf`, which describes the goals of the project and many of the algorithms used.
+
+The project can be build by using `make`. Along with the usual `make`, `make clean` and `make install`,
+there are several other build targets, see the makefile for details.
+There are also makefiles for certain specific platforms.
+
+## Testing
+
+Tests are located in `demo/` and can be built in two flavors.
+* `make test` creates a stand-alone test binary that executes several test routines.
+* `make mtest_opponent` creates a test binary that is intended to be run against `mtest`.
+  `mtest` can be built with `make mtest` and test execution is done like `./mtest/mtest | ./mtest_opponent`.
+  `mtest` is creating test vectors using an alternative MPI library and `test` is consuming these vectors to verify correct behavior of ltm
+
+## Building and Installing
+
+Building is straightforward for GNU Linux only, the section "Building LibTomMath" in the documentation in `doc/bn.pdf` has the details.

src/libtommath/mp_2expt.c

@@ -0,0 +1,31 @@
+#include "tommath_private.h"
+#ifdef MP_2EXPT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* computes a = 2**b
+ *
+ * Simple algorithm which zeroes the int, grows it then just sets one bit
+ * as required.
+ */
+mp_err mp_2expt(mp_int *a, int b)
+{
+   mp_err    err;
+
+   /* zero a as per default */
+   mp_zero(a);
+
+   /* grow a to accomodate the single bit */
+   if ((err = mp_grow(a, (b / MP_DIGIT_BIT) + 1)) != MP_OKAY) {
+      return err;
+   }
+
+   /* set the used count of where the bit will go */
+   a->used = (b / MP_DIGIT_BIT) + 1;
+
+   /* put the single bit in its place */
+   a->dp[b / MP_DIGIT_BIT] = (mp_digit)1 << (mp_digit)(b % MP_DIGIT_BIT);
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_abs.c

@@ -0,0 +1,24 @@
+#include "tommath_private.h"
+#ifdef MP_ABS_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* b = |a|
+ *
+ * Simple function copies the input and fixes the sign to positive
+ */
+mp_err mp_abs(const mp_int *a, mp_int *b)
+{
+   mp_err err;
+
+   /* copy a to b */
+   if ((err = mp_copy(a, b)) != MP_OKAY) {
+      return err;
+   }
+
+   /* force the sign of b to positive */
+   b->sign = MP_ZPOS;
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_add.c

@@ -0,0 +1,29 @@
+#include "tommath_private.h"
+#ifdef MP_ADD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* high level addition (handles signs) */
+mp_err mp_add(const mp_int *a, const mp_int *b, mp_int *c)
+{
+   /* handle two cases, not four */
+   if (a->sign == b->sign) {
+      /* both positive or both negative */
+      /* add their magnitudes, copy the sign */
+      c->sign = a->sign;
+      return s_mp_add(a, b, c);
+   }
+
+   /* one positive, the other negative */
+   /* subtract the one with the greater magnitude from */
+   /* the one of the lesser magnitude. The result gets */
+   /* the sign of the one with the greater magnitude. */
+   if (mp_cmp_mag(a, b) == MP_LT) {
+      MP_EXCH(const mp_int *, a, b);
+   }
+
+   c->sign = a->sign;
+   return s_mp_sub(a, b, c);
+}
+
+#endif

src/libtommath/mp_add_d.c

@@ -0,0 +1,86 @@
+#include "tommath_private.h"
+#ifdef MP_ADD_D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* single digit addition */
+mp_err mp_add_d(const mp_int *a, mp_digit b, mp_int *c)
+{
+   mp_err err;
+   int oldused;
+
+   /* fast path for a == c */
+   if (a == c) {
+      if (!mp_isneg(c) &&
+          !mp_iszero(c) &&
+          ((c->dp[0] + b) < MP_DIGIT_MAX)) {
+         c->dp[0] += b;
+         return MP_OKAY;
+      }
+      if (mp_isneg(c) &&
+          (c->dp[0] > b)) {
+         c->dp[0] -= b;
+         return MP_OKAY;
+      }
+   }
+
+   /* grow c as required */
+   if ((err = mp_grow(c, a->used + 1)) != MP_OKAY) {
+      return err;
+   }
+
+   /* if a is negative and |a| >= b, call c = |a| - b */
+   if (mp_isneg(a) && ((a->used > 1) || (a->dp[0] >= b))) {
+      mp_int a_ = *a;
+      /* temporarily fix sign of a */
+      a_.sign = MP_ZPOS;
+
+      /* c = |a| - b */
+      err = mp_sub_d(&a_, b, c);
+
+      /* fix sign  */
+      c->sign = MP_NEG;
+
+      /* clamp */
+      mp_clamp(c);
+
+      return err;
+   }
+
+   /* old number of used digits in c */
+   oldused = c->used;
+
+   /* if a is positive */
+   if (!mp_isneg(a)) {
+      /* add digits, mu is carry */
+      int i;
+      mp_digit mu = b;
+      for (i = 0; i < a->used; i++) {
+         c->dp[i] = a->dp[i] + mu;
+         mu = c->dp[i] >> MP_DIGIT_BIT;
+         c->dp[i] &= MP_MASK;
+      }
+      /* set final carry */
+      c->dp[i] = mu;
+
+      /* setup size */
+      c->used = a->used + 1;
+   } else {
+      /* a was negative and |a| < b */
+      c->used = 1;
+
+      /* the result is a single digit */
+      c->dp[0] = (a->used == 1) ? b - a->dp[0] : b;
+   }
+
+   /* sign always positive */
+   c->sign = MP_ZPOS;
+
+   /* now zero to oldused */
+   s_mp_zero_digs(c->dp + c->used, oldused - c->used);
+   mp_clamp(c);
+
+   return MP_OKAY;
+}
+
+#endif

src/libtommath/mp_addmod.c

@@ -0,0 +1,15 @@
+#include "tommath_private.h"
+#ifdef MP_ADDMOD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* d = a + b (mod c) */
+mp_err mp_addmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d)
+{
+   mp_err err;
+   if ((err = mp_add(a, b, d)) != MP_OKAY) {
+      return err;
+   }
+   return mp_mod(d, c, d);
+}
+#endif

src/libtommath/mp_and.c

@@ -0,0 +1,54 @@
+#include "tommath_private.h"
+#ifdef MP_AND_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 used = MP_MAX(a->used, b->used) + 1, i;
+   mp_err err;
+   mp_digit ac = 1, bc = 1, cc = 1;
+   bool neg = (mp_isneg(a) && mp_isneg(b));
+
+   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 (mp_isneg(a)) {
+         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 (mp_isneg(b)) {
+         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 (neg) {
+         cc += ~c->dp[i] & MP_MASK;
+         c->dp[i] = cc & MP_MASK;
+         cc >>= MP_DIGIT_BIT;
+      }
+   }
+
+   c->used = used;
+   c->sign = (neg ? MP_NEG : MP_ZPOS);
+   mp_clamp(c);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_clamp.c

@@ -0,0 +1,27 @@
+#include "tommath_private.h"
+#ifdef MP_CLAMP_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* trim unused digits
+ *
+ * This is used to ensure that leading zero digits are
+ * trimed and the leading "used" digit will be non-zero
+ * Typically very fast.  Also fixes the sign if there
+ * are no more leading digits
+ */
+void mp_clamp(mp_int *a)
+{
+   /* decrease used while the most significant digit is
+    * zero.
+    */
+   while ((a->used > 0) && (a->dp[a->used - 1] == 0u)) {
+      --(a->used);
+   }
+
+   /* reset the sign flag if zero */
+   if (mp_iszero(a)) {
+      a->sign = MP_ZPOS;
+   }
+}
+#endif

src/libtommath/mp_clear.c

@@ -0,0 +1,20 @@
+#include "tommath_private.h"
+#ifdef MP_CLEAR_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* clear one (frees)  */
+void mp_clear(mp_int *a)
+{
+   /* only do anything if a hasn't been freed previously */
+   if (a->dp != NULL) {
+      /* free ram */
+      MP_FREE_DIGS(a->dp, a->alloc);
+
+      /* reset members to make debugging easier */
+      a->dp    = NULL;
+      a->alloc = a->used = 0;
+      a->sign  = MP_ZPOS;
+   }
+}
+#endif

src/libtommath/mp_clear_multi.c

@@ -0,0 +1,18 @@
+#include "tommath_private.h"
+#ifdef MP_CLEAR_MULTI_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+#include <stdarg.h>
+
+void mp_clear_multi(mp_int *mp, ...)
+{
+   va_list args;
+   va_start(args, mp);
+   while (mp != NULL) {
+      mp_clear(mp);
+      mp = va_arg(args, mp_int *);
+   }
+   va_end(args);
+}
+#endif

src/libtommath/mp_cmp.c

@@ -0,0 +1,21 @@
+#include "tommath_private.h"
+#ifdef MP_CMP_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* compare two ints (signed)*/
+mp_ord mp_cmp(const mp_int *a, const mp_int *b)
+{
+   /* compare based on sign */
+   if (a->sign != b->sign) {
+      return mp_isneg(a) ? MP_LT : MP_GT;
+   }
+
+   /* if negative compare opposite direction */
+   if (mp_isneg(a)) {
+      MP_EXCH(const mp_int *, a, b);
+   }
+
+   return mp_cmp_mag(a, b);
+}
+#endif

src/libtommath/mp_cmp_d.c

@@ -0,0 +1,26 @@
+#include "tommath_private.h"
+#ifdef MP_CMP_D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* compare a digit */
+mp_ord mp_cmp_d(const mp_int *a, mp_digit b)
+{
+   /* compare based on sign */
+   if (mp_isneg(a)) {
+      return MP_LT;
+   }
+
+   /* 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 a->dp[0] > b ? MP_GT : MP_LT;
+   }
+
+   return MP_EQ;
+}
+#endif

src/libtommath/mp_cmp_mag.c

@@ -0,0 +1,25 @@
+#include "tommath_private.h"
+#ifdef MP_CMP_MAG_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* compare maginitude of two ints (unsigned) */
+mp_ord mp_cmp_mag(const mp_int *a, const mp_int *b)
+{
+   int n;
+
+   /* compare based on # of non-zero digits */
+   if (a->used != b->used) {
+      return a->used > b->used ? MP_GT : MP_LT;
+   }
+
+   /* compare based on digits  */
+   for (n = a->used; n --> 0;) {
+      if (a->dp[n] != b->dp[n]) {
+         return a->dp[n] > b->dp[n] ? MP_GT : MP_LT;
+      }
+   }
+
+   return MP_EQ;
+}
+#endif

src/libtommath/mp_cnt_lsb.c

@@ -0,0 +1,38 @@
+#include "tommath_private.h"
+#ifdef MP_CNT_LSB_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+static const char lnz[16] = {
+   4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
+};
+
+/* Counts the number of lsbs which are zero before the first zero bit */
+int mp_cnt_lsb(const mp_int *a)
+{
+   int x;
+   mp_digit q;
+
+   /* easy out */
+   if (mp_iszero(a)) {
+      return 0;
+   }
+
+   /* scan lower digits until non-zero */
+   for (x = 0; (x < a->used) && (a->dp[x] == 0u); x++) {}
+   q = a->dp[x];
+   x *= MP_DIGIT_BIT;
+
+   /* now scan this digit until a 1 is found */
+   if ((q & 1u) == 0u) {
+      mp_digit p;
+      do {
+         p = q & 15u;
+         x += lnz[p];
+         q >>= 4;
+      } while (p == 0u);
+   }
+   return x;
+}
+
+#endif

src/libtommath/mp_complement.c

@@ -0,0 +1,13 @@
+#include "tommath_private.h"
+#ifdef MP_COMPLEMENT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* b = ~a */
+mp_err mp_complement(const mp_int *a, mp_int *b)
+{
+   mp_int a_ = *a;
+   a_.sign = ((a_.sign == MP_ZPOS) && !mp_iszero(a)) ? MP_NEG : MP_ZPOS;
+   return mp_sub_d(&a_, 1uL, b);
+}
+#endif

src/libtommath/mp_copy.c

@@ -0,0 +1,29 @@
+#include "tommath_private.h"
+#ifdef MP_COPY_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* copy, b = a */
+mp_err mp_copy(const mp_int *a, mp_int *b)
+{
+   mp_err err;
+
+   /* if dst == src do nothing */
+   if (a == b) {
+      return MP_OKAY;
+   }
+
+   /* grow dest */
+   if ((err = mp_grow(b, a->used)) != MP_OKAY) {
+      return err;
+   }
+
+   /* copy everything over and zero high digits */
+   s_mp_copy_digs(b->dp, a->dp, a->used);
+   s_mp_zero_digs(b->dp + a->used, b->used - a->used);
+   b->used = a->used;
+   b->sign = a->sign;
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_count_bits.c

@@ -0,0 +1,28 @@
+#include "tommath_private.h"
+#ifdef MP_COUNT_BITS_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* returns the number of bits in an int */
+int mp_count_bits(const mp_int *a)
+{
+   int     r;
+   mp_digit q;
+
+   /* shortcut */
+   if (mp_iszero(a)) {
+      return 0;
+   }
+
+   /* get number of digits and add that */
+   r = (a->used - 1) * MP_DIGIT_BIT;
+
+   /* take the last digit and count the bits in it */
+   q = a->dp[a->used - 1];
+   while (q > 0u) {
+      ++r;
+      q >>= 1u;
+   }
+   return r;
+}
+#endif

src/libtommath/mp_cutoffs.c

@@ -0,0 +1,14 @@
+#include "tommath_private.h"
+#ifdef MP_CUTOFFS_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+#ifndef MP_FIXED_CUTOFFS
+#include "tommath_cutoffs.h"
+int MP_MUL_KARATSUBA_CUTOFF = MP_DEFAULT_MUL_KARATSUBA_CUTOFF,
+    MP_SQR_KARATSUBA_CUTOFF = MP_DEFAULT_SQR_KARATSUBA_CUTOFF,
+    MP_MUL_TOOM_CUTOFF = MP_DEFAULT_MUL_TOOM_CUTOFF,
+    MP_SQR_TOOM_CUTOFF = MP_DEFAULT_SQR_TOOM_CUTOFF;
+#endif
+
+#endif

src/libtommath/mp_div.c

@@ -0,0 +1,42 @@
+#include "tommath_private.h"
+#ifdef MP_DIV_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
+{
+   mp_err err;
+
+   /* is divisor zero ? */
+   if (mp_iszero(b)) {
+      return MP_VAL;
+   }
+
+   /* if a < b then q = 0, r = a */
+   if (mp_cmp_mag(a, b) == MP_LT) {
+      if (d != NULL) {
+         if ((err = mp_copy(a, d)) != MP_OKAY) {
+            return err;
+         }
+      }
+      if (c != NULL) {
+         mp_zero(c);
+      }
+      return MP_OKAY;
+   }
+
+   if (MP_HAS(S_MP_DIV_RECURSIVE)
+       && (b->used > (2 * MP_MUL_KARATSUBA_CUTOFF))
+       && (b->used <= ((a->used/3)*2))) {
+      err = s_mp_div_recursive(a, b, c, d);
+   } else if (MP_HAS(S_MP_DIV_SCHOOL)) {
+      err = s_mp_div_school(a, b, c, d);
+   } else if (MP_HAS(S_MP_DIV_SMALL)) {
+      err = s_mp_div_small(a, b, c, d);
+   } else {
+      err = MP_VAL;
+   }
+
+   return err;
+}
+#endif

src/libtommath/mp_div_2.c

@@ -0,0 +1,40 @@
+#include "tommath_private.h"
+#ifdef MP_DIV_2_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* b = a/2 */
+mp_err mp_div_2(const mp_int *a, mp_int *b)
+{
+   mp_err err;
+   int x, oldused;
+   mp_digit r;
+
+   if ((err = mp_grow(b, a->used)) != MP_OKAY) {
+      return err;
+   }
+
+   oldused = b->used;
+   b->used = a->used;
+
+   /* carry */
+   r = 0;
+   for (x = b->used; x --> 0;) {
+      /* get the carry for the next iteration */
+      mp_digit rr = a->dp[x] & 1u;
+
+      /* shift the current digit, add in carry and store */
+      b->dp[x] = (a->dp[x] >> 1) | (r << (MP_DIGIT_BIT - 1));
+
+      /* forward carry to next iteration */
+      r = rr;
+   }
+
+   /* zero excess digits */
+   s_mp_zero_digs(b->dp + b->used, oldused - b->used);
+
+   b->sign = a->sign;
+   mp_clamp(b);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_div_2d.c

@@ -0,0 +1,61 @@
+#include "tommath_private.h"
+#ifdef MP_DIV_2D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* shift right by a certain bit count (store quotient in c, optional remainder in d) */
+mp_err mp_div_2d(const mp_int *a, int b, mp_int *c, mp_int *d)
+{
+   mp_err err;
+
+   if (b < 0) {
+      return MP_VAL;
+   }
+
+   if ((err = mp_copy(a, c)) != MP_OKAY) {
+      return err;
+   }
+
+   /* 'a' should not be used after here - it might be the same as d */
+
+   /* get the remainder */
+   if (d != NULL) {
+      if ((err = mp_mod_2d(a, b, d)) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   /* shift by as many digits in the bit count */
+   if (b >= MP_DIGIT_BIT) {
+      mp_rshd(c, b / MP_DIGIT_BIT);
+   }
+
+   /* shift any bit count < MP_DIGIT_BIT */
+   b %= MP_DIGIT_BIT;
+   if (b != 0u) {
+      int x;
+      mp_digit r, mask, shift;
+
+      /* mask */
+      mask = ((mp_digit)1 << b) - 1uL;
+
+      /* shift for lsb */
+      shift = (mp_digit)(MP_DIGIT_BIT - b);
+
+      /* carry */
+      r = 0;
+      for (x = c->used; x --> 0;) {
+         /* get the lower  bits of this word in a temp */
+         mp_digit rr = c->dp[x] & mask;
+
+         /* shift the current word and mix in the carry bits from the previous word */
+         c->dp[x] = (c->dp[x] >> b) | (r << shift);
+
+         /* set the carry to the carry bits of the current word found above */
+         r = rr;
+      }
+   }
+   mp_clamp(c);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_div_d.c

@@ -0,0 +1,84 @@
+#include "tommath_private.h"
+#ifdef MP_DIV_D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* single digit division (based on routine from MPI) */
+mp_err mp_div_d(const mp_int *a, mp_digit b, mp_int *c, mp_digit *d)
+{
+   mp_int  q;
+   mp_word w;
+   mp_err err;
+   int ix;
+
+   /* cannot divide by zero */
+   if (b == 0u) {
+      return MP_VAL;
+   }
+
+   /* quick outs */
+   if ((b == 1u) || mp_iszero(a)) {
+      if (d != NULL) {
+         *d = 0;
+      }
+      if (c != NULL) {
+         return mp_copy(a, c);
+      }
+      return MP_OKAY;
+   }
+
+   /* power of two ? */
+   if (MP_HAS(MP_DIV_2) && (b == 2u)) {
+      if (d != NULL) {
+         *d = mp_isodd(a) ? 1u : 0u;
+      }
+      return (c == NULL) ? MP_OKAY : mp_div_2(a, c);
+   }
+   if (MP_HAS(MP_DIV_2D) && MP_IS_2EXPT(b)) {
+      ix = 1;
+      while ((ix < MP_DIGIT_BIT) && (b != (((mp_digit)1)<<ix))) {
+         ix++;
+      }
+      if (d != NULL) {
+         *d = a->dp[0] & (((mp_digit)1<<(mp_digit)ix) - 1uL);
+      }
+      return (c == NULL) ? MP_OKAY : mp_div_2d(a, ix, c, NULL);
+   }
+
+   /* three? */
+   if (MP_HAS(S_MP_DIV_3) && (b == 3u)) {
+      return s_mp_div_3(a, c, d);
+   }
+
+   /* no easy answer [c'est la vie].  Just division */
+   if ((err = mp_init_size(&q, a->used)) != MP_OKAY) {
+      return err;
+   }
+
+   q.used = a->used;
+   q.sign = a->sign;
+   w = 0;
+   for (ix = a->used; ix --> 0;) {
+      mp_digit t = 0;
+      w = (w << (mp_word)MP_DIGIT_BIT) | (mp_word)a->dp[ix];
+      if (w >= b) {
+         t = (mp_digit)(w / b);
+         w -= (mp_word)t * (mp_word)b;
+      }
+      q.dp[ix] = t;
+   }
+
+   if (d != NULL) {
+      *d = (mp_digit)w;
+   }
+
+   if (c != NULL) {
+      mp_clamp(&q);
+      mp_exch(&q, c);
+   }
+   mp_clear(&q);
+
+   return MP_OKAY;
+}
+
+#endif

src/libtommath/mp_dr_is_modulus.c

@@ -0,0 +1,27 @@
+#include "tommath_private.h"
+#ifdef MP_DR_IS_MODULUS_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* determines if a number is a valid DR modulus */
+bool mp_dr_is_modulus(const mp_int *a)
+{
+   int ix;
+
+   /* must be at least two digits */
+   if (a->used < 2) {
+      return false;
+   }
+
+   /* must be of the form b**k - a [a <= b] so all
+    * but the first digit must be equal to -1 (mod b).
+    */
+   for (ix = 1; ix < a->used; ix++) {
+      if (a->dp[ix] != MP_MASK) {
+         return false;
+      }
+   }
+   return true;
+}
+
+#endif

src/libtommath/mp_dr_reduce.c

@@ -0,0 +1,68 @@
+#include "tommath_private.h"
+#ifdef MP_DR_REDUCE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* reduce "x" in place modulo "n" using the Diminished Radix algorithm.
+ *
+ * Based on algorithm from the paper
+ *
+ * "Generating Efficient Primes for Discrete Log Cryptosystems"
+ *                 Chae Hoon Lim, Pil Joong Lee,
+ *          POSTECH Information Research Laboratories
+ *
+ * The modulus must be of a special format [see manual]
+ *
+ * Has been modified to use algorithm 7.10 from the LTM book instead
+ *
+ * Input x must be in the range 0 <= x <= (n-1)**2
+ */
+mp_err mp_dr_reduce(mp_int *x, const mp_int *n, mp_digit k)
+{
+   mp_err err;
+
+   /* m = digits in modulus */
+   int m = n->used;
+
+   /* ensure that "x" has at least 2m digits */
+   if ((err = mp_grow(x, m + m)) != MP_OKAY) {
+      return err;
+   }
+
+   /* top of loop, this is where the code resumes if
+    * another reduction pass is required.
+    */
+   for (;;) {
+      int i;
+      mp_digit mu = 0;
+
+      /* compute (x mod B**m) + k * [x/B**m] inline and inplace */
+      for (i = 0; i < m; i++) {
+         mp_word r         = ((mp_word)x->dp[i + m] * (mp_word)k) + x->dp[i] + mu;
+         x->dp[i]  = (mp_digit)(r & MP_MASK);
+         mu        = (mp_digit)(r >> ((mp_word)MP_DIGIT_BIT));
+      }
+
+      /* set final carry */
+      x->dp[i] = mu;
+
+      /* zero words above m */
+      s_mp_zero_digs(x->dp + m + 1, (x->used - m) - 1);
+
+      /* clamp, sub and return */
+      mp_clamp(x);
+
+      /* if x >= n then subtract and reduce again
+       * Each successive "recursion" makes the input smaller and smaller.
+       */
+      if (mp_cmp_mag(x, n) == MP_LT) {
+         break;
+      }
+
+      if ((err = s_mp_sub(x, n, x)) != MP_OKAY) {
+         return err;
+      }
+   }
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_dr_setup.c

@@ -0,0 +1,15 @@
+#include "tommath_private.h"
+#ifdef MP_DR_SETUP_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* determines the setup value */
+void mp_dr_setup(const mp_int *a, mp_digit *d)
+{
+   /* the casts are required if MP_DIGIT_BIT is one less than
+    * the number of bits in a mp_digit [e.g. MP_DIGIT_BIT==31]
+    */
+   *d = (mp_digit)(((mp_word)1 << (mp_word)MP_DIGIT_BIT) - (mp_word)a->dp[0]);
+}
+
+#endif

src/libtommath/mp_error_to_string.c

@@ -0,0 +1,29 @@
+#include "tommath_private.h"
+#ifdef MP_ERROR_TO_STRING_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* return a char * string for a given code */
+const char *mp_error_to_string(mp_err code)
+{
+   switch (code) {
+   case MP_OKAY:
+      return "Successful";
+   case MP_ERR:
+      return "Unknown error";
+   case MP_MEM:
+      return "Out of heap";
+   case MP_VAL:
+      return "Value out of range";
+   case MP_ITER:
+      return "Max. iterations reached";
+   case MP_BUF:
+      return "Buffer overflow";
+   case MP_OVF:
+      return "Integer overflow";
+   default:
+      return "Invalid error code";
+   }
+}
+
+#endif

src/libtommath/mp_exch.c

@@ -0,0 +1,13 @@
+#include "tommath_private.h"
+#ifdef MP_EXCH_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* swap the elements of two integers, for cases where you can't simply swap the
+ * mp_int pointers around
+ */
+void mp_exch(mp_int *a, mp_int *b)
+{
+   MP_EXCH(mp_int, *a, *b);
+}
+#endif

src/libtommath/mp_expt_n.c

@@ -0,0 +1,43 @@
+#include "tommath_private.h"
+#ifdef MP_EXPT_N_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* calculate c = a**b  using a square-multiply algorithm */
+mp_err mp_expt_n(const mp_int *a, int b, mp_int *c)
+{
+   mp_err err;
+   mp_int  g;
+
+   if ((err = mp_init_copy(&g, a)) != MP_OKAY) {
+      return err;
+   }
+
+   /* set initial result */
+   mp_set(c, 1uL);
+
+   while (b > 0) {
+      /* if the bit is set multiply */
+      if ((b & 1) != 0) {
+         if ((err = mp_mul(c, &g, c)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+      }
+
+      /* square */
+      if (b > 1) {
+         if ((err = mp_sqr(&g, &g)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+      }
+
+      /* shift to next bit */
+      b >>= 1;
+   }
+
+LBL_ERR:
+   mp_clear(&g);
+   return err;
+}
+
+#endif

src/libtommath/mp_exptmod.c

@@ -0,0 +1,78 @@
+#include "tommath_private.h"
+#ifdef MP_EXPTMOD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* this is a shell function that calls either the normal or Montgomery
+ * exptmod functions.  Originally the call to the montgomery code was
+ * embedded in the normal function but that wasted alot of stack space
+ * for nothing (since 99% of the time the Montgomery code would be called)
+ */
+mp_err mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y)
+{
+   int dr;
+
+   /* modulus P must be positive */
+   if (mp_isneg(P)) {
+      return MP_VAL;
+   }
+
+   /* if exponent X is negative we have to recurse */
+   if (mp_isneg(X)) {
+      mp_int tmpG, tmpX;
+      mp_err err;
+
+      if (!MP_HAS(MP_INVMOD)) {
+         return MP_VAL;
+      }
+
+      if ((err = mp_init_multi(&tmpG, &tmpX, NULL)) != MP_OKAY) {
+         return err;
+      }
+
+      /* first compute 1/G mod P */
+      if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* now get |X| */
+      if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* and now compute (1/G)**|X| instead of G**X [X < 0] */
+      err = mp_exptmod(&tmpG, &tmpX, P, Y);
+LBL_ERR:
+      mp_clear_multi(&tmpG, &tmpX, NULL);
+      return err;
+   }
+
+   /* modified diminished radix reduction */
+   if (MP_HAS(MP_REDUCE_IS_2K_L) && MP_HAS(MP_REDUCE_2K_L) && MP_HAS(S_MP_EXPTMOD) &&
+       mp_reduce_is_2k_l(P)) {
+      return s_mp_exptmod(G, X, P, Y, 1);
+   }
+
+   /* is it a DR modulus? default to no */
+   dr = (MP_HAS(MP_DR_IS_MODULUS) && mp_dr_is_modulus(P)) ? 1 : 0;
+
+   /* if not, is it a unrestricted DR modulus? */
+   if (MP_HAS(MP_REDUCE_IS_2K) && (dr == 0)) {
+      dr = (mp_reduce_is_2k(P)) ? 2 : 0;
+   }
+
+   /* if the modulus is odd or dr != 0 use the montgomery method */
+   if (MP_HAS(S_MP_EXPTMOD_FAST) && (mp_isodd(P) || (dr != 0))) {
+      return s_mp_exptmod_fast(G, X, P, Y, dr);
+   }
+
+   /* otherwise use the generic Barrett reduction technique */
+   if (MP_HAS(S_MP_EXPTMOD)) {
+      return s_mp_exptmod(G, X, P, Y, 0);
+   }
+
+   /* no exptmod for evens */
+   return MP_VAL;
+}
+
+#endif

src/libtommath/mp_exteuclid.c

@@ -0,0 +1,72 @@
+#include "tommath_private.h"
+#ifdef MP_EXTEUCLID_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Extended euclidean algorithm of (a, b) produces
+   a*u1 + b*u2 = u3
+ */
+mp_err mp_exteuclid(const mp_int *a, const mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3)
+{
+   mp_int u1, u2, u3, v1, v2, v3, t1, t2, t3, q, tmp;
+   mp_err err;
+
+   if ((err = mp_init_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL)) != MP_OKAY) {
+      return err;
+   }
+
+   /* initialize, (u1,u2,u3) = (1,0,a) */
+   mp_set(&u1, 1uL);
+   if ((err = mp_copy(a, &u3)) != MP_OKAY)                        goto LBL_ERR;
+
+   /* initialize, (v1,v2,v3) = (0,1,b) */
+   mp_set(&v2, 1uL);
+   if ((err = mp_copy(b, &v3)) != MP_OKAY)                        goto LBL_ERR;
+
+   /* loop while v3 != 0 */
+   while (!mp_iszero(&v3)) {
+      /* q = u3/v3 */
+      if ((err = mp_div(&u3, &v3, &q, NULL)) != MP_OKAY)          goto LBL_ERR;
+
+      /* (t1,t2,t3) = (u1,u2,u3) - (v1,v2,v3)q */
+      if ((err = mp_mul(&v1, &q, &tmp)) != MP_OKAY)               goto LBL_ERR;
+      if ((err = mp_sub(&u1, &tmp, &t1)) != MP_OKAY)              goto LBL_ERR;
+      if ((err = mp_mul(&v2, &q, &tmp)) != MP_OKAY)               goto LBL_ERR;
+      if ((err = mp_sub(&u2, &tmp, &t2)) != MP_OKAY)              goto LBL_ERR;
+      if ((err = mp_mul(&v3, &q, &tmp)) != MP_OKAY)               goto LBL_ERR;
+      if ((err = mp_sub(&u3, &tmp, &t3)) != MP_OKAY)              goto LBL_ERR;
+
+      /* (u1,u2,u3) = (v1,v2,v3) */
+      if ((err = mp_copy(&v1, &u1)) != MP_OKAY)                   goto LBL_ERR;
+      if ((err = mp_copy(&v2, &u2)) != MP_OKAY)                   goto LBL_ERR;
+      if ((err = mp_copy(&v3, &u3)) != MP_OKAY)                   goto LBL_ERR;
+
+      /* (v1,v2,v3) = (t1,t2,t3) */
+      if ((err = mp_copy(&t1, &v1)) != MP_OKAY)                   goto LBL_ERR;
+      if ((err = mp_copy(&t2, &v2)) != MP_OKAY)                   goto LBL_ERR;
+      if ((err = mp_copy(&t3, &v3)) != MP_OKAY)                   goto LBL_ERR;
+   }
+
+   /* make sure U3 >= 0 */
+   if (mp_isneg(&u3)) {
+      if ((err = mp_neg(&u1, &u1)) != MP_OKAY)                    goto LBL_ERR;
+      if ((err = mp_neg(&u2, &u2)) != MP_OKAY)                    goto LBL_ERR;
+      if ((err = mp_neg(&u3, &u3)) != MP_OKAY)                    goto LBL_ERR;
+   }
+
+   /* copy result out */
+   if (U1 != NULL) {
+      mp_exch(U1, &u1);
+   }
+   if (U2 != NULL) {
+      mp_exch(U2, &u2);
+   }
+   if (U3 != NULL) {
+      mp_exch(U3, &u3);
+   }
+
+LBL_ERR:
+   mp_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL);
+   return err;
+}
+#endif

src/libtommath/mp_fread.c

@@ -0,0 +1,66 @@
+#include "tommath_private.h"
+#ifdef MP_FREAD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+#ifndef MP_NO_FILE
+/* read a bigint from a file stream in ASCII */
+mp_err mp_fread(mp_int *a, int radix, FILE *stream)
+{
+   mp_err err;
+   mp_sign sign = MP_ZPOS;
+   int ch;
+
+   /* make sure the radix is ok */
+   if ((radix < 2) || (radix > 64)) {
+      return MP_VAL;
+   }
+
+   /* if first digit is - then set negative */
+   ch = fgetc(stream);
+   if (ch == (int)'-') {
+      sign = MP_NEG;
+      ch = fgetc(stream);
+   }
+
+   /* no digits, return error */
+   if (ch == EOF) {
+      return MP_ERR;
+   }
+
+   /* clear a */
+   mp_zero(a);
+
+   do {
+      uint8_t y;
+      unsigned pos;
+      ch = (radix <= 36) ? MP_TOUPPER(ch) : ch;
+      pos = (unsigned)(ch - (int)'+');
+      if (MP_RADIX_MAP_REVERSE_SIZE <= pos) {
+         break;
+      }
+
+      y = s_mp_radix_map_reverse[pos];
+
+      if (y >= radix) {
+         break;
+      }
+
+      /* shift up and add */
+      if ((err = mp_mul_d(a, (mp_digit)radix, a)) != MP_OKAY) {
+         return err;
+      }
+      if ((err = mp_add_d(a, y, a)) != MP_OKAY) {
+         return err;
+      }
+   } while ((ch = fgetc(stream)) != EOF);
+
+   if (!mp_iszero(a)) {
+      a->sign = sign;
+   }
+
+   return MP_OKAY;
+}
+#endif
+
+#endif

src/libtommath/mp_from_sbin.c

@@ -0,0 +1,21 @@
+#include "tommath_private.h"
+#ifdef MP_FROM_SBIN_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* read signed bin, big endian, first byte is 0==positive or 1==negative */
+mp_err mp_from_sbin(mp_int *a, const uint8_t *buf, size_t size)
+{
+   mp_err err;
+
+   /* read magnitude */
+   if ((err = mp_from_ubin(a, buf + 1, size - 1u)) != MP_OKAY) {
+      return err;
+   }
+
+   /* first byte is 0 for positive, non-zero for negative */
+   a->sign = (buf[0] != (uint8_t)0) ? MP_NEG : MP_ZPOS;
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_from_ubin.c

@@ -0,0 +1,30 @@
+#include "tommath_private.h"
+#ifdef MP_FROM_UBIN_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* reads a uint8_t array, assumes the msb is stored first [big endian] */
+mp_err mp_from_ubin(mp_int *a, const uint8_t *buf, size_t size)
+{
+   mp_err err;
+
+   /* make sure there are at least two digits */
+   if ((err = mp_grow(a, 2)) != MP_OKAY) {
+      return err;
+   }
+
+   /* zero the int */
+   mp_zero(a);
+
+   /* read the bytes in */
+   while (size-- > 0u) {
+      if ((err = mp_mul_2d(a, 8, a)) != MP_OKAY) {
+         return err;
+      }
+      a->dp[0] |= *buf++;
+      a->used += 1;
+   }
+   mp_clamp(a);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_fwrite.c

@@ -0,0 +1,33 @@
+#include "tommath_private.h"
+#ifdef MP_FWRITE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+#ifndef MP_NO_FILE
+mp_err mp_fwrite(const mp_int *a, int radix, FILE *stream)
+{
+   char *buf;
+   mp_err err;
+   size_t size, written;
+
+   if ((err = mp_radix_size_overestimate(a, radix, &size)) != MP_OKAY) {
+      return err;
+   }
+
+   buf = (char *) MP_MALLOC(size);
+   if (buf == NULL) {
+      return MP_MEM;
+   }
+
+   if ((err = mp_to_radix(a, buf, size, &written, radix)) == MP_OKAY) {
+      if (fwrite(buf, written, 1uL, stream) != 1uL) {
+         err = MP_ERR;
+      }
+   }
+
+   MP_FREE_BUF(buf, size);
+   return err;
+}
+#endif
+
+#endif

src/libtommath/mp_gcd.c

@@ -0,0 +1,92 @@
+#include "tommath_private.h"
+#ifdef MP_GCD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Greatest Common Divisor using the binary method */
+mp_err mp_gcd(const mp_int *a, const mp_int *b, mp_int *c)
+{
+   mp_int  u, v;
+   int     k, u_lsb, v_lsb;
+   mp_err err;
+
+   /* either zero than gcd is the largest */
+   if (mp_iszero(a)) {
+      return mp_abs(b, c);
+   }
+   if (mp_iszero(b)) {
+      return mp_abs(a, c);
+   }
+
+   /* get copies of a and b we can modify */
+   if ((err = mp_init_copy(&u, a)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = mp_init_copy(&v, b)) != MP_OKAY) {
+      goto LBL_U;
+   }
+
+   /* must be positive for the remainder of the algorithm */
+   u.sign = v.sign = MP_ZPOS;
+
+   /* B1.  Find the common power of two for u and v */
+   u_lsb = mp_cnt_lsb(&u);
+   v_lsb = mp_cnt_lsb(&v);
+   k     = MP_MIN(u_lsb, v_lsb);
+
+   if (k > 0) {
+      /* divide the power of two out */
+      if ((err = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
+         goto LBL_V;
+      }
+
+      if ((err = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
+         goto LBL_V;
+      }
+   }
+
+   /* divide any remaining factors of two out */
+   if (u_lsb != k) {
+      if ((err = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
+         goto LBL_V;
+      }
+   }
+
+   if (v_lsb != k) {
+      if ((err = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
+         goto LBL_V;
+      }
+   }
+
+   while (!mp_iszero(&v)) {
+      /* make sure v is the largest */
+      if (mp_cmp_mag(&u, &v) == MP_GT) {
+         /* swap u and v to make sure v is >= u */
+         mp_exch(&u, &v);
+      }
+
+      /* subtract smallest from largest */
+      if ((err = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
+         goto LBL_V;
+      }
+
+      /* Divide out all factors of two */
+      if ((err = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
+         goto LBL_V;
+      }
+   }
+
+   /* multiply by 2**k which we divided out at the beginning */
+   if ((err = mp_mul_2d(&u, k, c)) != MP_OKAY) {
+      goto LBL_V;
+   }
+   c->sign = MP_ZPOS;
+   err = MP_OKAY;
+LBL_V:
+   mp_clear(&u);
+LBL_U:
+   mp_clear(&v);
+   return err;
+}
+#endif

src/libtommath/mp_get_double.c

@@ -0,0 +1,18 @@
+#include "tommath_private.h"
+#ifdef MP_GET_DOUBLE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+double mp_get_double(const mp_int *a)
+{
+   int i;
+   double d = 0.0, fac = 1.0;
+   for (i = 0; i < MP_DIGIT_BIT; ++i) {
+      fac *= 2.0;
+   }
+   for (i = a->used; i --> 0;) {
+      d = (d * fac) + (double)a->dp[i];
+   }
+   return mp_isneg(a) ? -d : d;
+}
+#endif

src/libtommath/mp_get_i32.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_GET_I32_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_GET_SIGNED(mp_get_i32, mp_get_mag_u32, int32_t, uint32_t)
+#endif

src/libtommath/mp_get_i64.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_GET_I64_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_GET_SIGNED(mp_get_i64, mp_get_mag_u64, int64_t, uint64_t)
+#endif

src/libtommath/mp_get_l.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_GET_L_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_GET_SIGNED(mp_get_l, mp_get_mag_ul, long, unsigned long)
+#endif

src/libtommath/mp_get_mag_u32.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_GET_MAG_U32_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_GET_MAG(mp_get_mag_u32, uint32_t)
+#endif

src/libtommath/mp_get_mag_u64.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_GET_MAG_U64_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_GET_MAG(mp_get_mag_u64, uint64_t)
+#endif

src/libtommath/mp_get_mag_ul.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_GET_MAG_UL_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_GET_MAG(mp_get_mag_ul, unsigned long)
+#endif

src/libtommath/mp_grow.c

@@ -0,0 +1,40 @@
+#include "tommath_private.h"
+#ifdef MP_GROW_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* grow as required */
+mp_err mp_grow(mp_int *a, int size)
+{
+   /* if the alloc size is smaller alloc more ram */
+   if (a->alloc < size) {
+      mp_digit *dp;
+
+      if (size > MP_MAX_DIGIT_COUNT) {
+         return MP_OVF;
+      }
+
+      /* reallocate the array a->dp
+       *
+       * We store the return in a temporary variable
+       * in case the operation failed we don't want
+       * to overwrite the dp member of a.
+       */
+      dp = (mp_digit *) MP_REALLOC(a->dp,
+                                   (size_t)a->alloc * sizeof(mp_digit),
+                                   (size_t)size * sizeof(mp_digit));
+      if (dp == NULL) {
+         /* reallocation failed but "a" is still valid [can be freed] */
+         return MP_MEM;
+      }
+
+      /* reallocation succeeded so set a->dp */
+      a->dp = dp;
+
+      /* zero excess digits */
+      s_mp_zero_digs(a->dp + a->alloc, size - a->alloc);
+      a->alloc = size;
+   }
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_init.c

@@ -0,0 +1,23 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* init a new mp_int */
+mp_err mp_init(mp_int *a)
+{
+   /* allocate memory required and clear it */
+   a->dp = (mp_digit *) MP_CALLOC((size_t)MP_DEFAULT_DIGIT_COUNT, sizeof(mp_digit));
+   if (a->dp == NULL) {
+      return MP_MEM;
+   }
+
+   /* set the used to zero, allocated digits to the default precision
+    * and sign to positive */
+   a->used  = 0;
+   a->alloc = MP_DEFAULT_DIGIT_COUNT;
+   a->sign  = MP_ZPOS;
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_init_copy.c

@@ -0,0 +1,21 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_COPY_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* creates "a" then copies b into it */
+mp_err mp_init_copy(mp_int *a, const mp_int *b)
+{
+   mp_err     err;
+
+   if ((err = mp_init_size(a, b->used)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = mp_copy(b, a)) != MP_OKAY) {
+      mp_clear(a);
+   }
+
+   return err;
+}
+#endif

src/libtommath/mp_init_i32.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_I32_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_INIT_INT(mp_init_i32, mp_set_i32, int32_t)
+#endif

src/libtommath/mp_init_i64.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_I64_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_INIT_INT(mp_init_i64, mp_set_i64, int64_t)
+#endif

src/libtommath/mp_init_l.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_L_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_INIT_INT(mp_init_l, mp_set_l, long)
+#endif

src/libtommath/mp_init_multi.c

@@ -0,0 +1,41 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_MULTI_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+#include <stdarg.h>
+
+mp_err mp_init_multi(mp_int *mp, ...)
+{
+   mp_err err = MP_OKAY;
+   int n = 0;                 /* Number of ok inits */
+   mp_int *cur_arg = mp;
+   va_list args;
+
+   va_start(args, mp);        /* init args to next argument from caller */
+   while (cur_arg != NULL) {
+      err = mp_init(cur_arg);
+      if (err != MP_OKAY) {
+         /* Oops - error! Back-track and mp_clear what we already
+            succeeded in init-ing, then return error.
+         */
+         va_list clean_args;
+
+         /* now start cleaning up */
+         cur_arg = mp;
+         va_start(clean_args, mp);
+         while (n-- != 0) {
+            mp_clear(cur_arg);
+            cur_arg = va_arg(clean_args, mp_int *);
+         }
+         va_end(clean_args);
+         break;
+      }
+      n++;
+      cur_arg = va_arg(args, mp_int *);
+   }
+   va_end(args);
+   return err;
+}
+
+#endif

src/libtommath/mp_init_set.c

@@ -0,0 +1,16 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_SET_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* initialize and set a digit */
+mp_err mp_init_set(mp_int *a, mp_digit b)
+{
+   mp_err err;
+   if ((err = mp_init(a)) != MP_OKAY) {
+      return err;
+   }
+   mp_set(a, b);
+   return err;
+}
+#endif

src/libtommath/mp_init_size.c

@@ -0,0 +1,28 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_SIZE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* init an mp_init for a given size */
+mp_err mp_init_size(mp_int *a, int size)
+{
+   size = MP_MAX(MP_MIN_DIGIT_COUNT, size);
+
+   if (size > MP_MAX_DIGIT_COUNT) {
+      return MP_OVF;
+   }
+
+   /* alloc mem */
+   a->dp = (mp_digit *) MP_CALLOC((size_t)size, sizeof(mp_digit));
+   if (a->dp == NULL) {
+      return MP_MEM;
+   }
+
+   /* set the members */
+   a->used  = 0;
+   a->alloc = size;
+   a->sign  = MP_ZPOS;
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_init_u32.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_U32_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_INIT_INT(mp_init_u32, mp_set_u32, uint32_t)
+#endif

src/libtommath/mp_init_u64.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_U64_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_INIT_INT(mp_init_u64, mp_set_u64, uint64_t)
+#endif

src/libtommath/mp_init_ul.c

@@ -0,0 +1,7 @@
+#include "tommath_private.h"
+#ifdef MP_INIT_UL_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+MP_INIT_INT(mp_init_ul, mp_set_ul, unsigned long)
+#endif

src/libtommath/mp_invmod.c

@@ -0,0 +1,29 @@
+#include "tommath_private.h"
+#ifdef MP_INVMOD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* hac 14.61, pp608 */
+mp_err mp_invmod(const mp_int *a, const mp_int *b, mp_int *c)
+{
+   /* for all n in N and n > 0, n = 0 mod 1 */
+   if (!mp_isneg(a) && mp_cmp_d(b, 1uL) == MP_EQ) {
+      mp_zero(c);
+      return MP_OKAY;
+   }
+
+   /* b cannot be negative and has to be >1 */
+   if (mp_isneg(b) || (mp_cmp_d(b, 1uL) != MP_GT)) {
+      return MP_VAL;
+   }
+
+   /* if the modulus is odd we can use a faster routine instead */
+   if (MP_HAS(S_MP_INVMOD_ODD) && mp_isodd(b)) {
+      return s_mp_invmod_odd(a, b, c);
+   }
+
+   return MP_HAS(S_MP_INVMOD)
+          ? s_mp_invmod(a, b, c)
+          : MP_VAL;
+}
+#endif

src/libtommath/mp_is_square.c

@@ -0,0 +1,93 @@
+#include "tommath_private.h"
+#ifdef MP_IS_SQUARE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Check if remainders are possible squares - fast exclude non-squares */
+static const char rem_128[128] = {
+   0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
+};
+
+static const char rem_105[105] = {
+   0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
+   0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
+   0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
+   1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
+   0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
+   1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
+   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
+};
+
+/* Store non-zero to ret if arg is square, and zero if not */
+mp_err mp_is_square(const mp_int *arg, bool *ret)
+{
+   mp_err   err;
+   mp_digit c;
+   mp_int   t;
+   uint32_t r;
+
+   /* Default to Non-square :) */
+   *ret = false;
+
+   if (mp_isneg(arg)) {
+      return MP_VAL;
+   }
+
+   if (mp_iszero(arg)) {
+      return MP_OKAY;
+   }
+
+   /* First check mod 128 (suppose that MP_DIGIT_BIT is at least 7) */
+   if (rem_128[127u & arg->dp[0]] == (char)1) {
+      return MP_OKAY;
+   }
+
+   /* Next check mod 105 (3*5*7) */
+   if ((err = mp_mod_d(arg, 105uL, &c)) != MP_OKAY) {
+      return err;
+   }
+   if (rem_105[c] == (char)1) {
+      return MP_OKAY;
+   }
+
+
+   if ((err = mp_init_u32(&t, 11u*13u*17u*19u*23u*29u*31u)) != MP_OKAY) {
+      return err;
+   }
+   if ((err = mp_mod(arg, &t, &t)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   r = mp_get_u32(&t);
+   /* Check for other prime modules, note it's not an ERROR but we must
+    * free "t" so the easiest way is to goto LBL_ERR.  We know that err
+    * is already equal to MP_OKAY from the mp_mod call
+    */
+   if (((1uL<<(r%11uL)) & 0x5C4uL) != 0uL)         goto LBL_ERR;
+   if (((1uL<<(r%13uL)) & 0x9E4uL) != 0uL)         goto LBL_ERR;
+   if (((1uL<<(r%17uL)) & 0x5CE8uL) != 0uL)        goto LBL_ERR;
+   if (((1uL<<(r%19uL)) & 0x4F50CuL) != 0uL)       goto LBL_ERR;
+   if (((1uL<<(r%23uL)) & 0x7ACCA0uL) != 0uL)      goto LBL_ERR;
+   if (((1uL<<(r%29uL)) & 0xC2EDD0CuL) != 0uL)     goto LBL_ERR;
+   if (((1uL<<(r%31uL)) & 0x6DE2B848uL) != 0uL)    goto LBL_ERR;
+
+   /* Final check - is sqr(sqrt(arg)) == arg ? */
+   if ((err = mp_sqrt(arg, &t)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   if ((err = mp_sqr(&t, &t)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   *ret = (mp_cmp_mag(&t, arg) == MP_EQ);
+LBL_ERR:
+   mp_clear(&t);
+   return err;
+}
+#endif

src/libtommath/mp_kronecker.c

@@ -0,0 +1,129 @@
+#include "tommath_private.h"
+#ifdef MP_KRONECKER_C
+
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/*
+   Kronecker symbol (a|p)
+   Straightforward implementation of algorithm 1.4.10 in
+   Henri Cohen: "A Course in Computational Algebraic Number Theory"
+
+   @book{cohen2013course,
+     title={A course in computational algebraic number theory},
+     author={Cohen, Henri},
+     volume={138},
+     year={2013},
+     publisher={Springer Science \& Business Media}
+    }
+ */
+mp_err mp_kronecker(const mp_int *a, const mp_int *p, int *c)
+{
+   mp_int a1, p1, r;
+   mp_err err;
+   int v, k;
+
+   static const char table[] = {0, 1, 0, -1, 0, -1, 0, 1};
+
+   if (mp_iszero(p)) {
+      if ((a->used == 1) && (a->dp[0] == 1u)) {
+         *c = 1;
+      } else {
+         *c = 0;
+      }
+      return MP_OKAY;
+   }
+
+   if (mp_iseven(a) && mp_iseven(p)) {
+      *c = 0;
+      return MP_OKAY;
+   }
+
+   if ((err = mp_init_copy(&a1, a)) != MP_OKAY) {
+      return err;
+   }
+   if ((err = mp_init_copy(&p1, p)) != MP_OKAY) {
+      goto LBL_KRON_0;
+   }
+
+   v = mp_cnt_lsb(&p1);
+   if ((err = mp_div_2d(&p1, v, &p1, NULL)) != MP_OKAY) {
+      goto LBL_KRON_1;
+   }
+
+   if ((v & 1) == 0) {
+      k = 1;
+   } else {
+      k = table[a->dp[0] & 7u];
+   }
+
+   if (mp_isneg(&p1)) {
+      p1.sign = MP_ZPOS;
+      if (mp_isneg(&a1)) {
+         k = -k;
+      }
+   }
+
+   if ((err = mp_init(&r)) != MP_OKAY) {
+      goto LBL_KRON_1;
+   }
+
+   for (;;) {
+      if (mp_iszero(&a1)) {
+         if (mp_cmp_d(&p1, 1uL) == MP_EQ) {
+            *c = k;
+            goto LBL_KRON;
+         } else {
+            *c = 0;
+            goto LBL_KRON;
+         }
+      }
+
+      v = mp_cnt_lsb(&a1);
+      if ((err = mp_div_2d(&a1, v, &a1, NULL)) != MP_OKAY) {
+         goto LBL_KRON;
+      }
+
+      if ((v & 1) == 1) {
+         k = k * table[p1.dp[0] & 7u];
+      }
+
+      if (mp_isneg(&a1)) {
+         /*
+          * Compute k = (-1)^((a1)*(p1-1)/4) * k
+          * a1.dp[0] + 1 cannot overflow because the MSB
+          * of the type mp_digit is not set by definition
+          */
+         if (((a1.dp[0] + 1u) & p1.dp[0] & 2u) != 0u) {
+            k = -k;
+         }
+      } else {
+         /* compute k = (-1)^((a1-1)*(p1-1)/4) * k */
+         if ((a1.dp[0] & p1.dp[0] & 2u) != 0u) {
+            k = -k;
+         }
+      }
+
+      if ((err = mp_copy(&a1, &r)) != MP_OKAY) {
+         goto LBL_KRON;
+      }
+      r.sign = MP_ZPOS;
+      if ((err = mp_mod(&p1, &r, &a1)) != MP_OKAY) {
+         goto LBL_KRON;
+      }
+      if ((err = mp_copy(&r, &p1)) != MP_OKAY) {
+         goto LBL_KRON;
+      }
+   }
+
+LBL_KRON:
+   mp_clear(&r);
+LBL_KRON_1:
+   mp_clear(&p1);
+LBL_KRON_0:
+   mp_clear(&a1);
+
+   return err;
+}
+
+#endif

src/libtommath/mp_lcm.c

@@ -0,0 +1,44 @@
+#include "tommath_private.h"
+#ifdef MP_LCM_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* computes least common multiple as |a*b|/(a, b) */
+mp_err mp_lcm(const mp_int *a, const mp_int *b, mp_int *c)
+{
+   mp_err  err;
+   mp_int  t1, t2;
+
+
+   if ((err = mp_init_multi(&t1, &t2, NULL)) != MP_OKAY) {
+      return err;
+   }
+
+   /* t1 = get the GCD of the two inputs */
+   if ((err = mp_gcd(a, b, &t1)) != MP_OKAY) {
+      goto LBL_T;
+   }
+
+   /* divide the smallest by the GCD */
+   if (mp_cmp_mag(a, b) == MP_LT) {
+      /* store quotient in t2 such that t2 * b is the LCM */
+      if ((err = mp_div(a, &t1, &t2, NULL)) != MP_OKAY) {
+         goto LBL_T;
+      }
+      err = mp_mul(b, &t2, c);
+   } else {
+      /* store quotient in t2 such that t2 * a is the LCM */
+      if ((err = mp_div(b, &t1, &t2, NULL)) != MP_OKAY) {
+         goto LBL_T;
+      }
+      err = mp_mul(a, &t2, c);
+   }
+
+   /* fix the sign to positive */
+   c->sign = MP_ZPOS;
+
+LBL_T:
+   mp_clear_multi(&t1, &t2, NULL);
+   return err;
+}
+#endif

src/libtommath/mp_log_n.c

@@ -0,0 +1,29 @@
+#include "tommath_private.h"
+#ifdef MP_LOG_N_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+mp_err mp_log_n(const mp_int *a, int base, int *c)
+{
+   if (mp_isneg(a) || mp_iszero(a) || (base < 2) || (unsigned)base > (unsigned)MP_DIGIT_MAX) {
+      return MP_VAL;
+   }
+
+   if (MP_HAS(S_MP_LOG_2EXPT) && MP_IS_2EXPT((mp_digit)base)) {
+      *c = s_mp_log_2expt(a, (mp_digit)base);
+      return MP_OKAY;
+   }
+
+   if (MP_HAS(S_MP_LOG_D) && (a->used == 1)) {
+      *c = s_mp_log_d((mp_digit)base, a->dp[0]);
+      return MP_OKAY;
+   }
+
+   if (MP_HAS(S_MP_LOG)) {
+      return s_mp_log(a, (mp_digit)base, c);
+   }
+
+   return MP_VAL;
+}
+
+#endif

src/libtommath/mp_lshd.c

@@ -0,0 +1,42 @@
+#include "tommath_private.h"
+#ifdef MP_LSHD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* shift left a certain amount of digits */
+mp_err mp_lshd(mp_int *a, int b)
+{
+   mp_err err;
+   int x;
+
+   /* if its less than zero return */
+   if (b <= 0) {
+      return MP_OKAY;
+   }
+   /* no need to shift 0 around */
+   if (mp_iszero(a)) {
+      return MP_OKAY;
+   }
+
+   /* grow to fit the new digits */
+   if ((err = mp_grow(a, a->used + b)) != MP_OKAY) {
+      return err;
+   }
+
+   /* increment the used by the shift amount then copy upwards */
+   a->used += b;
+
+   /* much like mp_rshd this is implemented using a sliding window
+    * except the window goes the otherway around.  Copying from
+    * the bottom to the top.  see mp_rshd.c for more info.
+    */
+   for (x = a->used; x --> b;) {
+      a->dp[x] = a->dp[x - b];
+   }
+
+   /* zero the lower digits */
+   s_mp_zero_digs(a->dp, b);
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_mod.c

@@ -0,0 +1,15 @@
+#include "tommath_private.h"
+#ifdef MP_MOD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* c = a mod b, 0 <= c < b if b > 0, b < c <= 0 if b < 0 */
+mp_err mp_mod(const mp_int *a, const mp_int *b, mp_int *c)
+{
+   mp_err err;
+   if ((err = mp_div(a, b, NULL, c)) != MP_OKAY) {
+      return err;
+   }
+   return mp_iszero(c) || (c->sign == b->sign) ? MP_OKAY : mp_add(b, c, c);
+}
+#endif

src/libtommath/mp_mod_2d.c

@@ -0,0 +1,40 @@
+#include "tommath_private.h"
+#ifdef MP_MOD_2D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* calc a value mod 2**b */
+mp_err mp_mod_2d(const mp_int *a, int b, mp_int *c)
+{
+   int x;
+   mp_err err;
+
+   if (b < 0) {
+      return MP_VAL;
+   }
+
+   if (b == 0) {
+      mp_zero(c);
+      return MP_OKAY;
+   }
+
+   /* if the modulus is larger than the value than return */
+   if (b >= (a->used * MP_DIGIT_BIT)) {
+      return mp_copy(a, c);
+   }
+
+   if ((err = mp_copy(a, c)) != MP_OKAY) {
+      return err;
+   }
+
+   /* zero digits above the last digit of the modulus */
+   x = (b / MP_DIGIT_BIT) + (((b % MP_DIGIT_BIT) == 0) ? 0 : 1);
+   s_mp_zero_digs(c->dp + x, c->used - x);
+
+   /* clear the digit that is not completely outside/inside the modulus */
+   c->dp[b / MP_DIGIT_BIT] &=
+      ((mp_digit)1 << (mp_digit)(b % MP_DIGIT_BIT)) - (mp_digit)1;
+   mp_clamp(c);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_montgomery_calc_normalization.c

@@ -0,0 +1,43 @@
+#include "tommath_private.h"
+#ifdef MP_MONTGOMERY_CALC_NORMALIZATION_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/*
+ * shifts with subtractions when the result is greater than b.
+ *
+ * The method is slightly modified to shift B unconditionally upto just under
+ * the leading bit of b.  This saves alot of multiple precision shifting.
+ */
+mp_err mp_montgomery_calc_normalization(mp_int *a, const mp_int *b)
+{
+   int    x, bits;
+   mp_err err;
+
+   /* how many bits of last digit does b use */
+   bits = mp_count_bits(b) % MP_DIGIT_BIT;
+
+   if (b->used > 1) {
+      if ((err = mp_2expt(a, ((b->used - 1) * MP_DIGIT_BIT) + bits - 1)) != MP_OKAY) {
+         return err;
+      }
+   } else {
+      mp_set(a, 1uL);
+      bits = 1;
+   }
+
+   /* now compute C = A * B mod b */
+   for (x = bits - 1; x < (int)MP_DIGIT_BIT; x++) {
+      if ((err = mp_mul_2(a, a)) != MP_OKAY) {
+         return err;
+      }
+      if (mp_cmp_mag(a, b) != MP_LT) {
+         if ((err = s_mp_sub(a, b, a)) != MP_OKAY) {
+            return err;
+         }
+      }
+   }
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_montgomery_reduce.c

@@ -0,0 +1,89 @@
+#include "tommath_private.h"
+#ifdef MP_MONTGOMERY_REDUCE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* computes xR**-1 == x (mod N) via Montgomery Reduction */
+mp_err mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho)
+{
+   mp_err err;
+   int ix, digs;
+
+   /* can the fast reduction [comba] method be used?
+    *
+    * Note that unlike in mul you're safely allowed *less*
+    * than the available columns [255 per default] since carries
+    * are fixed up in the inner loop.
+    */
+   digs = (n->used * 2) + 1;
+   if ((digs < MP_WARRAY) &&
+       (x->used <= MP_WARRAY) &&
+       (n->used < MP_MAX_COMBA)) {
+      return s_mp_montgomery_reduce_comba(x, n, rho);
+   }
+
+   /* grow the input as required */
+   if ((err = mp_grow(x, digs)) != MP_OKAY) {
+      return err;
+   }
+   x->used = digs;
+
+   for (ix = 0; ix < n->used; ix++) {
+      int iy;
+      mp_digit u, mu;
+
+      /* mu = ai * rho mod b
+       *
+       * The value of rho must be precalculated via
+       * montgomery_setup() such that
+       * it equals -1/n0 mod b this allows the
+       * following inner loop to reduce the
+       * input one digit at a time
+       */
+      mu = (mp_digit)(((mp_word)x->dp[ix] * (mp_word)rho) & MP_MASK);
+
+      /* a = a + mu * m * b**i */
+
+      /* Multiply and add in place */
+      u = 0;
+      for (iy = 0; iy < n->used; iy++) {
+         /* compute product and sum */
+         mp_word r = ((mp_word)mu * (mp_word)n->dp[iy]) +
+                     (mp_word)u + (mp_word)x->dp[ix + iy];
+
+         /* get carry */
+         u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
+
+         /* fix digit */
+         x->dp[ix + iy] = (mp_digit)(r & (mp_word)MP_MASK);
+      }
+      /* At this point the ix'th digit of x should be zero */
+
+      /* propagate carries upwards as required*/
+      while (u != 0u) {
+         x->dp[ix + iy]   += u;
+         u        = x->dp[ix + iy] >> MP_DIGIT_BIT;
+         x->dp[ix + iy] &= MP_MASK;
+         ++iy;
+      }
+   }
+
+   /* at this point the n.used'th least
+    * significant digits of x are all zero
+    * which means we can shift x to the
+    * right by n.used digits and the
+    * residue is unchanged.
+    */
+
+   /* x = x/b**n.used */
+   mp_clamp(x);
+   mp_rshd(x, n->used);
+
+   /* if x >= n then x = x - n */
+   if (mp_cmp_mag(x, n) != MP_LT) {
+      return s_mp_sub(x, n, x);
+   }
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_montgomery_setup.c

@@ -0,0 +1,40 @@
+#include "tommath_private.h"
+#ifdef MP_MONTGOMERY_SETUP_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* setups the montgomery reduction stuff */
+mp_err mp_montgomery_setup(const mp_int *n, mp_digit *rho)
+{
+   mp_digit x, b;
+
+   /* fast inversion mod 2**k
+    *
+    * Based on the fact that
+    *
+    * XA = 1 (mod 2**n)  =>  (X(2-XA)) A = 1 (mod 2**2n)
+    *                    =>  2*X*A - X*X*A*A = 1
+    *                    =>  2*(1) - (1)     = 1
+    */
+   b = n->dp[0];
+
+   if ((b & 1u) == 0u) {
+      return MP_VAL;
+   }
+
+   x = (((b + 2u) & 4u) << 1) + b; /* here x*a==1 mod 2**4 */
+   x *= 2u - (b * x);              /* here x*a==1 mod 2**8 */
+   x *= 2u - (b * x);              /* here x*a==1 mod 2**16 */
+#if defined(MP_64BIT) || !(defined(MP_16BIT))
+   x *= 2u - (b * x);              /* here x*a==1 mod 2**32 */
+#endif
+#ifdef MP_64BIT
+   x *= 2u - (b * x);              /* here x*a==1 mod 2**64 */
+#endif
+
+   /* rho = -1/m mod b */
+   *rho = (mp_digit)(((mp_word)1 << (mp_word)MP_DIGIT_BIT) - x) & MP_MASK;
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_mul.c

@@ -0,0 +1,68 @@
+#include "tommath_private.h"
+#ifdef MP_MUL_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* high level multiplication (handles sign) */
+mp_err mp_mul(const mp_int *a, const mp_int *b, mp_int *c)
+{
+   mp_err err;
+   int min = MP_MIN(a->used, b->used),
+       max = MP_MAX(a->used, b->used),
+       digs = a->used + b->used + 1;
+   bool neg = (a->sign != b->sign);
+
+   if ((a == b) &&
+       MP_HAS(S_MP_SQR_TOOM) && /* use Toom-Cook? */
+       (a->used >= MP_SQR_TOOM_CUTOFF)) {
+      err = s_mp_sqr_toom(a, c);
+   } else if ((a == b) &&
+              MP_HAS(S_MP_SQR_KARATSUBA) &&  /* Karatsuba? */
+              (a->used >= MP_SQR_KARATSUBA_CUTOFF)) {
+      err = s_mp_sqr_karatsuba(a, c);
+   } else if ((a == b) &&
+              MP_HAS(S_MP_SQR_COMBA) && /* can we use the fast comba multiplier? */
+              (((a->used * 2) + 1) < MP_WARRAY) &&
+              (a->used < (MP_MAX_COMBA / 2))) {
+      err = s_mp_sqr_comba(a, c);
+   } else if ((a == b) &&
+              MP_HAS(S_MP_SQR)) {
+      err = s_mp_sqr(a, c);
+   } else if (MP_HAS(S_MP_MUL_BALANCE) &&
+              /* Check sizes. The smaller one needs to be larger than the Karatsuba cut-off.
+               * The bigger one needs to be at least about one MP_MUL_KARATSUBA_CUTOFF bigger
+               * to make some sense, but it depends on architecture, OS, position of the
+               * stars... so YMMV.
+               * Using it to cut the input into slices small enough for s_mp_mul_comba
+               * was actually slower on the author's machine, but YMMV.
+               */
+              (min >= MP_MUL_KARATSUBA_CUTOFF) &&
+              ((max / 2) >= MP_MUL_KARATSUBA_CUTOFF) &&
+              /* Not much effect was observed below a ratio of 1:2, but again: YMMV. */
+              (max >= (2 * min))) {
+      err = s_mp_mul_balance(a,b,c);
+   } else if (MP_HAS(S_MP_MUL_TOOM) &&
+              (min >= MP_MUL_TOOM_CUTOFF)) {
+      err = s_mp_mul_toom(a, b, c);
+   } else if (MP_HAS(S_MP_MUL_KARATSUBA) &&
+              (min >= MP_MUL_KARATSUBA_CUTOFF)) {
+      err = s_mp_mul_karatsuba(a, b, c);
+   } else if (MP_HAS(S_MP_MUL_COMBA) &&
+              /* can we use the fast multiplier?
+               *
+               * The fast multiplier can be used if the output will
+               * have less than MP_WARRAY digits and the number of
+               * digits won't affect carry propagation
+               */
+              (digs < MP_WARRAY) &&
+              (min <= MP_MAX_COMBA)) {
+      err = s_mp_mul_comba(a, b, c, digs);
+   } else if (MP_HAS(S_MP_MUL)) {
+      err = s_mp_mul(a, b, c, digs);
+   } else {
+      err = MP_VAL;
+   }
+   c->sign = ((c->used > 0) && neg) ? MP_NEG : MP_ZPOS;
+   return err;
+}
+#endif

src/libtommath/mp_mul_2.c

@@ -0,0 +1,53 @@
+#include "tommath_private.h"
+#ifdef MP_MUL_2_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* b = a*2 */
+mp_err mp_mul_2(const mp_int *a, mp_int *b)
+{
+   mp_err err;
+   int x, oldused;
+   mp_digit r;
+
+   /* grow to accomodate result */
+   if ((err = mp_grow(b, a->used + 1)) != MP_OKAY) {
+      return err;
+   }
+
+   oldused = b->used;
+   b->used = a->used;
+
+   /* carry */
+   r = 0;
+   for (x = 0; x < a->used; x++) {
+
+      /* get what will be the *next* carry bit from the
+       * MSB of the current digit
+       */
+      mp_digit rr = a->dp[x] >> (mp_digit)(MP_DIGIT_BIT - 1);
+
+      /* now shift up this digit, add in the carry [from the previous] */
+      b->dp[x] = ((a->dp[x] << 1uL) | r) & MP_MASK;
+
+      /* copy the carry that would be from the source
+       * digit into the next iteration
+       */
+      r = rr;
+   }
+
+   /* new leading digit? */
+   if (r != 0u) {
+      /* add a MSB which is always 1 at this point */
+      b->dp[b->used++] = 1;
+   }
+
+   /* now zero any excess digits on the destination
+    * that we didn't write to
+    */
+   s_mp_zero_digs(b->dp + b->used, oldused - b->used);
+
+   b->sign = a->sign;
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_mul_2d.c

@@ -0,0 +1,63 @@
+#include "tommath_private.h"
+#ifdef MP_MUL_2D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* shift left by a certain bit count */
+mp_err mp_mul_2d(const mp_int *a, int b, mp_int *c)
+{
+   mp_err err;
+
+   if (b < 0) {
+      return MP_VAL;
+   }
+
+   if ((err = mp_copy(a, c)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = mp_grow(c, c->used + (b / MP_DIGIT_BIT) + 1)) != MP_OKAY) {
+      return err;
+   }
+
+   /* shift by as many digits in the bit count */
+   if (b >= MP_DIGIT_BIT) {
+      if ((err = mp_lshd(c, b / MP_DIGIT_BIT)) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   /* shift any bit count < MP_DIGIT_BIT */
+   b %= MP_DIGIT_BIT;
+   if (b != 0u) {
+      mp_digit shift, mask, r;
+      int x;
+
+      /* bitmask for carries */
+      mask = ((mp_digit)1 << b) - (mp_digit)1;
+
+      /* shift for msbs */
+      shift = (mp_digit)(MP_DIGIT_BIT - b);
+
+      /* carry */
+      r    = 0;
+      for (x = 0; x < c->used; x++) {
+         /* get the higher bits of the current word */
+         mp_digit rr = (c->dp[x] >> shift) & mask;
+
+         /* shift the current word and OR in the carry */
+         c->dp[x] = ((c->dp[x] << b) | r) & MP_MASK;
+
+         /* set the carry to the carry bits of the current word */
+         r = rr;
+      }
+
+      /* set final carry */
+      if (r != 0u) {
+         c->dp[(c->used)++] = r;
+      }
+   }
+   mp_clamp(c);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_mul_d.c

@@ -0,0 +1,68 @@
+#include "tommath_private.h"
+#ifdef MP_MUL_D_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* multiply by a digit */
+mp_err mp_mul_d(const mp_int *a, mp_digit b, mp_int *c)
+{
+   mp_digit u;
+   mp_err   err;
+   int   ix, oldused;
+
+   if (b == 1u) {
+      return mp_copy(a, c);
+   }
+
+   /* power of two ? */
+   if (MP_HAS(MP_MUL_2) && (b == 2u)) {
+      return mp_mul_2(a, c);
+   }
+   if (MP_HAS(MP_MUL_2D) && MP_IS_2EXPT(b)) {
+      ix = 1;
+      while ((ix < MP_DIGIT_BIT) && (b != (((mp_digit)1)<<ix))) {
+         ix++;
+      }
+      return mp_mul_2d(a, ix, c);
+   }
+
+   /* make sure c is big enough to hold a*b */
+   if ((err = mp_grow(c, a->used + 1)) != MP_OKAY) {
+      return err;
+   }
+
+   /* get the original destinations used count */
+   oldused = c->used;
+
+   /* set the sign */
+   c->sign = a->sign;
+
+   /* zero carry */
+   u = 0;
+
+   /* compute columns */
+   for (ix = 0; ix < a->used; ix++) {
+      /* compute product and carry sum for this term */
+      mp_word r       = (mp_word)u + ((mp_word)a->dp[ix] * (mp_word)b);
+
+      /* mask off higher bits to get a single digit */
+      c->dp[ix] = (mp_digit)(r & (mp_word)MP_MASK);
+
+      /* send carry into next iteration */
+      u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
+   }
+
+   /* store final carry [if any] and increment ix offset  */
+   c->dp[ix] = u;
+
+   /* set used count */
+   c->used = a->used + 1;
+
+   /* now zero digits above the top */
+   s_mp_zero_digs(c->dp + c->used, oldused - c->used);
+
+   mp_clamp(c);
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_mulmod.c

@@ -0,0 +1,15 @@
+#include "tommath_private.h"
+#ifdef MP_MULMOD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* d = a * b (mod c) */
+mp_err mp_mulmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d)
+{
+   mp_err err;
+   if ((err = mp_mul(a, b, d)) != MP_OKAY) {
+      return err;
+   }
+   return mp_mod(d, c, d);
+}
+#endif

src/libtommath/mp_neg.c

@@ -0,0 +1,18 @@
+#include "tommath_private.h"
+#ifdef MP_NEG_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* b = -a */
+mp_err mp_neg(const mp_int *a, mp_int *b)
+{
+   mp_err err;
+   if ((err = mp_copy(a, b)) != MP_OKAY) {
+      return err;
+   }
+
+   b->sign = ((!mp_iszero(b) && !mp_isneg(b)) ? MP_NEG : MP_ZPOS);
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_or.c

@@ -0,0 +1,54 @@
+#include "tommath_private.h"
+#ifdef MP_OR_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 used = MP_MAX(a->used, b->used) + 1, i;
+   mp_err err;
+   mp_digit ac = 1, bc = 1, cc = 1;
+   bool neg = (mp_isneg(a) || mp_isneg(b));
+
+   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 (mp_isneg(a)) {
+         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 (mp_isneg(b)) {
+         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 (neg) {
+         cc += ~c->dp[i] & MP_MASK;
+         c->dp[i] = cc & MP_MASK;
+         cc >>= MP_DIGIT_BIT;
+      }
+   }
+
+   c->used = used;
+   c->sign = (neg ? MP_NEG : MP_ZPOS);
+   mp_clamp(c);
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_pack.c

@@ -0,0 +1,69 @@
+#include "tommath_private.h"
+#ifdef MP_PACK_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* based on gmp's mpz_export.
+ * see http://gmplib.org/manual/Integer-Import-and-Export.html
+ */
+mp_err mp_pack(void *rop, size_t maxcount, size_t *written, mp_order order, size_t size,
+               mp_endian endian, size_t nails, const mp_int *op)
+{
+   mp_err err;
+   size_t odd_nails, nail_bytes, i, j, count;
+   uint8_t odd_nail_mask;
+
+   mp_int t;
+
+   count = mp_pack_count(op, nails, size);
+
+   if (count > maxcount) {
+      return MP_BUF;
+   }
+
+   if ((err = mp_init_copy(&t, op)) != MP_OKAY) {
+      return err;
+   }
+
+   if (endian == MP_NATIVE_ENDIAN) {
+      MP_GET_ENDIANNESS(endian);
+   }
+
+   odd_nails = (nails % 8u);
+   odd_nail_mask = 0xff;
+   for (i = 0u; i < odd_nails; ++i) {
+      odd_nail_mask ^= (uint8_t)(1u << (7u - i));
+   }
+   nail_bytes = nails / 8u;
+
+   for (i = 0u; i < count; ++i) {
+      for (j = 0u; j < size; ++j) {
+         uint8_t *byte = (uint8_t *)rop +
+                         (((order == MP_LSB_FIRST) ? i : ((count - 1u) - i)) * size) +
+                         ((endian == MP_LITTLE_ENDIAN) ? j : ((size - 1u) - j));
+
+         if (j >= (size - nail_bytes)) {
+            *byte = 0;
+            continue;
+         }
+
+         *byte = (uint8_t)((j == ((size - nail_bytes) - 1u)) ? (t.dp[0] & odd_nail_mask) : (t.dp[0] & 0xFFuL));
+
+         if ((err = mp_div_2d(&t, (j == ((size - nail_bytes) - 1u)) ? (int)(8u - odd_nails) : 8, &t, NULL)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+
+      }
+   }
+
+   if (written != NULL) {
+      *written = count;
+   }
+   err = MP_OKAY;
+
+LBL_ERR:
+   mp_clear(&t);
+   return err;
+}
+
+#endif

src/libtommath/mp_pack_count.c

@@ -0,0 +1,12 @@
+#include "tommath_private.h"
+#ifdef MP_PACK_COUNT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+size_t mp_pack_count(const mp_int *a, size_t nails, size_t size)
+{
+   size_t bits = (size_t)mp_count_bits(a);
+   return ((bits / ((size * 8u) - nails)) + (((bits % ((size * 8u) - nails)) != 0u) ? 1u : 0u));
+}
+
+#endif

src/libtommath/mp_prime_fermat.c

@@ -0,0 +1,41 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_FERMAT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* performs one Fermat test.
+ *
+ * If "a" were prime then b**a == b (mod a) since the order of
+ * the multiplicative sub-group would be phi(a) = a-1.  That means
+ * it would be the same as b**(a mod (a-1)) == b**1 == b (mod a).
+ *
+ * Sets result to 1 if the congruence holds, or zero otherwise.
+ */
+mp_err mp_prime_fermat(const mp_int *a, const mp_int *b, bool *result)
+{
+   mp_int  t;
+   mp_err  err;
+
+   /* ensure b > 1 */
+   if (mp_cmp_d(b, 1uL) != MP_GT) {
+      return MP_VAL;
+   }
+
+   /* init t */
+   if ((err = mp_init(&t)) != MP_OKAY) {
+      return err;
+   }
+
+   /* compute t = b**a mod a */
+   if ((err = mp_exptmod(b, a, a, &t)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   /* is it equal to b? */
+   *result = mp_cmp(&t, b) == MP_EQ;
+
+LBL_ERR:
+   mp_clear(&t);
+   return err;
+}
+#endif

src/libtommath/mp_prime_frobenius_underwood.c

@@ -0,0 +1,127 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_FROBENIUS_UNDERWOOD_C
+
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/*
+ *  See file mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
+ */
+#ifndef LTM_USE_ONLY_MR
+
+/*
+ * floor of positive solution of
+ * (2^16)-1 = (a+4)*(2*a+5)
+ * TODO: Both values are smaller than N^(1/4), would have to use a bigint
+ *       for a instead but any a biger than about 120 are already so rare that
+ *       it is possible to ignore them and still get enough pseudoprimes.
+ *       But it is still a restriction of the set of available pseudoprimes
+ *       which makes this implementation less secure if used stand-alone.
+ */
+#define LTM_FROBENIUS_UNDERWOOD_A 32764
+
+mp_err mp_prime_frobenius_underwood(const mp_int *N, bool *result)
+{
+   mp_int T1z, T2z, Np1z, sz, tz;
+   int a, ap2, i;
+   mp_err err;
+
+   if ((err = mp_init_multi(&T1z, &T2z, &Np1z, &sz, &tz, NULL)) != MP_OKAY) {
+      return err;
+   }
+
+   for (a = 0; a < LTM_FROBENIUS_UNDERWOOD_A; a++) {
+      int j;
+
+      /* TODO: That's ugly! No, really, it is! */
+      if ((a==2) || (a==4) || (a==7) || (a==8) || (a==10) ||
+          (a==14) || (a==18) || (a==23) || (a==26) || (a==28)) {
+         continue;
+      }
+
+      mp_set_i32(&T1z, (int32_t)((a * a) - 4));
+
+      if ((err = mp_kronecker(&T1z, N, &j)) != MP_OKAY)           goto LBL_END;
+
+      if (j == -1) {
+         break;
+      }
+
+      if (j == 0) {
+         /* composite */
+         *result = false;
+         goto LBL_END;
+      }
+   }
+   /* Tell it a composite and set return value accordingly */
+   if (a >= LTM_FROBENIUS_UNDERWOOD_A) {
+      err = MP_ITER;
+      goto LBL_END;
+   }
+   /* Composite if N and (a+4)*(2*a+5) are not coprime */
+   mp_set_u32(&T1z, (uint32_t)((a+4)*((2*a)+5)));
+
+   if ((err = mp_gcd(N, &T1z, &T1z)) != MP_OKAY)                  goto LBL_END;
+
+   if (!((T1z.used == 1) && (T1z.dp[0] == 1u))) {
+      /* composite */
+      *result = false;
+      goto LBL_END;
+   }
+
+   ap2 = a + 2;
+   if ((err = mp_add_d(N, 1uL, &Np1z)) != MP_OKAY)                goto LBL_END;
+
+   mp_set(&sz, 1uL);
+   mp_set(&tz, 2uL);
+
+   for (i = mp_count_bits(&Np1z) - 2; i >= 0; i--) {
+      /*
+       * temp = (sz*(a*sz+2*tz))%N;
+       * tz   = ((tz-sz)*(tz+sz))%N;
+       * sz   = temp;
+       */
+      if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY)                 goto LBL_END;
+
+      /* a = 0 at about 50% of the cases (non-square and odd input) */
+      if (a != 0) {
+         if ((err = mp_mul_d(&sz, (mp_digit)a, &T1z)) != MP_OKAY) goto LBL_END;
+         if ((err = mp_add(&T1z, &T2z, &T2z)) != MP_OKAY)         goto LBL_END;
+      }
+
+      if ((err = mp_mul(&T2z, &sz, &T1z)) != MP_OKAY)             goto LBL_END;
+      if ((err = mp_sub(&tz, &sz, &T2z)) != MP_OKAY)              goto LBL_END;
+      if ((err = mp_add(&sz, &tz, &sz)) != MP_OKAY)               goto LBL_END;
+      if ((err = mp_mul(&sz, &T2z, &tz)) != MP_OKAY)              goto LBL_END;
+      if ((err = mp_mod(&tz, N, &tz)) != MP_OKAY)                 goto LBL_END;
+      if ((err = mp_mod(&T1z, N, &sz)) != MP_OKAY)                goto LBL_END;
+      if (s_mp_get_bit(&Np1z, i)) {
+         /*
+          *  temp = (a+2) * sz + tz
+          *  tz   = 2 * tz - sz
+          *  sz   = temp
+          */
+         if (a == 0) {
+            if ((err = mp_mul_2(&sz, &T1z)) != MP_OKAY)           goto LBL_END;
+         } else {
+            if ((err = mp_mul_d(&sz, (mp_digit)ap2, &T1z)) != MP_OKAY) goto LBL_END;
+         }
+         if ((err = mp_add(&T1z, &tz, &T1z)) != MP_OKAY)          goto LBL_END;
+         if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY)              goto LBL_END;
+         if ((err = mp_sub(&T2z, &sz, &tz)) != MP_OKAY)           goto LBL_END;
+         mp_exch(&sz, &T1z);
+      }
+   }
+
+   mp_set_u32(&T1z, (uint32_t)((2 * a) + 5));
+   if ((err = mp_mod(&T1z, N, &T1z)) != MP_OKAY)                  goto LBL_END;
+
+   *result = mp_iszero(&sz) && (mp_cmp(&tz, &T1z) == MP_EQ);
+
+LBL_END:
+   mp_clear_multi(&tz, &sz, &Np1z, &T2z, &T1z, NULL);
+   return err;
+}
+
+#endif
+#endif

src/libtommath/mp_prime_is_prime.c

@@ -0,0 +1,282 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_IS_PRIME_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* portable integer log of two with small footprint */
+static unsigned int s_floor_ilog2(int value)
+{
+   unsigned int r = 0;
+   while ((value >>= 1) != 0) {
+      r++;
+   }
+   return r;
+}
+
+mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
+{
+   mp_int  b;
+   int     ix;
+   bool    res;
+   mp_err  err;
+
+   /* default to no */
+   *result = false;
+
+   /* Some shortcuts */
+   /* N > 3 */
+   if (a->used == 1) {
+      if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
+         *result = false;
+         return MP_OKAY;
+      }
+      if (a->dp[0] == 2u) {
+         *result = true;
+         return MP_OKAY;
+      }
+   }
+
+   /* N must be odd */
+   if (mp_iseven(a)) {
+      return MP_OKAY;
+   }
+   /* N is not a perfect square: floor(sqrt(N))^2 != N */
+   if ((err = mp_is_square(a, &res)) != MP_OKAY) {
+      return err;
+   }
+   if (res) {
+      return MP_OKAY;
+   }
+
+   /* is the input equal to one of the primes in the table? */
+   for (ix = 0; ix < MP_PRIME_TAB_SIZE; ix++) {
+      if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
+         *result = true;
+         return MP_OKAY;
+      }
+   }
+   /* first perform trial division */
+   if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
+      return err;
+   }
+
+   /* return if it was trivially divisible */
+   if (res) {
+      return MP_OKAY;
+   }
+
+   /*
+       Run the Miller-Rabin test with base 2 for the BPSW test.
+    */
+   if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+      goto LBL_B;
+   }
+   if (!res) {
+      goto LBL_B;
+   }
+   /*
+      Rumours have it that Mathematica does a second M-R test with base 3.
+      Other rumours have it that their strong L-S test is slightly different.
+      It does not hurt, though, beside a bit of extra runtime.
+   */
+   b.dp[0]++;
+   if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+      goto LBL_B;
+   }
+   if (!res) {
+      goto LBL_B;
+   }
+
+   /*
+    * Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
+    * slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
+    * bases 2, 3 and t random bases.
+    */
+#ifndef LTM_USE_ONLY_MR
+   if (t >= 0) {
+#ifdef LTM_USE_FROBENIUS_TEST
+      err = mp_prime_frobenius_underwood(a, &res);
+      if ((err != MP_OKAY) && (err != MP_ITER)) {
+         goto LBL_B;
+      }
+      if (!res) {
+         goto LBL_B;
+      }
+#else
+      if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
+         goto LBL_B;
+      }
+      if (!res) {
+         goto LBL_B;
+      }
+#endif
+   }
+#endif
+
+   /* run at least one Miller-Rabin test with a random base */
+   if (t == 0) {
+      t = 1;
+   }
+
+   /*
+      Only recommended if the input range is known to be < 3317044064679887385961981
+
+      It uses the bases necessary for a deterministic M-R test if the input is
+      smaller than  3317044064679887385961981
+      The caller has to check the size.
+      TODO: can be made a bit finer grained but comparing is not free.
+   */
+   if (t < 0) {
+      int p_max = 0;
+
+      /*
+          Sorenson, Jonathan; Webster, Jonathan (2015).
+           "Strong Pseudoprimes to Twelve Prime Bases".
+       */
+      /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
+      if ((err =   mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
+         goto LBL_B;
+      }
+
+      if (mp_cmp(a, &b) == MP_LT) {
+         p_max = 12;
+      } else {
+         /* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
+         if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
+            goto LBL_B;
+         }
+
+         if (mp_cmp(a, &b) == MP_LT) {
+            p_max = 13;
+         } else {
+            err = MP_VAL;
+            goto LBL_B;
+         }
+      }
+
+      /* we did bases 2 and 3  already, skip them */
+      for (ix = 2; ix < p_max; ix++) {
+         mp_set(&b, s_mp_prime_tab[ix]);
+         if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+            goto LBL_B;
+         }
+         if (!res) {
+            goto LBL_B;
+         }
+      }
+   }
+   /*
+       Do "t" M-R tests with random bases between 3 and "a".
+       See Fips 186.4 p. 126ff
+   */
+   else if (t > 0) {
+      unsigned int mask;
+      int size_a;
+
+      /*
+       * The mp_digit's have a defined bit-size but the size of the
+       * array a.dp is a simple 'int' and this library can not assume full
+       * compliance to the current C-standard (ISO/IEC 9899:2011) because
+       * it gets used for small embeded processors, too. Some of those MCUs
+       * have compilers that one cannot call standard compliant by any means.
+       * Hence the ugly type-fiddling in the following code.
+       */
+      size_a = mp_count_bits(a);
+      mask = (1u << s_floor_ilog2(size_a)) - 1u;
+      /*
+         Assuming the General Rieman hypothesis (never thought to write that in a
+         comment) the upper bound can be lowered to  2*(log a)^2.
+         E. Bach, "Explicit bounds for primality testing and related problems,"
+         Math. Comp. 55 (1990), 355-380.
+
+            size_a = (size_a/10) * 7;
+            len = 2 * (size_a * size_a);
+
+         E.g.: a number of size 2^2048 would be reduced to the upper limit
+
+            floor(2048/10)*7 = 1428
+            2 * 1428^2       = 4078368
+
+         (would have been ~4030331.9962 with floats and natural log instead)
+         That number is smaller than 2^28, the default bit-size of mp_digit.
+      */
+
+      /*
+        How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
+        does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
+        Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.
+
+        The function mp_rand() goes to some length to use a cryptographically
+        good PRNG. That also means that the chance to always get the same base
+        in the loop is non-zero, although very low.
+        If the BPSW test and/or the addtional Frobenious test have been
+        performed instead of just the Miller-Rabin test with the bases 2 and 3,
+        a single extra test should suffice, so such a very unlikely event
+        will not do much harm.
+
+        To preemptivly answer the dangling question: no, a witness does not
+        need to be prime.
+      */
+      for (ix = 0; ix < t; ix++) {
+         unsigned int fips_rand;
+         int len;
+
+         /* mp_rand() guarantees the first digit to be non-zero */
+         if ((err = mp_rand(&b, 1)) != MP_OKAY) {
+            goto LBL_B;
+         }
+         /*
+          * Reduce digit before casting because mp_digit might be bigger than
+          * an unsigned int and "mask" on the other side is most probably not.
+          */
+         fips_rand = (unsigned int)(b.dp[0] & (mp_digit) mask);
+         if (fips_rand > (unsigned int)(INT_MAX - MP_DIGIT_BIT)) {
+            len = INT_MAX / MP_DIGIT_BIT;
+         } else {
+            len = (((int)fips_rand + MP_DIGIT_BIT) / MP_DIGIT_BIT);
+         }
+         /*  Unlikely. */
+         if (len < 0) {
+            ix--;
+            continue;
+         }
+         if ((err = mp_rand(&b, len)) != MP_OKAY) {
+            goto LBL_B;
+         }
+         /*
+          * That number might got too big and the witness has to be
+          * smaller than "a"
+          */
+         len = mp_count_bits(&b);
+         if (len >= size_a) {
+            len = (len - size_a) + 1;
+            if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
+               goto LBL_B;
+            }
+         }
+         /* Although the chance for b <= 3 is miniscule, try again. */
+         if (mp_cmp_d(&b, 3uL) != MP_GT) {
+            ix--;
+            continue;
+         }
+         if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+            goto LBL_B;
+         }
+         if (!res) {
+            goto LBL_B;
+         }
+      }
+   }
+
+   /* passed the test */
+   *result = true;
+LBL_B:
+   mp_clear(&b);
+   return err;
+}
+
+#endif

src/libtommath/mp_prime_miller_rabin.c

@@ -0,0 +1,91 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_MILLER_RABIN_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Miller-Rabin test of "a" to the base of "b" as described in
+ * HAC pp. 139 Algorithm 4.24
+ *
+ * Sets result to 0 if definitely composite or 1 if probably prime.
+ * Randomly the chance of error is no more than 1/4 and often
+ * very much lower.
+ */
+mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
+{
+   mp_int  n1, y, r;
+   mp_err  err;
+   int     s, j;
+
+   /* ensure b > 1 */
+   if (mp_cmp_d(b, 1uL) != MP_GT) {
+      return MP_VAL;
+   }
+
+   /* get n1 = a - 1 */
+   if ((err = mp_init_copy(&n1, a)) != MP_OKAY) {
+      return err;
+   }
+   if ((err = mp_sub_d(&n1, 1uL, &n1)) != MP_OKAY) {
+      goto LBL_ERR1;
+   }
+
+   /* set 2**s * r = n1 */
+   if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
+      goto LBL_ERR1;
+   }
+
+   /* count the number of least significant bits
+    * which are zero
+    */
+   s = mp_cnt_lsb(&r);
+
+   /* now divide n - 1 by 2**s */
+   if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) {
+      goto LBL_ERR2;
+   }
+
+   /* compute y = b**r mod a */
+   if ((err = mp_init(&y)) != MP_OKAY) {
+      goto LBL_ERR2;
+   }
+   if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
+      goto LBL_END;
+   }
+
+   /* if y != 1 and y != n1 do */
+   if ((mp_cmp_d(&y, 1uL) != MP_EQ) && (mp_cmp(&y, &n1) != MP_EQ)) {
+      j = 1;
+      /* while j <= s-1 and y != n1 */
+      while ((j <= (s - 1)) && (mp_cmp(&y, &n1) != MP_EQ)) {
+         if ((err = mp_sqrmod(&y, a, &y)) != MP_OKAY) {
+            goto LBL_END;
+         }
+
+         /* if y == 1 then composite */
+         if (mp_cmp_d(&y, 1uL) == MP_EQ) {
+            *result = false;
+            goto LBL_END;
+         }
+
+         ++j;
+      }
+
+      /* if y != n1 then composite */
+      if (mp_cmp(&y, &n1) != MP_EQ) {
+         *result = false;
+         goto LBL_END;
+      }
+   }
+
+   /* probably prime now */
+   *result = true;
+
+LBL_END:
+   mp_clear(&y);
+LBL_ERR2:
+   mp_clear(&r);
+LBL_ERR1:
+   mp_clear(&n1);
+   return err;
+}
+#endif

src/libtommath/mp_prime_next_prime.c

@@ -0,0 +1,127 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_NEXT_PRIME_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* finds the next prime after the number "a" using "t" trials
+ * of Miller-Rabin.
+ *
+ * bbs_style = true means the prime must be congruent to 3 mod 4
+ */
+mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
+{
+   int      x;
+   mp_err   err;
+   bool  res = false;
+   mp_digit res_tab[MP_PRIME_TAB_SIZE], kstep;
+   mp_int   b;
+
+   /* force positive */
+   a->sign = MP_ZPOS;
+
+   /* simple algo if a is less than the largest prime in the table */
+   if (mp_cmp_d(a, s_mp_prime_tab[MP_PRIME_TAB_SIZE-1]) == MP_LT) {
+      /* find which prime it is bigger than "a" */
+      for (x = 0; x < MP_PRIME_TAB_SIZE; x++) {
+         mp_ord cmp = mp_cmp_d(a, s_mp_prime_tab[x]);
+         if (cmp == MP_EQ) {
+            continue;
+         }
+         if (cmp != MP_GT) {
+            if ((bbs_style) && ((s_mp_prime_tab[x] & 3u) != 3u)) {
+               /* try again until we get a prime congruent to 3 mod 4 */
+               continue;
+            } else {
+               mp_set(a, s_mp_prime_tab[x]);
+               return MP_OKAY;
+            }
+         }
+      }
+      /* fall through to the sieve */
+   }
+
+   /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
+   kstep = bbs_style ? 4 : 2;
+
+   /* at this point we will use a combination of a sieve and Miller-Rabin */
+
+   if (bbs_style) {
+      /* if a mod 4 != 3 subtract the correct value to make it so */
+      if ((a->dp[0] & 3u) != 3u) {
+         if ((err = mp_sub_d(a, (a->dp[0] & 3u) + 1u, a)) != MP_OKAY) {
+            return err;
+         }
+      }
+   } else {
+      if (mp_iseven(a)) {
+         /* force odd */
+         if ((err = mp_sub_d(a, 1uL, a)) != MP_OKAY) {
+            return err;
+         }
+      }
+   }
+
+   /* generate the restable */
+   for (x = 1; x < MP_PRIME_TAB_SIZE; x++) {
+      if ((err = mp_mod_d(a, s_mp_prime_tab[x], res_tab + x)) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   /* init temp used for Miller-Rabin Testing */
+   if ((err = mp_init(&b)) != MP_OKAY) {
+      return err;
+   }
+
+   for (;;) {
+      mp_digit step = 0;
+      bool y;
+      /* skip to the next non-trivially divisible candidate */
+      do {
+         /* y == true if any residue was zero [e.g. cannot be prime] */
+         y     = false;
+
+         /* increase step to next candidate */
+         step += kstep;
+
+         /* compute the new residue without using division */
+         for (x = 1; x < MP_PRIME_TAB_SIZE; x++) {
+            /* add the step to each residue */
+            res_tab[x] += kstep;
+
+            /* subtract the modulus [instead of using division] */
+            if (res_tab[x] >= s_mp_prime_tab[x]) {
+               res_tab[x]  -= s_mp_prime_tab[x];
+            }
+
+            /* set flag if zero */
+            if (res_tab[x] == 0u) {
+               y = true;
+            }
+         }
+      } while (y && (step < (((mp_digit)1 << MP_DIGIT_BIT) - kstep)));
+
+      /* add the step */
+      if ((err = mp_add_d(a, step, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* if didn't pass sieve and step == MP_MAX then skip test */
+      if (y && (step >= (((mp_digit)1 << MP_DIGIT_BIT) - kstep))) {
+         continue;
+      }
+
+      if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (res) {
+         break;
+      }
+   }
+
+LBL_ERR:
+   mp_clear(&b);
+   return err;
+}
+
+#endif

src/libtommath/mp_prime_rabin_miller_trials.c

@@ -0,0 +1,48 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_RABIN_MILLER_TRIALS_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+static const struct {
+   int k, t;
+} sizes[] = {
+   {    80, -1 }, /* Use deterministic algorithm for size <= 80 bits */
+   {    81, 37 }, /* max. error = 2^(-96)*/
+   {    96, 32 }, /* max. error = 2^(-96)*/
+   {   128, 40 }, /* max. error = 2^(-112)*/
+   {   160, 35 }, /* max. error = 2^(-112)*/
+   {   256, 27 }, /* max. error = 2^(-128)*/
+   {   384, 16 }, /* max. error = 2^(-128)*/
+   {   512, 18 }, /* max. error = 2^(-160)*/
+   {   768, 11 }, /* max. error = 2^(-160)*/
+   {   896, 10 }, /* max. error = 2^(-160)*/
+   {  1024, 12 }, /* max. error = 2^(-192)*/
+   {  1536, 8  }, /* max. error = 2^(-192)*/
+   {  2048, 6  }, /* max. error = 2^(-192)*/
+   {  3072, 4  }, /* max. error = 2^(-192)*/
+   {  4096, 5  }, /* max. error = 2^(-256)*/
+   {  5120, 4  }, /* max. error = 2^(-256)*/
+   {  6144, 4  }, /* max. error = 2^(-256)*/
+   {  8192, 3  }, /* max. error = 2^(-256)*/
+   {  9216, 3  }, /* max. error = 2^(-256)*/
+   { 10240, 2  }  /* For bigger keysizes use always at least 2 Rounds */
+};
+
+/* returns # of RM trials required for a given bit size */
+int mp_prime_rabin_miller_trials(int size)
+{
+   int x;
+
+   for (x = 0; x < (int)(sizeof(sizes)/(sizeof(sizes[0]))); x++) {
+      if (sizes[x].k == size) {
+         return sizes[x].t;
+      }
+      if (sizes[x].k > size) {
+         return (x == 0) ? sizes[0].t : sizes[x - 1].t;
+      }
+   }
+   return sizes[x-1].t;
+}
+
+
+#endif

src/libtommath/mp_prime_rand.c

@@ -0,0 +1,123 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_RAND_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* makes a truly random prime of a given size (bits),
+ *
+ * Flags are as follows:
+ *
+ *   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
+ *
+ */
+
+/* This is possibly the mother of all prime generation functions, muahahahahaha! */
+mp_err mp_prime_rand(mp_int *a, int t, int size, int flags)
+{
+   uint8_t *tmp, maskAND, maskOR_msb, maskOR_lsb;
+   int bsize, maskOR_msb_offset;
+   bool res;
+   mp_err err;
+
+   /* sanity check the input */
+   if ((size <= 1) || (t <= 0)) {
+      return MP_VAL;
+   }
+
+   /* MP_PRIME_SAFE implies MP_PRIME_BBS */
+   if ((flags & MP_PRIME_SAFE) != 0) {
+      flags |= MP_PRIME_BBS;
+   }
+
+   /* calc the byte size */
+   bsize = (size>>3) + ((size&7)?1:0);
+
+   /* we need a buffer of bsize bytes */
+   tmp = (uint8_t *) MP_MALLOC((size_t)bsize);
+   if (tmp == NULL) {
+      return MP_MEM;
+   }
+
+   /* calc the maskAND value for the MSbyte*/
+   maskAND = ((size&7) == 0) ? 0xFFu : (uint8_t)(0xFFu >> (8 - (size & 7)));
+
+   /* calc the maskOR_msb */
+   maskOR_msb        = 0;
+   maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
+   if ((flags & MP_PRIME_2MSB_ON) != 0) {
+      maskOR_msb       |= (uint8_t)(0x80 >> ((9 - size) & 7));
+   }
+
+   /* get the maskOR_lsb */
+   maskOR_lsb         = 1u;
+   if ((flags & MP_PRIME_BBS) != 0) {
+      maskOR_lsb     |= 3u;
+   }
+
+   do {
+      /* read the bytes */
+      if ((err = s_mp_rand_source(tmp, (size_t)bsize)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* work over the MSbyte */
+      tmp[0]    &= maskAND;
+      tmp[0]    |= (uint8_t)(1 << ((size - 1) & 7));
+
+      /* mix in the maskORs */
+      tmp[maskOR_msb_offset]   |= maskOR_msb;
+      tmp[bsize-1]             |= maskOR_lsb;
+
+      /* read it in */
+      /* TODO: casting only for now until all lengths have been changed to the type "size_t"*/
+      if ((err = mp_from_ubin(a, tmp, (size_t)bsize)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* is it prime? */
+      if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (!res) {
+         continue;
+      }
+
+      if ((flags & MP_PRIME_SAFE) != 0) {
+         /* see if (a-1)/2 is prime */
+         if ((err = mp_sub_d(a, 1uL, a)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_div_2(a, a)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+
+         /* is it prime? */
+         if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+      }
+   } while (!res);
+
+   if ((flags & MP_PRIME_SAFE) != 0) {
+      /* restore a to the original value */
+      if ((err = mp_mul_2(a, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_add_d(a, 1uL, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+
+   err = MP_OKAY;
+LBL_ERR:
+   MP_FREE_BUF(tmp, (size_t)bsize);
+   return err;
+}
+
+#endif

src/libtommath/mp_prime_strong_lucas_selfridge.c

@@ -0,0 +1,281 @@
+#include "tommath_private.h"
+#ifdef MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
+
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/*
+ *  See file mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
+ */
+#ifndef LTM_USE_ONLY_MR
+
+/*
+ * multiply bigint a with int d and put the result in c
+ * Like mp_mul_d() but with a signed long as the small input
+ */
+static mp_err s_mul_si(const mp_int *a, int32_t d, mp_int *c)
+{
+   mp_int t;
+   mp_err err;
+
+   if ((err = mp_init(&t)) != MP_OKAY) {
+      return err;
+   }
+
+   /*
+    * mp_digit might be smaller than a long, which excludes
+    * the use of mp_mul_d() here.
+    */
+   mp_set_i32(&t, d);
+   err = mp_mul(a, &t, c);
+   mp_clear(&t);
+   return err;
+}
+/*
+    Strong Lucas-Selfridge test.
+    returns true if it is a strong L-S prime, false if it is composite
+
+    Code ported from  Thomas Ray Nicely's implementation of the BPSW test
+    at http://www.trnicely.net/misc/bpsw.html
+
+    Freeware copyright (C) 2016 Thomas R. Nicely <http://www.trnicely.net>.
+    Released into the public domain by the author, who disclaims any legal
+    liability arising from its use
+
+    The multi-line comments are made by Thomas R. Nicely and are copied verbatim.
+    Additional comments marked "CZ" (without the quotes) are by the code-portist.
+
+    (If that name sounds familiar, he is the guy who found the fdiv bug in the
+     Pentium (P5x, I think) Intel processor)
+*/
+mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, bool *result)
+{
+   /* CZ TODO: choose better variable names! */
+   mp_int Dz, gcd, Np1, Uz, Vz, U2mz, V2mz, Qmz, Q2mz, Qkdz, T1z, T2z, T3z, T4z, Q2kdz;
+   int32_t D, Ds, J, sign, P, Q, r, s, u, Nbits;
+   mp_err err;
+   bool oddness;
+
+   *result = false;
+   /*
+   Find the first element D in the sequence {5, -7, 9, -11, 13, ...}
+   such that Jacobi(D,N) = -1 (Selfridge's algorithm). Theory
+   indicates that, if N is not a perfect square, D will "nearly
+   always" be "small." Just in case, an overflow trap for D is
+   included.
+   */
+
+   if ((err = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz,
+                            NULL)) != MP_OKAY) {
+      return err;
+   }
+
+   D = 5;
+   sign = 1;
+
+   for (;;) {
+      Ds   = sign * D;
+      sign = -sign;
+      mp_set_u32(&Dz, (uint32_t)D);
+      if ((err = mp_gcd(a, &Dz, &gcd)) != MP_OKAY)                goto LBL_LS_ERR;
+
+      /* if 1 < GCD < N then N is composite with factor "D", and
+         Jacobi(D,N) is technically undefined (but often returned
+         as zero). */
+      if ((mp_cmp_d(&gcd, 1uL) == MP_GT) && (mp_cmp(&gcd, a) == MP_LT)) {
+         goto LBL_LS_ERR;
+      }
+      if (Ds < 0) {
+         Dz.sign = MP_NEG;
+      }
+      if ((err = mp_kronecker(&Dz, a, &J)) != MP_OKAY)            goto LBL_LS_ERR;
+
+      if (J == -1) {
+         break;
+      }
+      D += 2;
+
+      if (D > (INT_MAX - 2)) {
+         err = MP_VAL;
+         goto LBL_LS_ERR;
+      }
+   }
+
+
+
+   P = 1;              /* Selfridge's choice */
+   Q = (1 - Ds) / 4;   /* Required so D = P*P - 4*Q */
+
+   /* NOTE: The conditions (a) N does not divide Q, and
+      (b) D is square-free or not a perfect square, are included by
+      some authors; e.g., "Prime numbers and computer methods for
+      factorization," Hans Riesel (2nd ed., 1994, Birkhauser, Boston),
+      p. 130. For this particular application of Lucas sequences,
+      these conditions were found to be immaterial. */
+
+   /* Now calculate N - Jacobi(D,N) = N + 1 (even), and calculate the
+      odd positive integer d and positive integer s for which
+      N + 1 = 2^s*d (similar to the step for N - 1 in Miller's test).
+      The strong Lucas-Selfridge test then returns N as a strong
+      Lucas probable prime (slprp) if any of the following
+      conditions is met: U_d=0, V_d=0, V_2d=0, V_4d=0, V_8d=0,
+      V_16d=0, ..., etc., ending with V_{2^(s-1)*d}=V_{(N+1)/2}=0
+      (all equalities mod N). Thus d is the highest index of U that
+      must be computed (since V_2m is independent of U), compared
+      to U_{N+1} for the standard Lucas-Selfridge test; and no
+      index of V beyond (N+1)/2 is required, just as in the
+      standard Lucas-Selfridge test. However, the quantity Q^d must
+      be computed for use (if necessary) in the latter stages of
+      the test. The result is that the strong Lucas-Selfridge test
+      has a running time only slightly greater (order of 10 %) than
+      that of the standard Lucas-Selfridge test, while producing
+      only (roughly) 30 % as many pseudoprimes (and every strong
+      Lucas pseudoprime is also a standard Lucas pseudoprime). Thus
+      the evidence indicates that the strong Lucas-Selfridge test is
+      more effective than the standard Lucas-Selfridge test, and a
+      Baillie-PSW test based on the strong Lucas-Selfridge test
+      should be more reliable. */
+
+   if ((err = mp_add_d(a, 1uL, &Np1)) != MP_OKAY)                 goto LBL_LS_ERR;
+   s = mp_cnt_lsb(&Np1);
+
+   /* CZ
+    * This should round towards zero because
+    * Thomas R. Nicely used GMP's mpz_tdiv_q_2exp()
+    * and mp_div_2d() is equivalent. Additionally:
+    * dividing an even number by two does not produce
+    * any leftovers.
+    */
+   if ((err = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY)          goto LBL_LS_ERR;
+   /* We must now compute U_d and V_d. Since d is odd, the accumulated
+      values U and V are initialized to U_1 and V_1 (if the target
+      index were even, U and V would be initialized instead to U_0=0
+      and V_0=2). The values of U_2m and V_2m are also initialized to
+      U_1 and V_1; the FOR loop calculates in succession U_2 and V_2,
+      U_4 and V_4, U_8 and V_8, etc. If the corresponding bits
+      (1, 2, 3, ...) of t are on (the zero bit having been accounted
+      for in the initialization of U and V), these values are then
+      combined with the previous totals for U and V, using the
+      composition formulas for addition of indices. */
+
+   mp_set(&Uz, 1uL);    /* U=U_1 */
+   mp_set(&Vz, (mp_digit)P);    /* V=V_1 */
+   mp_set(&U2mz, 1uL);  /* U_1 */
+   mp_set(&V2mz, (mp_digit)P);  /* V_1 */
+
+   mp_set_i32(&Qmz, Q);
+   if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY)                  goto LBL_LS_ERR;
+   /* Initializes calculation of Q^d */
+   mp_set_i32(&Qkdz, Q);
+
+   Nbits = mp_count_bits(&Dz);
+
+   for (u = 1; u < Nbits; u++) { /* zero bit off, already accounted for */
+      /* Formulas for doubling of indices (carried out mod N). Note that
+       * the indices denoted as "2m" are actually powers of 2, specifically
+       * 2^(ul-1) beginning each loop and 2^ul ending each loop.
+       *
+       * U_2m = U_m*V_m
+       * V_2m = V_m*V_m - 2*Q^m
+       */
+
+      if ((err = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY)         goto LBL_LS_ERR;
+      if ((err = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY)             goto LBL_LS_ERR;
+      if ((err = mp_sqr(&V2mz, &V2mz)) != MP_OKAY)                goto LBL_LS_ERR;
+      if ((err = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY)         goto LBL_LS_ERR;
+      if ((err = mp_mod(&V2mz, a, &V2mz)) != MP_OKAY)             goto LBL_LS_ERR;
+
+      /* Must calculate powers of Q for use in V_2m, also for Q^d later */
+      if ((err = mp_sqr(&Qmz, &Qmz)) != MP_OKAY)                  goto LBL_LS_ERR;
+
+      /* prevents overflow */ /* CZ  still necessary without a fixed prealloc'd mem.? */
+      if ((err = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY)               goto LBL_LS_ERR;
+      if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY)               goto LBL_LS_ERR;
+
+      if (s_mp_get_bit(&Dz, u)) {
+         /* Formulas for addition of indices (carried out mod N);
+          *
+          * U_(m+n) = (U_m*V_n + U_n*V_m)/2
+          * V_(m+n) = (V_m*V_n + D*U_m*U_n)/2
+          *
+          * Be careful with division by 2 (mod N)!
+          */
+         if ((err = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY)         goto LBL_LS_ERR;
+         if ((err = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY)         goto LBL_LS_ERR;
+         if ((err = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY)         goto LBL_LS_ERR;
+         if ((err = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY)         goto LBL_LS_ERR;
+         if ((err = s_mul_si(&T4z, Ds, &T4z)) != MP_OKAY)      goto LBL_LS_ERR;
+         if ((err = mp_add(&T1z, &T2z, &Uz)) != MP_OKAY)          goto LBL_LS_ERR;
+         if (mp_isodd(&Uz)) {
+            if ((err = mp_add(&Uz, a, &Uz)) != MP_OKAY)           goto LBL_LS_ERR;
+         }
+         /* CZ
+          * This should round towards negative infinity because
+          * Thomas R. Nicely used GMP's mpz_fdiv_q_2exp().
+          * But mp_div_2() does not do so, it is truncating instead.
+          */
+         oddness = mp_isodd(&Uz);
+         if ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY)               goto LBL_LS_ERR;
+         if (mp_isneg(&Uz) && oddness) {
+            if ((err = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY)       goto LBL_LS_ERR;
+         }
+         if ((err = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY)          goto LBL_LS_ERR;
+         if (mp_isodd(&Vz)) {
+            if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY)           goto LBL_LS_ERR;
+         }
+         oddness = mp_isodd(&Vz);
+         if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY)               goto LBL_LS_ERR;
+         if (mp_isneg(&Vz) && oddness) {
+            if ((err = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY)       goto LBL_LS_ERR;
+         }
+         if ((err = mp_mod(&Uz, a, &Uz)) != MP_OKAY)              goto LBL_LS_ERR;
+         if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY)              goto LBL_LS_ERR;
+
+         /* Calculating Q^d for later use */
+         if ((err = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY)       goto LBL_LS_ERR;
+         if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY)          goto LBL_LS_ERR;
+      }
+   }
+
+   /* If U_d or V_d is congruent to 0 mod N, then N is a prime or a
+      strong Lucas pseudoprime. */
+   if (mp_iszero(&Uz) || mp_iszero(&Vz)) {
+      *result = true;
+      goto LBL_LS_ERR;
+   }
+
+   /* NOTE: Ribenboim ("The new book of prime number records," 3rd ed.,
+      1995/6) omits the condition V0 on p.142, but includes it on
+      p. 130. The condition is NECESSARY; otherwise the test will
+      return false negatives---e.g., the primes 29 and 2000029 will be
+      returned as composite. */
+
+   /* Otherwise, we must compute V_2d, V_4d, V_8d, ..., V_{2^(s-1)*d}
+      by repeated use of the formula V_2m = V_m*V_m - 2*Q^m. If any of
+      these are congruent to 0 mod N, then N is a prime or a strong
+      Lucas pseudoprime. */
+
+   /* Initialize 2*Q^(d*2^r) for V_2m */
+   if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY)                goto LBL_LS_ERR;
+
+   for (r = 1; r < s; r++) {
+      if ((err = mp_sqr(&Vz, &Vz)) != MP_OKAY)                    goto LBL_LS_ERR;
+      if ((err = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY)            goto LBL_LS_ERR;
+      if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY)                 goto LBL_LS_ERR;
+      if (mp_iszero(&Vz)) {
+         *result = true;
+         goto LBL_LS_ERR;
+      }
+      /* Calculate Q^{d*2^r} for next r (final iteration irrelevant). */
+      if (r < (s - 1)) {
+         if ((err = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY)             goto LBL_LS_ERR;
+         if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY)          goto LBL_LS_ERR;
+         if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY)          goto LBL_LS_ERR;
+      }
+   }
+LBL_LS_ERR:
+   mp_clear_multi(&Q2kdz, &T4z, &T3z, &T2z, &T1z, &Qkdz, &Q2mz, &Qmz, &V2mz, &U2mz, &Vz, &Uz, &Np1, &gcd, &Dz, NULL);
+   return err;
+}
+#endif
+#endif

src/libtommath/mp_radix_size.c

@@ -0,0 +1,34 @@
+#include "tommath_private.h"
+#ifdef MP_RADIX_SIZE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* returns size of ASCII representation */
+mp_err mp_radix_size(const mp_int *a, int radix, size_t *size)
+{
+   mp_err err;
+   mp_int a_;
+   int b;
+
+   /* make sure the radix is in range */
+   if ((radix < 2) || (radix > 64)) {
+      return MP_VAL;
+   }
+
+   if (mp_iszero(a)) {
+      *size = 2;
+      return MP_OKAY;
+   }
+
+   a_ = *a;
+   a_.sign = MP_ZPOS;
+   if ((err = mp_log_n(&a_, radix, &b)) != MP_OKAY) {
+      return err;
+   }
+
+   /* mp_ilogb truncates to zero, hence we need one extra put on top and one for `\0`. */
+   *size = (size_t)b + 2U + (mp_isneg(a) ? 1U : 0U);
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_radix_size_overestimate.c

@@ -0,0 +1,17 @@
+#include "tommath_private.h"
+#ifdef MP_RADIX_SIZE_OVERESTIMATE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+mp_err mp_radix_size_overestimate(const mp_int *a, const int radix, size_t *size)
+{
+   if (MP_HAS(S_MP_RADIX_SIZE_OVERESTIMATE)) {
+      return s_mp_radix_size_overestimate(a, radix, size);
+   }
+   if (MP_HAS(MP_RADIX_SIZE)) {
+      return mp_radix_size(a, radix, size);
+   }
+   return MP_ERR;
+}
+
+#endif

src/libtommath/mp_rand.c

@@ -0,0 +1,39 @@
+#include "tommath_private.h"
+#ifdef MP_RAND_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+mp_err mp_rand(mp_int *a, int digits)
+{
+   int i;
+   mp_err err;
+
+   mp_zero(a);
+
+   if (digits <= 0) {
+      return MP_OKAY;
+   }
+
+   if ((err = mp_grow(a, digits)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = s_mp_rand_source(a->dp, (size_t)digits * sizeof(mp_digit))) != MP_OKAY) {
+      return err;
+   }
+
+   /* TODO: We ensure that the highest digit is nonzero. Should this be removed? */
+   while ((a->dp[digits - 1] & MP_MASK) == 0u) {
+      if ((err = s_mp_rand_source(a->dp + digits - 1, sizeof(mp_digit))) != MP_OKAY) {
+         return err;
+      }
+   }
+
+   a->used = digits;
+   for (i = 0; i < digits; ++i) {
+      a->dp[i] &= MP_MASK;
+   }
+
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_rand_source.c

@@ -0,0 +1,12 @@
+#include "tommath_private.h"
+#ifdef MP_RAND_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+mp_err(*s_mp_rand_source)(void *out, size_t size) = s_mp_rand_platform;
+
+void mp_rand_source(mp_err(*source)(void *out, size_t size))
+{
+   s_mp_rand_source = (source == NULL) ? s_mp_rand_platform : source;
+}
+#endif

src/libtommath/mp_read_radix.c

@@ -0,0 +1,69 @@
+#include "tommath_private.h"
+#ifdef MP_READ_RADIX_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* read a string [ASCII] in a given radix */
+mp_err mp_read_radix(mp_int *a, const char *str, int radix)
+{
+   mp_err   err;
+   mp_sign  sign = MP_ZPOS;
+
+   /* make sure the radix is ok */
+   if ((radix < 2) || (radix > 64)) {
+      return MP_VAL;
+   }
+
+   /* if the leading digit is a
+    * minus set the sign to negative.
+    */
+   if (*str == '-') {
+      ++str;
+      sign = MP_NEG;
+   }
+
+   /* set the integer to the default of zero */
+   mp_zero(a);
+
+   /* process each digit of the string */
+   while (*str != '\0') {
+      /* if the radix <= 36 the conversion is case insensitive
+       * this allows numbers like 1AB and 1ab to represent the same  value
+       * [e.g. in hex]
+       */
+      uint8_t y;
+      char ch = (radix <= 36) ? (char)MP_TOUPPER((int)*str) : *str;
+      unsigned pos = (unsigned)(ch - '+');
+      if (MP_RADIX_MAP_REVERSE_SIZE <= pos) {
+         break;
+      }
+      y = s_mp_radix_map_reverse[pos];
+
+      /* if the char was found in the map
+       * and is less than the given radix add it
+       * to the number, otherwise exit the loop.
+       */
+      if (y >= radix) {
+         break;
+      }
+      if ((err = mp_mul_d(a, (mp_digit)radix, a)) != MP_OKAY) {
+         return err;
+      }
+      if ((err = mp_add_d(a, y, a)) != MP_OKAY) {
+         return err;
+      }
+      ++str;
+   }
+
+   /* if an illegal character was found, fail. */
+   if ((*str != '\0') && (*str != '\r') && (*str != '\n')) {
+      return MP_VAL;
+   }
+
+   /* set the sign only if a != 0 */
+   if (!mp_iszero(a)) {
+      a->sign = sign;
+   }
+   return MP_OKAY;
+}
+#endif

src/libtommath/mp_reduce.c

@@ -0,0 +1,83 @@
+#include "tommath_private.h"
+#ifdef MP_REDUCE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* reduces x mod m, assumes 0 < x < m**2, mu is
+ * precomputed via mp_reduce_setup.
+ * From HAC pp.604 Algorithm 14.42
+ */
+mp_err mp_reduce(mp_int *x, const mp_int *m, const mp_int *mu)
+{
+   mp_int  q;
+   mp_err  err;
+   int     um = m->used;
+
+   /* q = x */
+   if ((err = mp_init_copy(&q, x)) != MP_OKAY) {
+      return err;
+   }
+
+   /* q1 = x / b**(k-1)  */
+   mp_rshd(&q, um - 1);
+
+   /* according to HAC this optimization is ok */
+   if ((mp_digit)um > ((mp_digit)1 << (MP_DIGIT_BIT - 1))) {
+      if ((err = mp_mul(&q, mu, &q)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   } else if (MP_HAS(S_MP_MUL_HIGH)) {
+      if ((err = s_mp_mul_high(&q, mu, &q, um)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   } else if (MP_HAS(S_MP_MUL_HIGH_COMBA)) {
+      if ((err = s_mp_mul_high_comba(&q, mu, &q, um)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   } else {
+      err = MP_VAL;
+      goto LBL_ERR;
+   }
+
+   /* q3 = q2 / b**(k+1) */
+   mp_rshd(&q, um + 1);
+
+   /* x = x mod b**(k+1), quick (no division) */
+   if ((err = mp_mod_2d(x, MP_DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   /* q = q * m mod b**(k+1), quick (no division) */
+   if ((err = s_mp_mul(&q, m, &q, um + 1)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   /* x = x - q */
+   if ((err = mp_sub(x, &q, x)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   /* If x < 0, add b**(k+1) to it */
+   if (mp_cmp_d(x, 0uL) == MP_LT) {
+      mp_set(&q, 1uL);
+      if ((err = mp_lshd(&q, um + 1)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_add(x, &q, x)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+
+   /* Back off if it's too big */
+   while (mp_cmp(x, m) != MP_LT) {
+      if ((err = s_mp_sub(x, m, x)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+
+LBL_ERR:
+   mp_clear(&q);
+
+   return err;
+}
+#endif

src/libtommath/mp_reduce_2k.c

@@ -0,0 +1,49 @@
+#include "tommath_private.h"
+#ifdef MP_REDUCE_2K_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* reduces a modulo n where n is of the form 2**p - d */
+mp_err mp_reduce_2k(mp_int *a, const mp_int *n, mp_digit d)
+{
+   mp_int q;
+   mp_err err;
+   int p;
+
+   if ((err = mp_init(&q)) != MP_OKAY) {
+      return err;
+   }
+
+   p = mp_count_bits(n);
+   for (;;) {
+      /* q = a/2**p, a = a mod 2**p */
+      if ((err = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      if (d != 1u) {
+         /* q = q * d */
+         if ((err = mp_mul_d(&q, d, &q)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+      }
+
+      /* a = a + q */
+      if ((err = s_mp_add(a, &q, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      if (mp_cmp_mag(a, n) == MP_LT) {
+         break;
+      }
+      if ((err = s_mp_sub(a, n, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+
+LBL_ERR:
+   mp_clear(&q);
+   return err;
+}
+
+#endif

src/libtommath/mp_reduce_2k_l.c

@@ -0,0 +1,52 @@
+#include "tommath_private.h"
+#ifdef MP_REDUCE_2K_L_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* reduces a modulo n where n is of the form 2**p - d
+   This differs from reduce_2k since "d" can be larger
+   than a single digit.
+*/
+mp_err mp_reduce_2k_l(mp_int *a, const mp_int *n, const mp_int *d)
+{
+   mp_int q;
+   mp_err err;
+   int    p;
+
+   if ((err = mp_init(&q)) != MP_OKAY) {
+      return err;
+   }
+
+   p = mp_count_bits(n);
+
+   for (;;) {
+      /* q = a/2**p, a = a mod 2**p */
+      if ((err = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* q = q * d */
+      if ((err = mp_mul(&q, d, &q)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* a = a + q */
+      if ((err = s_mp_add(a, &q, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      if (mp_cmp_mag(a, n) == MP_LT) {
+         break;
+      }
+      if ((err = s_mp_sub(a, n, a)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+   }
+
+LBL_ERR:
+   mp_clear(&q);
+   return err;
+}
+
+#endif

src/libtommath/mp_reduce_2k_setup.c

@@ -0,0 +1,30 @@
+#include "tommath_private.h"
+#ifdef MP_REDUCE_2K_SETUP_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* determines the setup value */
+mp_err mp_reduce_2k_setup(const mp_int *a, mp_digit *d)
+{
+   mp_err err;
+   mp_int tmp;
+
+   if ((err = mp_init(&tmp)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   if ((err = s_mp_sub(&tmp, a, &tmp)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   *d = tmp.dp[0];
+
+LBL_ERR:
+   mp_clear(&tmp);
+   return err;
+}
+#endif

src/libtommath/mp_reduce_2k_setup_l.c

@@ -0,0 +1,28 @@
+#include "tommath_private.h"
+#ifdef MP_REDUCE_2K_SETUP_L_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* determines the setup value */
+mp_err mp_reduce_2k_setup_l(const mp_int *a, mp_int *d)
+{
+   mp_err err;
+   mp_int tmp;
+
+   if ((err = mp_init(&tmp)) != MP_OKAY) {
+      return err;
+   }
+
+   if ((err = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   if ((err = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+LBL_ERR:
+   mp_clear(&tmp);
+   return err;
+}
+#endif