Remove all libtommath source-files, which are not used in Tcl. Don't bother about them any more

91 files changed, 0 insertions, 5678 deletions

diff --git a/libtommath/bn_cutoffs.c b/libtommath/bn_cutoffs.c
deleted file mode 100644
index b02ab71..0000000
--- a/libtommath/bn_cutoffs.c
+++ /dev/null
@@ -1,14 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_CUTOFFS_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-#ifndef MP_FIXED_CUTOFFS
-#include "tommath_cutoffs.h"
-int KARATSUBA_MUL_CUTOFF = MP_DEFAULT_KARATSUBA_MUL_CUTOFF,
-    KARATSUBA_SQR_CUTOFF = MP_DEFAULT_KARATSUBA_SQR_CUTOFF,
-    TOOM_MUL_CUTOFF = MP_DEFAULT_TOOM_MUL_CUTOFF,
-    TOOM_SQR_CUTOFF = MP_DEFAULT_TOOM_SQR_CUTOFF;
-#endif
-
-#endif
diff --git a/libtommath/bn_deprecated.c b/libtommath/bn_deprecated.c
deleted file mode 100644
index a4004f6..0000000
--- a/libtommath/bn_deprecated.c
+++ /dev/null
@@ -1,321 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_DEPRECATED_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-#ifdef BN_MP_GET_BIT_C
-int mp_get_bit(const mp_int *a, int b)
-{
-   if (b < 0) {
-      return MP_VAL;
-   }
-   return (s_mp_get_bit(a, (unsigned int)b) == MP_YES) ? MP_YES : MP_NO;
-}
-#endif
-#ifdef BN_MP_JACOBI_C
-mp_err mp_jacobi(const mp_int *a, const mp_int *n, int *c)
-{
-   if (a->sign == MP_NEG) {
-      return MP_VAL;
-   }
-   if (mp_cmp_d(n, 0uL) != MP_GT) {
-      return MP_VAL;
-   }
-   return mp_kronecker(a, n, c);
-}
-#endif
-#ifdef BN_MP_PRIME_RANDOM_EX_C
-mp_err mp_prime_random_ex(mp_int *a, int t, int size, int flags, private_mp_prime_callback cb, void *dat)
-{
-   return s_mp_prime_random_ex(a, t, size, flags, cb, dat);
-}
-#endif
-#ifdef BN_MP_RAND_DIGIT_C
-mp_err mp_rand_digit(mp_digit *r)
-{
-   mp_err err = s_mp_rand_source(r, sizeof(mp_digit));
-   *r &= MP_MASK;
-   return err;
-}
-#endif
-#ifdef BN_FAST_MP_INVMOD_C
-mp_err fast_mp_invmod(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return s_mp_invmod_fast(a, b, c);
-}
-#endif
-#ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
-mp_err fast_mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho)
-{
-   return s_mp_montgomery_reduce_fast(x, n, rho);
-}
-#endif
-#ifdef BN_FAST_S_MP_MUL_DIGS_C
-mp_err fast_s_mp_mul_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
-{
-   return s_mp_mul_digs_fast(a, b, c, digs);
-}
-#endif
-#ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
-mp_err fast_s_mp_mul_high_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
-{
-   return s_mp_mul_high_digs_fast(a, b, c, digs);
-}
-#endif
-#ifdef BN_FAST_S_MP_SQR_C
-mp_err fast_s_mp_sqr(const mp_int *a, mp_int *b)
-{
-   return s_mp_sqr_fast(a, b);
-}
-#endif
-#ifdef BN_MP_BALANCE_MUL_C
-mp_err mp_balance_mul(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return s_mp_balance_mul(a, b, c);
-}
-#endif
-#ifdef BN_MP_EXPTMOD_FAST_C
-mp_err mp_exptmod_fast(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
-{
-   return s_mp_exptmod_fast(G, X, P, Y, redmode);
-}
-#endif
-#ifdef BN_MP_INVMOD_SLOW_C
-mp_err mp_invmod_slow(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return s_mp_invmod_slow(a, b, c);
-}
-#endif
-#ifdef BN_MP_KARATSUBA_MUL_C
-mp_err mp_karatsuba_mul(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return s_mp_karatsuba_mul(a, b, c);
-}
-#endif
-#ifdef BN_MP_KARATSUBA_SQR_C
-mp_err mp_karatsuba_sqr(const mp_int *a, mp_int *b)
-{
-   return s_mp_karatsuba_sqr(a, b);
-}
-#endif
-#ifdef BN_MP_TOOM_MUL_C
-mp_err mp_toom_mul(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return s_mp_toom_mul(a, b, c);
-}
-#endif
-#ifdef BN_MP_TOOM_SQR_C
-mp_err mp_toom_sqr(const mp_int *a, mp_int *b)
-{
-   return s_mp_toom_sqr(a, b);
-}
-#endif
-#ifdef S_MP_REVERSE_C
-void bn_reverse(unsigned char *s, int len)
-{
-   if (len > 0) {
-      s_mp_reverse(s, (size_t)len);
-   }
-}
-#endif
-#ifdef BN_MP_TC_AND_C
-mp_err mp_tc_and(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return mp_and(a, b, c);
-}
-#endif
-#ifdef BN_MP_TC_OR_C
-mp_err mp_tc_or(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return mp_or(a, b, c);
-}
-#endif
-#ifdef BN_MP_TC_XOR_C
-mp_err mp_tc_xor(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   return mp_xor(a, b, c);
-}
-#endif
-#ifdef BN_MP_TC_DIV_2D_C
-mp_err mp_tc_div_2d(const mp_int *a, int b, mp_int *c)
-{
-   return mp_signed_rsh(a, b, c);
-}
-#endif
-#ifdef BN_MP_INIT_SET_INT_C
-mp_err mp_init_set_int(mp_int *a, unsigned long b)
-{
-   return mp_init_u32(a, (uint32_t)b);
-}
-#endif
-#ifdef BN_MP_SET_INT_C
-mp_err mp_set_int(mp_int *a, unsigned long b)
-{
-   mp_set_u32(a, (uint32_t)b);
-   return MP_OKAY;
-}
-#endif
-#ifdef BN_MP_SET_LONG_C
-mp_err mp_set_long(mp_int *a, unsigned long b)
-{
-   mp_set_u64(a, b);
-   return MP_OKAY;
-}
-#endif
-#ifdef BN_MP_SET_LONG_LONG_C
-mp_err mp_set_long_long(mp_int *a, unsigned long long b)
-{
-   mp_set_u64(a, b);
-   return MP_OKAY;
-}
-#endif
-#ifdef BN_MP_GET_INT_C
-unsigned long mp_get_int(const mp_int *a)
-{
-   return (unsigned long)mp_get_mag_u32(a);
-}
-#endif
-#ifdef BN_MP_GET_LONG_C
-unsigned long mp_get_long(const mp_int *a)
-{
-   return (unsigned long)mp_get_mag_ul(a);
-}
-#endif
-#ifdef BN_MP_GET_LONG_LONG_C
-unsigned long long mp_get_long_long(const mp_int *a)
-{
-   return mp_get_mag_ull(a);
-}
-#endif
-#ifdef BN_MP_PRIME_IS_DIVISIBLE_C
-mp_err mp_prime_is_divisible(const mp_int *a, mp_bool *result)
-{
-   return s_mp_prime_is_divisible(a, result);
-}
-#endif
-#ifdef BN_MP_EXPT_D_EX_C
-mp_err mp_expt_d_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
-{
-   (void)fast;
-   if (b > MP_MIN(MP_DIGIT_MAX, UINT32_MAX)) {
-      return MP_VAL;
-   }
-   return mp_expt_u32(a, (uint32_t)b, c);
-}
-#endif
-#ifdef BN_MP_EXPT_D_C
-mp_err mp_expt_d(const mp_int *a, mp_digit b, mp_int *c)
-{
-   if (b > MP_MIN(MP_DIGIT_MAX, UINT32_MAX)) {
-      return MP_VAL;
-   }
-   return mp_expt_u32(a, (uint32_t)b, c);
-}
-#endif
-#ifdef BN_MP_N_ROOT_EX_C
-mp_err mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
-{
-   (void)fast;
-   if (b > MP_MIN(MP_DIGIT_MAX, UINT32_MAX)) {
-      return MP_VAL;
-   }
-   return mp_root_u32(a, (unsigned int)b, c);
-}
-#endif
-#ifdef BN_MP_N_ROOT_C
-mp_err mp_n_root(const mp_int *a, mp_digit b, mp_int *c)
-{
-   if (b > MP_MIN(MP_DIGIT_MAX, UINT32_MAX)) {
-      return MP_VAL;
-   }
-   return mp_root_u32(a, (unsigned int)b, c);
-}
-#endif
-#ifdef BN_MP_UNSIGNED_BIN_SIZE_C
-int mp_unsigned_bin_size(const mp_int *a)
-{
-   return (int)mp_ubin_size(a);
-}
-#endif
-#ifdef BN_MP_READ_UNSIGNED_BIN_C
-mp_err mp_read_unsigned_bin(mp_int *a, const unsigned char *b, int c)
-{
-   return mp_from_ubin(a, b, (size_t) c);
-}
-#endif
-#ifdef BN_MP_TO_UNSIGNED_BIN_C
-mp_err mp_to_unsigned_bin(const mp_int *a, unsigned char *b)
-{
-   return mp_to_ubin(a, b, SIZE_MAX, NULL);
-}
-#endif
-#ifdef BN_MP_TO_UNSIGNED_BIN_N_C
-mp_err mp_to_unsigned_bin_n(const mp_int *a, unsigned char *b, unsigned long *outlen)
-{
-   size_t n = mp_ubin_size(a);
-   if (*outlen < (unsigned long)n) {
-      return MP_VAL;
-   }
-   *outlen = (unsigned long)n;
-   return mp_to_ubin(a, b, n, NULL);
-}
-#endif
-#ifdef BN_MP_SIGNED_BIN_SIZE_C
-int mp_signed_bin_size(const mp_int *a)
-{
-   return (int)mp_sbin_size(a);
-}
-#endif
-#ifdef BN_MP_READ_SIGNED_BIN_C
-mp_err mp_read_signed_bin(mp_int *a, const unsigned char *b, int c)
-{
-   return mp_from_sbin(a, b, (size_t) c);
-}
-#endif
-#ifdef BN_MP_TO_SIGNED_BIN_C
-mp_err mp_to_signed_bin(const mp_int *a, unsigned char *b)
-{
-   return mp_to_sbin(a, b, SIZE_MAX, NULL);
-}
-#endif
-#ifdef BN_MP_TO_SIGNED_BIN_N_C
-mp_err mp_to_signed_bin_n(const mp_int *a, unsigned char *b, unsigned long *outlen)
-{
-   size_t n = mp_sbin_size(a);
-   if (*outlen < (unsigned long)n) {
-      return MP_VAL;
-   }
-   *outlen = (unsigned long)n;
-   return mp_to_sbin(a, b, n, NULL);
-}
-#endif
-#ifdef BN_MP_TORADIX_N_C
-mp_err mp_toradix_n(const mp_int *a, char *str, int radix, int maxlen)
-{
-   if (maxlen < 0) {
-      return MP_VAL;
-   }
-   return mp_to_radix(a, str, (size_t)maxlen, NULL, radix);
-}
-#endif
-#ifdef BN_MP_TORADIX_C
-mp_err mp_toradix(const mp_int *a, char *str, int radix)
-{
-   return mp_to_radix(a, str, SIZE_MAX, NULL, radix);
-}
-#endif
-#ifdef BN_MP_IMPORT_C
-mp_err mp_import(mp_int *rop, size_t count, int order, size_t size, int endian, size_t nails,
-                 const void *op)
-{
-   return mp_unpack(rop, count, order, size, endian, nails, op);
-}
-#endif
-#ifdef BN_MP_EXPORT_C
-mp_err mp_export(void *rop, size_t *countp, int order, size_t size,
-                 int endian, size_t nails, const mp_int *op)
-{
-   return mp_pack(rop, SIZE_MAX, countp, order, size, endian, nails, op);
-}
-#endif
-#endif
diff --git a/libtommath/bn_mp_2expt.c b/libtommath/bn_mp_2expt.c
deleted file mode 100644
index 23de0c3..0000000
--- a/libtommath/bn_mp_2expt.c
+++ /dev/null
@@ -1,35 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-
-   if (b < 0) {
-      return MP_VAL;
-   }
-
-   /* 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
diff --git a/libtommath/bn_mp_abs.c b/libtommath/bn_mp_abs.c
deleted file mode 100644
index 00900bb..0000000
--- a/libtommath/bn_mp_abs.c
+++ /dev/null
@@ -1,26 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 (a != 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
diff --git a/libtommath/bn_mp_addmod.c b/libtommath/bn_mp_addmod.c
deleted file mode 100644
index 1dcfb67..0000000
--- a/libtommath/bn_mp_addmod.c
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-   mp_int  t;
-
-   if ((err = mp_init(&t)) != MP_OKAY) {
-      return err;
-   }
-
-   if ((err = mp_add(a, b, &t)) != MP_OKAY) {
-      goto LBL_ERR;
-   }
-   err = mp_mod(&t, c, d);
-
-LBL_ERR:
-   mp_clear(&t);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_complement.c b/libtommath/bn_mp_complement.c
deleted file mode 100644
index fef1423..0000000
--- a/libtommath/bn_mp_complement.c
+++ /dev/null
@@ -1,12 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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_err err = mp_neg(a, b);
-   return (err == MP_OKAY) ? mp_sub_d(b, 1uL, b) : err;
-}
-#endif
diff --git a/libtommath/bn_mp_decr.c b/libtommath/bn_mp_decr.c
deleted file mode 100644
index c6a1572..0000000
--- a/libtommath/bn_mp_decr.c
+++ /dev/null
@@ -1,34 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_DECR_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* Decrement "a" by one like "a--". Changes input! */
-mp_err mp_decr(mp_int *a)
-{
-   if (MP_IS_ZERO(a)) {
-      mp_set(a,1uL);
-      a->sign = MP_NEG;
-      return MP_OKAY;
-   } else if (a->sign == MP_NEG) {
-      mp_err err;
-      a->sign = MP_ZPOS;
-      if ((err = mp_incr(a)) != MP_OKAY) {
-         return err;
-      }
-      /* There is no -0 in LTM */
-      if (!MP_IS_ZERO(a)) {
-         a->sign = MP_NEG;
-      }
-      return MP_OKAY;
-   } else if (a->dp[0] > 1uL) {
-      a->dp[0]--;
-      if (a->dp[0] == 0u) {
-         mp_zero(a);
-      }
-      return MP_OKAY;
-   } else {
-      return mp_sub_d(a, 1uL,a);
-   }
-}
-#endif
diff --git a/libtommath/bn_mp_dr_is_modulus.c b/libtommath/bn_mp_dr_is_modulus.c
deleted file mode 100644
index 83760ea..0000000
--- a/libtommath/bn_mp_dr_is_modulus.c
+++ /dev/null
@@ -1,27 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 */
-mp_bool mp_dr_is_modulus(const mp_int *a)
-{
-   int ix;
-
-   /* must be at least two digits */
-   if (a->used < 2) {
-      return MP_NO;
-   }
-
-   /* 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 MP_NO;
-      }
-   }
-   return MP_YES;
-}
-
-#endif
diff --git a/libtommath/bn_mp_dr_reduce.c b/libtommath/bn_mp_dr_reduce.c
deleted file mode 100644
index ffc33a6..0000000
--- a/libtommath/bn_mp_dr_reduce.c
+++ /dev/null
@@ -1,78 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-   int i, m;
-   mp_word  r;
-   mp_digit mu, *tmpx1, *tmpx2;
-
-   /* m = digits in modulus */
-   m = n->used;
-
-   /* ensure that "x" has at least 2m digits */
-   if (x->alloc < (m + m)) {
-      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.
-    */
-top:
-   /* aliases for digits */
-   /* alias for lower half of x */
-   tmpx1 = x->dp;
-
-   /* alias for upper half of x, or x/B**m */
-   tmpx2 = x->dp + m;
-
-   /* set carry to zero */
-   mu = 0;
-
-   /* compute (x mod B**m) + k * [x/B**m] inline and inplace */
-   for (i = 0; i < m; i++) {
-      r         = ((mp_word)*tmpx2++ * (mp_word)k) + *tmpx1 + mu;
-      *tmpx1++  = (mp_digit)(r & MP_MASK);
-      mu        = (mp_digit)(r >> ((mp_word)MP_DIGIT_BIT));
-   }
-
-   /* set final carry */
-   *tmpx1++ = mu;
-
-   /* zero words above m */
-   MP_ZERO_DIGITS(tmpx1, (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) {
-      if ((err = s_mp_sub(x, n, x)) != MP_OKAY) {
-         return err;
-      }
-      goto top;
-   }
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_mp_dr_setup.c b/libtommath/bn_mp_dr_setup.c
deleted file mode 100644
index 32d5f38..0000000
--- a/libtommath/bn_mp_dr_setup.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_error_to_string.c b/libtommath/bn_mp_error_to_string.c
deleted file mode 100644
index 2e2adb0..0000000
--- a/libtommath/bn_mp_error_to_string.c
+++ /dev/null
@@ -1,27 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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";
-   default:
-      return "Invalid error code";
-   }
-}
-
-#endif
diff --git a/libtommath/bn_mp_exptmod.c b/libtommath/bn_mp_exptmod.c
deleted file mode 100644
index 5f811eb..0000000
--- a/libtommath/bn_mp_exptmod.c
+++ /dev/null
@@ -1,76 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 (P->sign == MP_NEG) {
-      return MP_VAL;
-   }
-
-   /* if exponent X is negative we have to recurse */
-   if (X->sign == MP_NEG) {
-      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) == MP_YES)) {
-      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) == MP_YES)) ? 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) == MP_YES) ? 2 : 0;
-   }
-
-   /* if the modulus is odd or dr != 0 use the montgomery method */
-   if (MP_HAS(S_MP_EXPTMOD_FAST) && (MP_IS_ODD(P) || (dr != 0))) {
-      return s_mp_exptmod_fast(G, X, P, Y, dr);
-   } else if (MP_HAS(S_MP_EXPTMOD)) {
-      /* otherwise use the generic Barrett reduction technique */
-      return s_mp_exptmod(G, X, P, Y, 0);
-   } else {
-      /* no exptmod for evens */
-      return MP_VAL;
-   }
-}
-
-#endif
diff --git a/libtommath/bn_mp_exteuclid.c b/libtommath/bn_mp_exteuclid.c
deleted file mode 100644
index faf47ba..0000000
--- a/libtommath/bn_mp_exteuclid.c
+++ /dev/null
@@ -1,73 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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_IS_ZERO(&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 (u3.sign == MP_NEG) {
-      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);
-   }
-
-   err = MP_OKAY;
-LBL_ERR:
-   mp_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_fread.c b/libtommath/bn_mp_fread.c
deleted file mode 100644
index 52ea773..0000000
--- a/libtommath/bn_mp_fread.c
+++ /dev/null
@@ -1,60 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 neg;
-
-   /* if first digit is - then set negative */
-   int ch = fgetc(stream);
-   if (ch == (int)'-') {
-      neg = MP_NEG;
-      ch = fgetc(stream);
-   } else {
-      neg = MP_ZPOS;
-   }
-
-   /* no digits, return error */
-   if (ch == EOF) {
-      return MP_ERR;
-   }
-
-   /* clear a */
-   mp_zero(a);
-
-   do {
-      int y;
-      unsigned pos = (unsigned)(ch - (int)'(');
-      if (mp_s_rmap_reverse_sz < pos) {
-         break;
-      }
-
-      y = (int)mp_s_rmap_reverse[pos];
-
-      if ((y == 0xff) || (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, (mp_digit)y, a)) != MP_OKAY) {
-         return err;
-      }
-   } while ((ch = fgetc(stream)) != EOF);
-
-   if (a->used != 0) {
-      a->sign = neg;
-   }
-
-   return MP_OKAY;
-}
-#endif
-
-#endif
diff --git a/libtommath/bn_mp_from_sbin.c b/libtommath/bn_mp_from_sbin.c
deleted file mode 100644
index 20e4597..0000000
--- a/libtommath/bn_mp_from_sbin.c
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 unsigned char *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 */
-   if (buf[0] == (unsigned char)0) {
-      a->sign = MP_ZPOS;
-   } else {
-      a->sign = MP_NEG;
-   }
-
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_mp_from_ubin.c b/libtommath/bn_mp_from_ubin.c
deleted file mode 100644
index 7f73cbc..0000000
--- a/libtommath/bn_mp_from_ubin.c
+++ /dev/null
@@ -1,39 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_FROM_UBIN_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* reads a unsigned char array, assumes the msb is stored first [big endian] */
-mp_err mp_from_ubin(mp_int *a, const unsigned char *buf, size_t size)
-{
-   mp_err err;
-
-   /* make sure there are at least two digits */
-   if (a->alloc < 2) {
-      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;
-      }
-
-#ifndef MP_8BIT
-      a->dp[0] |= *buf++;
-      a->used += 1;
-#else
-      a->dp[0] = (*buf & MP_MASK);
-      a->dp[1] |= ((*buf++ >> 7) & 1u);
-      a->used += 2;
-#endif
-   }
-   mp_clamp(a);
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_mp_fwrite.c b/libtommath/bn_mp_fwrite.c
deleted file mode 100644
index abe2e67..0000000
--- a/libtommath/bn_mp_fwrite.c
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-   int len;
-   size_t written;
-
-   /* TODO: this function is not in this PR */
-   if (MP_HAS(MP_RADIX_SIZE_OVERESTIMATE)) {
-      /* if ((err = mp_radix_size_overestimate(&t, base, &len)) != MP_OKAY)      goto LBL_ERR; */
-   } else {
-      if ((err = mp_radix_size(a, radix, &len)) != MP_OKAY) {
-         return err;
-      }
-   }
-
-   buf = (char *) MP_MALLOC((size_t)len);
-   if (buf == NULL) {
-      return MP_MEM;
-   }
-
-   if ((err = mp_to_radix(a, buf, (size_t)len, &written, radix)) != MP_OKAY) {
-      goto LBL_ERR;
-   }
-
-   if (fwrite(buf, written, 1uL, stream) != 1uL) {
-      err = MP_ERR;
-      goto LBL_ERR;
-   }
-   err = MP_OKAY;
-
-
-LBL_ERR:
-   MP_FREE_BUFFER(buf, (size_t)len);
-   return err;
-}
-#endif
-
-#endif
diff --git a/libtommath/bn_mp_gcd.c b/libtommath/bn_mp_gcd.c
deleted file mode 100644
index 53029ba..0000000
--- a/libtommath/bn_mp_gcd.c
+++ /dev/null
@@ -1,92 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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_IS_ZERO(a)) {
-      return mp_abs(b, c);
-   }
-   if (MP_IS_ZERO(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_IS_ZERO(&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
diff --git a/libtommath/bn_mp_get_double.c b/libtommath/bn_mp_get_double.c
deleted file mode 100644
index c9b1b19..0000000
--- a/libtommath/bn_mp_get_double.c
+++ /dev/null
@@ -1,18 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 (a->sign == MP_NEG) ? -d : d;
-}
-#endif
diff --git a/libtommath/bn_mp_get_i32.c b/libtommath/bn_mp_get_i32.c
deleted file mode 100644
index 030b657..0000000
--- a/libtommath/bn_mp_get_i32.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_get_i64.c b/libtommath/bn_mp_get_i64.c
deleted file mode 100644
index 969c8d2..0000000
--- a/libtommath/bn_mp_get_i64.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_get_l.c b/libtommath/bn_mp_get_l.c
deleted file mode 100644
index 55d78ec..0000000
--- a/libtommath/bn_mp_get_l.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_get_mag_u32.c b/libtommath/bn_mp_get_mag_u32.c
deleted file mode 100644
index d77189b..0000000
--- a/libtommath/bn_mp_get_mag_u32.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_get_mag_u64.c b/libtommath/bn_mp_get_mag_u64.c
deleted file mode 100644
index 36dd73f..0000000
--- a/libtommath/bn_mp_get_mag_u64.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_get_mag_ul.c b/libtommath/bn_mp_get_mag_ul.c
deleted file mode 100644
index e8819ae..0000000
--- a/libtommath/bn_mp_get_mag_ul.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_incr.c b/libtommath/bn_mp_incr.c
deleted file mode 100644
index 7695ac7..0000000
--- a/libtommath/bn_mp_incr.c
+++ /dev/null
@@ -1,30 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_INCR_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* Increment "a" by one like "a++". Changes input! */
-mp_err mp_incr(mp_int *a)
-{
-   if (MP_IS_ZERO(a)) {
-      mp_set(a,1uL);
-      return MP_OKAY;
-   } else if (a->sign == MP_NEG) {
-      mp_err err;
-      a->sign = MP_ZPOS;
-      if ((err = mp_decr(a)) != MP_OKAY) {
-         return err;
-      }
-      /* There is no -0 in LTM */
-      if (!MP_IS_ZERO(a)) {
-         a->sign = MP_NEG;
-      }
-      return MP_OKAY;
-   } else if (a->dp[0] < MP_DIGIT_MAX) {
-      a->dp[0]++;
-      return MP_OKAY;
-   } else {
-      return mp_add_d(a, 1uL,a);
-   }
-}
-#endif
diff --git a/libtommath/bn_mp_init_i32.c b/libtommath/bn_mp_init_i32.c
deleted file mode 100644
index bc4de8d..0000000
--- a/libtommath/bn_mp_init_i32.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_init_i64.c b/libtommath/bn_mp_init_i64.c
deleted file mode 100644
index 2fa1516..0000000
--- a/libtommath/bn_mp_init_i64.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_init_l.c b/libtommath/bn_mp_init_l.c
deleted file mode 100644
index bc380b5..0000000
--- a/libtommath/bn_mp_init_l.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_init_u32.c b/libtommath/bn_mp_init_u32.c
deleted file mode 100644
index 015d89b..0000000
--- a/libtommath/bn_mp_init_u32.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_init_u64.c b/libtommath/bn_mp_init_u64.c
deleted file mode 100644
index 2b35f7e..0000000
--- a/libtommath/bn_mp_init_u64.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_init_ul.c b/libtommath/bn_mp_init_ul.c
deleted file mode 100644
index 5164f72..0000000
--- a/libtommath/bn_mp_init_ul.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_init_ull.c b/libtommath/bn_mp_init_ull.c
deleted file mode 100644
index 84110c0..0000000
--- a/libtommath/bn_mp_init_ull.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_INIT_ULL_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_INIT_INT(mp_init_ull, mp_set_ull, unsigned long long)
-#endif
diff --git a/libtommath/bn_mp_invmod.c b/libtommath/bn_mp_invmod.c
deleted file mode 100644
index 7b35a24..0000000
--- a/libtommath/bn_mp_invmod.c
+++ /dev/null
@@ -1,23 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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)
-{
-   /* b cannot be negative and has to be >1 */
-   if ((b->sign == MP_NEG) || (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_FAST) && MP_IS_ODD(b)) {
-      return s_mp_invmod_fast(a, b, c);
-   }
-
-   return MP_HAS(S_MP_INVMOD_SLOW)
-          ? s_mp_invmod_slow(a, b, c)
-          : MP_VAL;
-}
-#endif
diff --git a/libtommath/bn_mp_is_square.c b/libtommath/bn_mp_is_square.c
deleted file mode 100644
index 69e77a2..0000000
--- a/libtommath/bn_mp_is_square.c
+++ /dev/null
@@ -1,93 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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, mp_bool *ret)
-{
-   mp_err        err;
-   mp_digit      c;
-   mp_int        t;
-   unsigned long r;
-
-   /* Default to Non-square :) */
-   *ret = MP_NO;
-
-   if (arg->sign == MP_NEG) {
-      return MP_VAL;
-   }
-
-   if (MP_IS_ZERO(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) ? MP_YES : MP_NO;
-LBL_ERR:
-   mp_clear(&t);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_iseven.c b/libtommath/bn_mp_iseven.c
deleted file mode 100644
index 5cb9622..0000000
--- a/libtommath/bn_mp_iseven.c
+++ /dev/null
@@ -1,10 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_ISEVEN_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-mp_bool mp_iseven(const mp_int *a)
-{
-   return MP_IS_EVEN(a) ? MP_YES : MP_NO;
-}
-#endif
diff --git a/libtommath/bn_mp_isodd.c b/libtommath/bn_mp_isodd.c
deleted file mode 100644
index bf17646..0000000
--- a/libtommath/bn_mp_isodd.c
+++ /dev/null
@@ -1,10 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_ISODD_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-mp_bool mp_isodd(const mp_int *a)
-{
-   return MP_IS_ODD(a) ? MP_YES : MP_NO;
-}
-#endif
diff --git a/libtommath/bn_mp_kronecker.c b/libtommath/bn_mp_kronecker.c
deleted file mode 100644
index 525a820..0000000
--- a/libtommath/bn_mp_kronecker.c
+++ /dev/null
@@ -1,129 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 int table[8] = {0, 1, 0, -1, 0, -1, 0, 1};
-
-   if (MP_IS_ZERO(p)) {
-      if ((a->used == 1) && (a->dp[0] == 1u)) {
-         *c = 1;
-      } else {
-         *c = 0;
-      }
-      return MP_OKAY;
-   }
-
-   if (MP_IS_EVEN(a) && MP_IS_EVEN(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 (p1.sign == MP_NEG) {
-      p1.sign = MP_ZPOS;
-      if (a1.sign == MP_NEG) {
-         k = -k;
-      }
-   }
-
-   if ((err = mp_init(&r)) != MP_OKAY) {
-      goto LBL_KRON_1;
-   }
-
-   for (;;) {
-      if (MP_IS_ZERO(&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 (a1.sign == MP_NEG) {
-         /*
-          * 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
diff --git a/libtommath/bn_mp_lcm.c b/libtommath/bn_mp_lcm.c
deleted file mode 100644
index c32b269..0000000
--- a/libtommath/bn_mp_lcm.c
+++ /dev/null
@@ -1,44 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_log_u32.c b/libtommath/bn_mp_log_u32.c
deleted file mode 100644
index 2531cd8..0000000
--- a/libtommath/bn_mp_log_u32.c
+++ /dev/null
@@ -1,180 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_LOG_U32_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* Compute log_{base}(a) */
-static mp_word s_pow(mp_word base, mp_word exponent)
-{
-   mp_word result = 1u;
-   while (exponent != 0u) {
-      if ((exponent & 1u) == 1u) {
-         result *= base;
-      }
-      exponent >>= 1;
-      base *= base;
-   }
-
-   return result;
-}
-
-static mp_digit s_digit_ilogb(mp_digit base, mp_digit n)
-{
-   mp_word bracket_low = 1u, bracket_mid, bracket_high, N;
-   mp_digit ret, high = 1u, low = 0uL, mid;
-
-   if (n < base) {
-      return 0uL;
-   }
-   if (n == base) {
-      return 1uL;
-   }
-
-   bracket_high = (mp_word) base ;
-   N = (mp_word) n;
-
-   while (bracket_high < N) {
-      low = high;
-      bracket_low = bracket_high;
-      high <<= 1;
-      bracket_high *= bracket_high;
-   }
-
-   while (((mp_digit)(high - low)) > 1u) {
-      mid = (low + high) >> 1;
-      bracket_mid = bracket_low * s_pow(base, (mp_word)(mid - low));
-
-      if (N < bracket_mid) {
-         high = mid ;
-         bracket_high = bracket_mid ;
-      }
-      if (N > bracket_mid) {
-         low = mid ;
-         bracket_low = bracket_mid ;
-      }
-      if (N == bracket_mid) {
-         return (mp_digit) mid;
-      }
-   }
-
-   if (bracket_high == N) {
-      ret = high;
-   } else {
-      ret = low;
-   }
-
-   return ret;
-}
-
-/* TODO: output could be "int" because the output of mp_radix_size is int, too,
-         as is the output of mp_bitcount.
-         With the same problem: max size is INT_MAX * MP_DIGIT not INT_MAX only!
-*/
-mp_err mp_log_u32(const mp_int *a, unsigned int base, unsigned int *c)
-{
-   mp_err err;
-   mp_ord cmp;
-   unsigned int high, low, mid;
-   mp_int bracket_low, bracket_high, bracket_mid, t, bi_base;
-
-   err = MP_OKAY;
-
-   if (a->sign == MP_NEG) {
-      return MP_VAL;
-   }
-
-   if (MP_IS_ZERO(a)) {
-      return MP_VAL;
-   }
-
-   if (base < 2u) {
-      return MP_VAL;
-   }
-
-   /* A small shortcut for bases that are powers of two. */
-   if ((base & (base - 1u)) == 0u) {
-      int y, bit_count;
-      for (y=0; (y < 7) && ((base & 1u) == 0u); y++) {
-         base >>= 1;
-      }
-      bit_count = mp_count_bits(a) - 1;
-      *c = (unsigned int)(bit_count/y);
-      return MP_OKAY;
-   }
-
-   if (a->used == 1) {
-      *c = (unsigned int)s_digit_ilogb(base, a->dp[0]);
-      return err;
-   }
-
-   cmp = mp_cmp_d(a, base);
-   if ((cmp == MP_LT) || (cmp == MP_EQ)) {
-      *c = cmp == MP_EQ;
-      return err;
-   }
-
-   if ((err =
-           mp_init_multi(&bracket_low, &bracket_high,
-                         &bracket_mid, &t, &bi_base, NULL)) != MP_OKAY) {
-      return err;
-   }
-
-   low = 0u;
-   mp_set(&bracket_low, 1uL);
-   high = 1u;
-
-   mp_set(&bracket_high, base);
-
-   /*
-       A kind of Giant-step/baby-step algorithm.
-       Idea shamelessly stolen from https://programmingpraxis.com/2010/05/07/integer-logarithms/2/
-       The effect is asymptotic, hence needs benchmarks to test if the Giant-step should be skipped
-       for small n.
-    */
-   while (mp_cmp(&bracket_high, a) == MP_LT) {
-      low = high;
-      if ((err = mp_copy(&bracket_high, &bracket_low)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-      high <<= 1;
-      if ((err = mp_sqr(&bracket_high, &bracket_high)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-   }
-   mp_set(&bi_base, base);
-
-   while ((high - low) > 1u) {
-      mid = (high + low) >> 1;
-
-      if ((err = mp_expt_u32(&bi_base, mid - low, &t)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-      if ((err = mp_mul(&bracket_low, &t, &bracket_mid)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-      cmp = mp_cmp(a, &bracket_mid);
-      if (cmp == MP_LT) {
-         high = mid;
-         mp_exch(&bracket_mid, &bracket_high);
-      }
-      if (cmp == MP_GT) {
-         low = mid;
-         mp_exch(&bracket_mid, &bracket_low);
-      }
-      if (cmp == MP_EQ) {
-         *c = mid;
-         goto LBL_END;
-      }
-   }
-
-   *c = (mp_cmp(&bracket_high, a) == MP_EQ) ? high : low;
-
-LBL_END:
-LBL_ERR:
-   mp_clear_multi(&bracket_low, &bracket_high, &bracket_mid,
-                  &t, &bi_base, NULL);
-   return err;
-}
-
-
-#endif
diff --git a/libtommath/bn_mp_mod_d.c b/libtommath/bn_mp_mod_d.c
deleted file mode 100644
index 0b6c12a..0000000
--- a/libtommath/bn_mp_mod_d.c
+++ /dev/null
@@ -1,10 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_MOD_D_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-mp_err mp_mod_d(const mp_int *a, mp_digit b, mp_digit *c)
-{
-   return mp_div_d(a, b, NULL, c);
-}
-#endif
diff --git a/libtommath/bn_mp_montgomery_calc_normalization.c b/libtommath/bn_mp_montgomery_calc_normalization.c
deleted file mode 100644
index 8379789..0000000
--- a/libtommath/bn_mp_montgomery_calc_normalization.c
+++ /dev/null
@@ -1,44 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_montgomery_reduce.c b/libtommath/bn_mp_montgomery_reduce.c
deleted file mode 100644
index ffe8341..0000000
--- a/libtommath/bn_mp_montgomery_reduce.c
+++ /dev/null
@@ -1,102 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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)
-{
-   int      ix, digs;
-   mp_err   err;
-   mp_digit mu;
-
-   /* 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_MAXFAST)) {
-      return s_mp_montgomery_reduce_fast(x, n, rho);
-   }
-
-   /* grow the input as required */
-   if (x->alloc < digs) {
-      if ((err = mp_grow(x, digs)) != MP_OKAY) {
-         return err;
-      }
-   }
-   x->used = digs;
-
-   for (ix = 0; ix < n->used; ix++) {
-      /* 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 */
-      {
-         int iy;
-         mp_digit *tmpn, *tmpx, u;
-         mp_word r;
-
-         /* alias for digits of the modulus */
-         tmpn = n->dp;
-
-         /* alias for the digits of x [the input] */
-         tmpx = x->dp + ix;
-
-         /* set the carry to zero */
-         u = 0;
-
-         /* Multiply and add in place */
-         for (iy = 0; iy < n->used; iy++) {
-            /* compute product and sum */
-            r       = ((mp_word)mu * (mp_word)*tmpn++) +
-                      (mp_word)u + (mp_word)*tmpx;
-
-            /* get carry */
-            u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
-
-            /* fix digit */
-            *tmpx++ = (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) {
-            *tmpx   += u;
-            u        = *tmpx >> MP_DIGIT_BIT;
-            *tmpx++ &= MP_MASK;
-         }
-      }
-   }
-
-   /* 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
diff --git a/libtommath/bn_mp_montgomery_setup.c b/libtommath/bn_mp_montgomery_setup.c
deleted file mode 100644
index 39f6e9d..0000000
--- a/libtommath/bn_mp_montgomery_setup.c
+++ /dev/null
@@ -1,42 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 */
-#if !defined(MP_8BIT)
-   x *= 2u - (b * x);              /* here x*a==1 mod 2**16 */
-#endif
-#if defined(MP_64BIT) || !(defined(MP_8BIT) || 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
diff --git a/libtommath/bn_mp_mulmod.c b/libtommath/bn_mp_mulmod.c
deleted file mode 100644
index 160d162..0000000
--- a/libtommath/bn_mp_mulmod.c
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-   mp_int t;
-
-   if ((err = mp_init_size(&t, c->used)) != MP_OKAY) {
-      return err;
-   }
-
-   if ((err = mp_mul(a, b, &t)) != MP_OKAY) {
-      goto LBL_ERR;
-   }
-   err = mp_mod(&t, c, d);
-
-LBL_ERR:
-   mp_clear(&t);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_prime_fermat.c b/libtommath/bn_mp_prime_fermat.c
deleted file mode 100644
index af3e884..0000000
--- a/libtommath/bn_mp_prime_fermat.c
+++ /dev/null
@@ -1,47 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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, mp_bool *result)
-{
-   mp_int  t;
-   mp_err  err;
-
-   /* default to composite  */
-   *result = MP_NO;
-
-   /* 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_T;
-   }
-
-   /* is it equal to b? */
-   if (mp_cmp(&t, b) == MP_EQ) {
-      *result = MP_YES;
-   }
-
-   err = MP_OKAY;
-LBL_T:
-   mp_clear(&t);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_prime_frobenius_underwood.c b/libtommath/bn_mp_prime_frobenius_underwood.c
deleted file mode 100644
index 253e8d5..0000000
--- a/libtommath/bn_mp_prime_frobenius_underwood.c
+++ /dev/null
@@ -1,132 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_PRIME_FROBENIUS_UNDERWOOD_C
-
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/*
- *  See file bn_mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
- */
-#ifndef LTM_USE_ONLY_MR
-
-#ifdef MP_8BIT
-/*
- * 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 177
-#else
-#define LTM_FROBENIUS_UNDERWOOD_A 32764
-#endif
-mp_err mp_prime_frobenius_underwood(const mp_int *N, mp_bool *result)
-{
-   mp_int T1z, T2z, Np1z, sz, tz;
-
-   int a, ap2, length, i, j;
-   mp_err err;
-
-   *result = MP_NO;
-
-   if ((err = mp_init_multi(&T1z, &T2z, &Np1z, &sz, &tz, NULL)) != MP_OKAY) {
-      return err;
-   }
-
-   for (a = 0; a < LTM_FROBENIUS_UNDERWOOD_A; a++) {
-      /* 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;
-      }
-      /* (32764^2 - 4) < 2^31, no bigint for >MP_8BIT needed) */
-      mp_set_u32(&T1z, (uint32_t)a);
-
-      if ((err = mp_sqr(&T1z, &T1z)) != MP_OKAY)                  goto LBL_FU_ERR;
-
-      if ((err = mp_sub_d(&T1z, 4uL, &T1z)) != MP_OKAY)           goto LBL_FU_ERR;
-
-      if ((err = mp_kronecker(&T1z, N, &j)) != MP_OKAY)           goto LBL_FU_ERR;
-
-      if (j == -1) {
-         break;
-      }
-
-      if (j == 0) {
-         /* composite */
-         goto LBL_FU_ERR;
-      }
-   }
-   /* Tell it a composite and set return value accordingly */
-   if (a >= LTM_FROBENIUS_UNDERWOOD_A) {
-      err = MP_ITER;
-      goto LBL_FU_ERR;
-   }
-   /* 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_FU_ERR;
-
-   if (!((T1z.used == 1) && (T1z.dp[0] == 1u)))                   goto LBL_FU_ERR;
-
-   ap2 = a + 2;
-   if ((err = mp_add_d(N, 1uL, &Np1z)) != MP_OKAY)                goto LBL_FU_ERR;
-
-   mp_set(&sz, 1uL);
-   mp_set(&tz, 2uL);
-   length = mp_count_bits(&Np1z);
-
-   for (i = length - 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_FU_ERR;
-
-      /* 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_FU_ERR;
-         if ((err = mp_add(&T1z, &T2z, &T2z)) != MP_OKAY)         goto LBL_FU_ERR;
-      }
-
-      if ((err = mp_mul(&T2z, &sz, &T1z)) != MP_OKAY)             goto LBL_FU_ERR;
-      if ((err = mp_sub(&tz, &sz, &T2z)) != MP_OKAY)              goto LBL_FU_ERR;
-      if ((err = mp_add(&sz, &tz, &sz)) != MP_OKAY)               goto LBL_FU_ERR;
-      if ((err = mp_mul(&sz, &T2z, &tz)) != MP_OKAY)              goto LBL_FU_ERR;
-      if ((err = mp_mod(&tz, N, &tz)) != MP_OKAY)                 goto LBL_FU_ERR;
-      if ((err = mp_mod(&T1z, N, &sz)) != MP_OKAY)                goto LBL_FU_ERR;
-      if (s_mp_get_bit(&Np1z, (unsigned int)i) == MP_YES) {
-         /*
-          *  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_FU_ERR;
-         } else {
-            if ((err = mp_mul_d(&sz, (mp_digit)ap2, &T1z)) != MP_OKAY) goto LBL_FU_ERR;
-         }
-         if ((err = mp_add(&T1z, &tz, &T1z)) != MP_OKAY)          goto LBL_FU_ERR;
-         if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY)              goto LBL_FU_ERR;
-         if ((err = mp_sub(&T2z, &sz, &tz)) != MP_OKAY)           goto LBL_FU_ERR;
-         mp_exch(&sz, &T1z);
-      }
-   }
-
-   mp_set_u32(&T1z, (uint32_t)((2 * a) + 5));
-   if ((err = mp_mod(&T1z, N, &T1z)) != MP_OKAY)                  goto LBL_FU_ERR;
-   if (MP_IS_ZERO(&sz) && (mp_cmp(&tz, &T1z) == MP_EQ)) {
-      *result = MP_YES;
-   }
-
-LBL_FU_ERR:
-   mp_clear_multi(&tz, &sz, &Np1z, &T2z, &T1z, NULL);
-   return err;
-}
-
-#endif
-#endif
diff --git a/libtommath/bn_mp_prime_is_prime.c b/libtommath/bn_mp_prime_is_prime.c
deleted file mode 100644
index 7f9fc0b..0000000
--- a/libtommath/bn_mp_prime_is_prime.c
+++ /dev/null
@@ -1,314 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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, mp_bool *result)
-{
-   mp_int  b;
-   int     ix, p_max = 0, size_a, len;
-   mp_bool res;
-   mp_err  err;
-   unsigned int fips_rand, mask;
-
-   /* default to no */
-   *result = MP_NO;
-
-   /* Some shortcuts */
-   /* N > 3 */
-   if (a->used == 1) {
-      if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
-         *result = MP_NO;
-         return MP_OKAY;
-      }
-      if (a->dp[0] == 2u) {
-         *result = MP_YES;
-         return MP_OKAY;
-      }
-   }
-
-   /* N must be odd */
-   if (MP_IS_EVEN(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 != MP_NO) {
-      return MP_OKAY;
-   }
-
-   /* is the input equal to one of the primes in the table? */
-   for (ix = 0; ix < PRIVATE_MP_PRIME_TAB_SIZE; ix++) {
-      if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
-         *result = MP_YES;
-         return MP_OKAY;
-      }
-   }
-#ifdef MP_8BIT
-   /* The search in the loop above was exhaustive in this case */
-   if ((a->used == 1) && (PRIVATE_MP_PRIME_TAB_SIZE >= 31)) {
-      return MP_OKAY;
-   }
-#endif
-
-   /* 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 == MP_YES) {
-      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 == MP_NO) {
-      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 == MP_NO) {
-      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) {
-      /*
-       * Use a Frobenius-Underwood test instead of the Lucas-Selfridge test for
-       * MP_8BIT (It is unknown if the Lucas-Selfridge test works with 16-bit
-       * integers but the necesssary analysis is on the todo-list).
-       */
-#if defined (MP_8BIT) || defined (LTM_USE_FROBENIUS_TEST)
-      err = mp_prime_frobenius_underwood(a, &res);
-      if ((err != MP_OKAY) && (err != MP_ITER)) {
-         goto LBL_B;
-      }
-      if (res == MP_NO) {
-         goto LBL_B;
-      }
-#else
-      if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
-         goto LBL_B;
-      }
-      if (res == MP_NO) {
-         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) {
-      /*
-          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 == MP_NO) {
-            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) {
-      /*
-       * 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++) {
-         /* 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);
-#ifdef MP_8BIT
-         /*
-          * One 8-bit digit is too small, so concatenate two if the size of
-          * unsigned int allows for it.
-          */
-         if ((MP_SIZEOF_BITS(unsigned int)/2) >= MP_SIZEOF_BITS(mp_digit)) {
-            if ((err = mp_rand(&b, 1)) != MP_OKAY) {
-               goto LBL_B;
-            }
-            fips_rand <<= MP_SIZEOF_BITS(mp_digit);
-            fips_rand |= (unsigned int) b.dp[0];
-            fips_rand &= mask;
-         }
-#endif
-         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;
-         }
-         /*
-          * As mentioned above, one 8-bit digit is too small and
-          * although it can only happen in the unlikely case that
-          * an "unsigned int" is smaller than 16 bit a simple test
-          * is cheap and the correction even cheaper.
-          */
-#ifdef MP_8BIT
-         /* All "a" < 2^8 have been caught before */
-         if (len == 1) {
-            len++;
-         }
-#endif
-         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 == MP_NO) {
-            goto LBL_B;
-         }
-      }
-   }
-
-   /* passed the test */
-   *result = MP_YES;
-LBL_B:
-   mp_clear(&b);
-   return err;
-}
-
-#endif
diff --git a/libtommath/bn_mp_prime_miller_rabin.c b/libtommath/bn_mp_prime_miller_rabin.c
deleted file mode 100644
index 96470db..0000000
--- a/libtommath/bn_mp_prime_miller_rabin.c
+++ /dev/null
@@ -1,91 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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, mp_bool *result)
-{
-   mp_int  n1, y, r;
-   mp_err  err;
-   int     s, j;
-
-   /* default */
-   *result = MP_NO;
-
-   /* 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_N1;
-   }
-
-   /* set 2**s * r = n1 */
-   if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
-      goto LBL_N1;
-   }
-
-   /* 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_R;
-   }
-
-   /* compute y = b**r mod a */
-   if ((err = mp_init(&y)) != MP_OKAY) {
-      goto LBL_R;
-   }
-   if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
-      goto LBL_Y;
-   }
-
-   /* 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_Y;
-         }
-
-         /* if y == 1 then composite */
-         if (mp_cmp_d(&y, 1uL) == MP_EQ) {
-            goto LBL_Y;
-         }
-
-         ++j;
-      }
-
-      /* if y != n1 then composite */
-      if (mp_cmp(&y, &n1) != MP_EQ) {
-         goto LBL_Y;
-      }
-   }
-
-   /* probably prime now */
-   *result = MP_YES;
-LBL_Y:
-   mp_clear(&y);
-LBL_R:
-   mp_clear(&r);
-LBL_N1:
-   mp_clear(&n1);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_prime_next_prime.c b/libtommath/bn_mp_prime_next_prime.c
deleted file mode 100644
index d656565..0000000
--- a/libtommath/bn_mp_prime_next_prime.c
+++ /dev/null
@@ -1,132 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 = 1 means the prime must be congruent to 3 mod 4
- */
-mp_err mp_prime_next_prime(mp_int *a, int t, int bbs_style)
-{
-   int      x, y;
-   mp_ord   cmp;
-   mp_err   err;
-   mp_bool  res = MP_NO;
-   mp_digit res_tab[PRIVATE_MP_PRIME_TAB_SIZE], step, 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[PRIVATE_MP_PRIME_TAB_SIZE-1]) == MP_LT) {
-      /* find which prime it is bigger than "a" */
-      for (x = 0; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
-         cmp = mp_cmp_d(a, s_mp_prime_tab[x]);
-         if (cmp == MP_EQ) {
-            continue;
-         }
-         if (cmp != MP_GT) {
-            if ((bbs_style == 1) && ((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? */
-   if (bbs_style == 1) {
-      kstep   = 4;
-   } else {
-      kstep   = 2;
-   }
-
-   /* at this point we will use a combination of a sieve and Miller-Rabin */
-
-   if (bbs_style == 1) {
-      /* 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_IS_EVEN(a)) {
-         /* force odd */
-         if ((err = mp_sub_d(a, 1uL, a)) != MP_OKAY) {
-            return err;
-         }
-      }
-   }
-
-   /* generate the restable */
-   for (x = 1; x < PRIVATE_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 (;;) {
-      /* skip to the next non-trivially divisible candidate */
-      step = 0;
-      do {
-         /* y == 1 if any residue was zero [e.g. cannot be prime] */
-         y     =  0;
-
-         /* increase step to next candidate */
-         step += kstep;
-
-         /* compute the new residue without using division */
-         for (x = 1; x < PRIVATE_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 = 1;
-            }
-         }
-      } while ((y == 1) && (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 == 1) && (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 == MP_YES) {
-         break;
-      }
-   }
-
-   err = MP_OKAY;
-LBL_ERR:
-   mp_clear(&b);
-   return err;
-}
-
-#endif
diff --git a/libtommath/bn_mp_prime_rabin_miller_trials.c b/libtommath/bn_mp_prime_rabin_miller_trials.c
deleted file mode 100644
index 8bbaf6c..0000000
--- a/libtommath/bn_mp_prime_rabin_miller_trials.c
+++ /dev/null
@@ -1,47 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-      } else if (sizes[x].k > size) {
-         return (x == 0) ? sizes[0].t : sizes[x - 1].t;
-      }
-   }
-   return sizes[x-1].t;
-}
-
-
-#endif
diff --git a/libtommath/bn_mp_prime_rand.c b/libtommath/bn_mp_prime_rand.c
deleted file mode 100644
index 4530e9a..0000000
--- a/libtommath/bn_mp_prime_rand.c
+++ /dev/null
@@ -1,141 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 s_mp_prime_random_ex(mp_int *a, int t, int size, int flags, private_mp_prime_callback cb, void *dat)
-{
-   unsigned char *tmp, maskAND, maskOR_msb, maskOR_lsb;
-   int bsize, maskOR_msb_offset;
-   mp_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 = (unsigned char *) MP_MALLOC((size_t)bsize);
-   if (tmp == NULL) {
-      return MP_MEM;
-   }
-
-   /* calc the maskAND value for the MSbyte*/
-   maskAND = ((size&7) == 0) ? 0xFFu : (unsigned char)(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       |= (unsigned char)(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 (cb(tmp, bsize, dat) != bsize) {
-         err = MP_VAL;
-         goto error;
-      }
-
-      /* work over the MSbyte */
-      tmp[0]    &= maskAND;
-      tmp[0]    |= (unsigned char)(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 error;
-      }
-
-      /* is it prime? */
-      if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
-         goto error;
-      }
-      if (res == MP_NO) {
-         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 error;
-         }
-         if ((err = mp_div_2(a, a)) != MP_OKAY) {
-            goto error;
-         }
-
-         /* is it prime? */
-         if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
-            goto error;
-         }
-      }
-   } while (res == MP_NO);
-
-   if ((flags & MP_PRIME_SAFE) != 0) {
-      /* restore a to the original value */
-      if ((err = mp_mul_2(a, a)) != MP_OKAY) {
-         goto error;
-      }
-      if ((err = mp_add_d(a, 1uL, a)) != MP_OKAY) {
-         goto error;
-      }
-   }
-
-   err = MP_OKAY;
-error:
-   MP_FREE_BUFFER(tmp, (size_t)bsize);
-   return err;
-}
-
-static int s_mp_rand_cb(unsigned char *dst, int len, void *dat)
-{
-   (void)dat;
-   if (len <= 0) {
-      return len;
-   }
-   if (s_mp_rand_source(dst, (size_t)len) != MP_OKAY) {
-      return 0;
-   }
-   return len;
-}
-
-mp_err mp_prime_rand(mp_int *a, int t, int size, int flags)
-{
-   return s_mp_prime_random_ex(a, t, size, flags, s_mp_rand_cb, NULL);
-}
-
-#endif
diff --git a/libtommath/bn_mp_prime_strong_lucas_selfridge.c b/libtommath/bn_mp_prime_strong_lucas_selfridge.c
deleted file mode 100644
index b50bbcd..0000000
--- a/libtommath/bn_mp_prime_strong_lucas_selfridge.c
+++ /dev/null
@@ -1,289 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
-
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/*
- *  See file bn_mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
- */
-#ifndef LTM_USE_ONLY_MR
-
-/*
- *  8-bit is just too small. You can try the Frobenius test
- *  but that frobenius test can fail, too, for the same reason.
- */
-#ifndef MP_8BIT
-
-/*
- * 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_mp_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 MP_YES if it is a strong L-S prime, MP_NO 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, mp_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;
-   /* CZ TODO: Some of them need the full 32 bit, hence the (temporary) exclusion of MP_8BIT */
-   int32_t D, Ds, J, sign, P, Q, r, s, u, Nbits;
-   mp_err err;
-   mp_bool oddness;
-
-   *result = MP_NO;
-   /*
-   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, (unsigned int)u) == MP_YES) {
-         /* 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_mp_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_IS_ODD(&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_IS_ODD(&Uz) ? MP_YES : MP_NO;
-         if ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY)               goto LBL_LS_ERR;
-         if ((Uz.sign == MP_NEG) && (oddness != MP_NO)) {
-            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_IS_ODD(&Vz)) {
-            if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY)           goto LBL_LS_ERR;
-         }
-         oddness = MP_IS_ODD(&Vz) ? MP_YES : MP_NO;
-         if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY)               goto LBL_LS_ERR;
-         if ((Vz.sign == MP_NEG) && (oddness != MP_NO)) {
-            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_IS_ZERO(&Uz) || MP_IS_ZERO(&Vz)) {
-      *result = MP_YES;
-      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_IS_ZERO(&Vz)) {
-         *result = MP_YES;
-         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
-#endif
diff --git a/libtommath/bn_mp_rand.c b/libtommath/bn_mp_rand.c
deleted file mode 100644
index 7e9052c..0000000
--- a/libtommath/bn_mp_rand.c
+++ /dev/null
@@ -1,46 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-}
-
-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
diff --git a/libtommath/bn_mp_reduce.c b/libtommath/bn_mp_reduce.c
deleted file mode 100644
index 3c669d4..0000000
--- a/libtommath/bn_mp_reduce.c
+++ /dev/null
@@ -1,83 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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 CLEANUP;
-      }
-   } else if (MP_HAS(S_MP_MUL_HIGH_DIGS)) {
-      if ((err = s_mp_mul_high_digs(&q, mu, &q, um)) != MP_OKAY) {
-         goto CLEANUP;
-      }
-   } else if (MP_HAS(S_MP_MUL_HIGH_DIGS_FAST)) {
-      if ((err = s_mp_mul_high_digs_fast(&q, mu, &q, um)) != MP_OKAY) {
-         goto CLEANUP;
-      }
-   } else {
-      err = MP_VAL;
-      goto CLEANUP;
-   }
-
-   /* 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 CLEANUP;
-   }
-
-   /* q = q * m mod b**(k+1), quick (no division) */
-   if ((err = s_mp_mul_digs(&q, m, &q, um + 1)) != MP_OKAY) {
-      goto CLEANUP;
-   }
-
-   /* x = x - q */
-   if ((err = mp_sub(x, &q, x)) != MP_OKAY) {
-      goto CLEANUP;
-   }
-
-   /* 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 CLEANUP;
-      }
-      if ((err = mp_add(x, &q, x)) != MP_OKAY) {
-         goto CLEANUP;
-      }
-   }
-
-   /* 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 CLEANUP;
-      }
-   }
-
-CLEANUP:
-   mp_clear(&q);
-
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_reduce_2k.c b/libtommath/bn_mp_reduce_2k.c
deleted file mode 100644
index 1cea6cb..0000000
--- a/libtommath/bn_mp_reduce_2k.c
+++ /dev/null
@@ -1,48 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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);
-top:
-   /* 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) {
-      if ((err = s_mp_sub(a, n, a)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-      goto top;
-   }
-
-LBL_ERR:
-   mp_clear(&q);
-   return err;
-}
-
-#endif
diff --git a/libtommath/bn_mp_reduce_2k_l.c b/libtommath/bn_mp_reduce_2k_l.c
deleted file mode 100644
index 6a9f3d3..0000000
--- a/libtommath/bn_mp_reduce_2k_l.c
+++ /dev/null
@@ -1,49 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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);
-top:
-   /* 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) {
-      if ((err = s_mp_sub(a, n, a)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-      goto top;
-   }
-
-LBL_ERR:
-   mp_clear(&q);
-   return err;
-}
-
-#endif
diff --git a/libtommath/bn_mp_reduce_2k_setup.c b/libtommath/bn_mp_reduce_2k_setup.c
deleted file mode 100644
index 2eaf7ad..0000000
--- a/libtommath/bn_mp_reduce_2k_setup.c
+++ /dev/null
@@ -1,32 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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;
-   int    p;
-
-   if ((err = mp_init(&tmp)) != MP_OKAY) {
-      return err;
-   }
-
-   p = mp_count_bits(a);
-   if ((err = mp_2expt(&tmp, p)) != MP_OKAY) {
-      mp_clear(&tmp);
-      return err;
-   }
-
-   if ((err = s_mp_sub(&tmp, a, &tmp)) != MP_OKAY) {
-      mp_clear(&tmp);
-      return err;
-   }
-
-   *d = tmp.dp[0];
-   mp_clear(&tmp);
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_mp_reduce_2k_setup_l.c b/libtommath/bn_mp_reduce_2k_setup_l.c
deleted file mode 100644
index 4f9aa14..0000000
--- a/libtommath/bn_mp_reduce_2k_setup_l.c
+++ /dev/null
@@ -1,28 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_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
diff --git a/libtommath/bn_mp_reduce_is_2k.c b/libtommath/bn_mp_reduce_is_2k.c
deleted file mode 100644
index a9f4f9f..0000000
--- a/libtommath/bn_mp_reduce_is_2k.c
+++ /dev/null
@@ -1,38 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_REDUCE_IS_2K_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* determines if mp_reduce_2k can be used */
-mp_bool mp_reduce_is_2k(const mp_int *a)
-{
-   int ix, iy, iw;
-   mp_digit iz;
-
-   if (a->used == 0) {
-      return MP_NO;
-   } else if (a->used == 1) {
-      return MP_YES;
-   } else if (a->used > 1) {
-      iy = mp_count_bits(a);
-      iz = 1;
-      iw = 1;
-
-      /* Test every bit from the second digit up, must be 1 */
-      for (ix = MP_DIGIT_BIT; ix < iy; ix++) {
-         if ((a->dp[iw] & iz) == 0u) {
-            return MP_NO;
-         }
-         iz <<= 1;
-         if (iz > MP_DIGIT_MAX) {
-            ++iw;
-            iz = 1;
-         }
-      }
-      return MP_YES;
-   } else {
-      return MP_YES;
-   }
-}
-
-#endif
diff --git a/libtommath/bn_mp_reduce_is_2k_l.c b/libtommath/bn_mp_reduce_is_2k_l.c
deleted file mode 100644
index 4bc69be..0000000
--- a/libtommath/bn_mp_reduce_is_2k_l.c
+++ /dev/null
@@ -1,28 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_REDUCE_IS_2K_L_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* determines if reduce_2k_l can be used */
-mp_bool mp_reduce_is_2k_l(const mp_int *a)
-{
-   int ix, iy;
-
-   if (a->used == 0) {
-      return MP_NO;
-   } else if (a->used == 1) {
-      return MP_YES;
-   } else if (a->used > 1) {
-      /* if more than half of the digits are -1 we're sold */
-      for (iy = ix = 0; ix < a->used; ix++) {
-         if (a->dp[ix] == MP_DIGIT_MAX) {
-            ++iy;
-         }
-      }
-      return (iy >= (a->used/2)) ? MP_YES : MP_NO;
-   } else {
-      return MP_NO;
-   }
-}
-
-#endif
diff --git a/libtommath/bn_mp_reduce_setup.c b/libtommath/bn_mp_reduce_setup.c
deleted file mode 100644
index f02160f..0000000
--- a/libtommath/bn_mp_reduce_setup.c
+++ /dev/null
@@ -1,17 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_REDUCE_SETUP_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* pre-calculate the value required for Barrett reduction
- * For a given modulus "b" it calulates the value required in "a"
- */
-mp_err mp_reduce_setup(mp_int *a, const mp_int *b)
-{
-   mp_err err;
-   if ((err = mp_2expt(a, b->used * 2 * MP_DIGIT_BIT)) != MP_OKAY) {
-      return err;
-   }
-   return mp_div(a, b, a, NULL);
-}
-#endif
diff --git a/libtommath/bn_mp_root_u32.c b/libtommath/bn_mp_root_u32.c
deleted file mode 100644
index b60cf26..0000000
--- a/libtommath/bn_mp_root_u32.c
+++ /dev/null
@@ -1,139 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_ROOT_U32_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* find the n'th root of an integer
- *
- * Result found such that (c)**b <= a and (c+1)**b > a
- *
- * This algorithm uses Newton's approximation
- * x[i+1] = x[i] - f(x[i])/f'(x[i])
- * which will find the root in log(N) time where
- * each step involves a fair bit.
- */
-mp_err mp_root_u32(const mp_int *a, unsigned int b, mp_int *c)
-{
-   mp_int t1, t2, t3, a_;
-   mp_ord cmp;
-   int    ilog2;
-   mp_err err;
-
-   /* input must be positive if b is even */
-   if (((b & 1u) == 0u) && (a->sign == MP_NEG)) {
-      return MP_VAL;
-   }
-
-   if ((err = mp_init_multi(&t1, &t2, &t3, NULL)) != MP_OKAY) {
-      return err;
-   }
-
-   /* if a is negative fudge the sign but keep track */
-   a_ = *a;
-   a_.sign = MP_ZPOS;
-
-   /* Compute seed: 2^(log_2(n)/b + 2)*/
-   ilog2 = mp_count_bits(a);
-
-   /*
-     If "b" is larger than INT_MAX it is also larger than
-     log_2(n) because the bit-length of the "n" is measured
-     with an int and hence the root is always < 2 (two).
-   */
-   if (b > (unsigned int)(INT_MAX/2)) {
-      mp_set(c, 1uL);
-      c->sign = a->sign;
-      err = MP_OKAY;
-      goto LBL_ERR;
-   }
-
-   /* "b" is smaller than INT_MAX, we can cast safely */
-   if (ilog2 < (int)b) {
-      mp_set(c, 1uL);
-      c->sign = a->sign;
-      err = MP_OKAY;
-      goto LBL_ERR;
-   }
-   ilog2 =  ilog2 / ((int)b);
-   if (ilog2 == 0) {
-      mp_set(c, 1uL);
-      c->sign = a->sign;
-      err = MP_OKAY;
-      goto LBL_ERR;
-   }
-   /* Start value must be larger than root */
-   ilog2 += 2;
-   if ((err = mp_2expt(&t2,ilog2)) != MP_OKAY)                    goto LBL_ERR;
-   do {
-      /* t1 = t2 */
-      if ((err = mp_copy(&t2, &t1)) != MP_OKAY)                   goto LBL_ERR;
-
-      /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
-
-      /* t3 = t1**(b-1) */
-      if ((err = mp_expt_u32(&t1, b - 1u, &t3)) != MP_OKAY)       goto LBL_ERR;
-
-      /* numerator */
-      /* t2 = t1**b */
-      if ((err = mp_mul(&t3, &t1, &t2)) != MP_OKAY)               goto LBL_ERR;
-
-      /* t2 = t1**b - a */
-      if ((err = mp_sub(&t2, &a_, &t2)) != MP_OKAY)               goto LBL_ERR;
-
-      /* denominator */
-      /* t3 = t1**(b-1) * b  */
-      if ((err = mp_mul_d(&t3, b, &t3)) != MP_OKAY)               goto LBL_ERR;
-
-      /* t3 = (t1**b - a)/(b * t1**(b-1)) */
-      if ((err = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY)         goto LBL_ERR;
-
-      if ((err = mp_sub(&t1, &t3, &t2)) != MP_OKAY)               goto LBL_ERR;
-
-      /*
-          Number of rounds is at most log_2(root). If it is more it
-          got stuck, so break out of the loop and do the rest manually.
-       */
-      if (ilog2-- == 0) {
-         break;
-      }
-   }  while (mp_cmp(&t1, &t2) != MP_EQ);
-
-   /* result can be off by a few so check */
-   /* Loop beneath can overshoot by one if found root is smaller than actual root */
-   for (;;) {
-      if ((err = mp_expt_u32(&t1, b, &t2)) != MP_OKAY)            goto LBL_ERR;
-      cmp = mp_cmp(&t2, &a_);
-      if (cmp == MP_EQ) {
-         err = MP_OKAY;
-         goto LBL_ERR;
-      }
-      if (cmp == MP_LT) {
-         if ((err = mp_add_d(&t1, 1uL, &t1)) != MP_OKAY)          goto LBL_ERR;
-      } else {
-         break;
-      }
-   }
-   /* correct overshoot from above or from recurrence */
-   for (;;) {
-      if ((err = mp_expt_u32(&t1, b, &t2)) != MP_OKAY)            goto LBL_ERR;
-      if (mp_cmp(&t2, &a_) == MP_GT) {
-         if ((err = mp_sub_d(&t1, 1uL, &t1)) != MP_OKAY)          goto LBL_ERR;
-      } else {
-         break;
-      }
-   }
-
-   /* set the result */
-   mp_exch(&t1, c);
-
-   /* set the sign of the result */
-   c->sign = a->sign;
-
-   err = MP_OKAY;
-
-LBL_ERR:
-   mp_clear_multi(&t1, &t2, &t3, NULL);
-   return err;
-}
-
-#endif
diff --git a/libtommath/bn_mp_sbin_size.c b/libtommath/bn_mp_sbin_size.c
deleted file mode 100644
index e0993d6..0000000
--- a/libtommath/bn_mp_sbin_size.c
+++ /dev/null
@@ -1,11 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SBIN_SIZE_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* get the size for an signed equivalent */
-size_t mp_sbin_size(const mp_int *a)
-{
-   return 1u + mp_ubin_size(a);
-}
-#endif
diff --git a/libtommath/bn_mp_set_double.c b/libtommath/bn_mp_set_double.c
deleted file mode 100644
index 7f1ab75..0000000
--- a/libtommath/bn_mp_set_double.c
+++ /dev/null
@@ -1,47 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_DOUBLE_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-#if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
-mp_err mp_set_double(mp_int *a, double b)
-{
-   uint64_t frac;
-   int exp;
-   mp_err err;
-   union {
-      double   dbl;
-      uint64_t bits;
-   } cast;
-   cast.dbl = b;
-
-   exp = (int)((unsigned)(cast.bits >> 52) & 0x7FFu);
-   frac = (cast.bits & (((uint64_t)1 << 52) - (uint64_t)1)) | ((uint64_t)1 << 52);
-
-   if (exp == 0x7FF) { /* +-inf, NaN */
-      return MP_VAL;
-   }
-   exp -= 1023 + 52;
-
-   mp_set_u64(a, frac);
-
-   err = (exp < 0) ? mp_div_2d(a, -exp, a, NULL) : mp_mul_2d(a, exp, a);
-   if (err != MP_OKAY) {
-      return err;
-   }
-
-   if (((cast.bits >> 63) != 0u) && !MP_IS_ZERO(a)) {
-      a->sign = MP_NEG;
-   }
-
-   return MP_OKAY;
-}
-#else
-/* pragma message() not supported by several compilers (in mostly older but still used versions) */
-#  ifdef _MSC_VER
-#    pragma message("mp_set_double implementation is only available on platforms with IEEE754 floating point format")
-#  else
-#    warning "mp_set_double implementation is only available on platforms with IEEE754 floating point format"
-#  endif
-#endif
-#endif
diff --git a/libtommath/bn_mp_set_i32.c b/libtommath/bn_mp_set_i32.c
deleted file mode 100644
index df4513d..0000000
--- a/libtommath/bn_mp_set_i32.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_I32_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_SIGNED(mp_set_i32, mp_set_u32, int32_t, uint32_t)
-#endif
diff --git a/libtommath/bn_mp_set_i64.c b/libtommath/bn_mp_set_i64.c
deleted file mode 100644
index 395103b..0000000
--- a/libtommath/bn_mp_set_i64.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_I64_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_SIGNED(mp_set_i64, mp_set_u64, int64_t, uint64_t)
-#endif
diff --git a/libtommath/bn_mp_set_l.c b/libtommath/bn_mp_set_l.c
deleted file mode 100644
index 1e445fb..0000000
--- a/libtommath/bn_mp_set_l.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_L_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_SIGNED(mp_set_l, mp_set_ul, long, unsigned long)
-#endif
diff --git a/libtommath/bn_mp_set_ll.c b/libtommath/bn_mp_set_ll.c
deleted file mode 100644
index 3e2324f..0000000
--- a/libtommath/bn_mp_set_ll.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_LL_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_SIGNED(mp_set_ll, mp_set_ull, long long, unsigned long long)
-#endif
diff --git a/libtommath/bn_mp_set_u32.c b/libtommath/bn_mp_set_u32.c
deleted file mode 100644
index 18ba5e1..0000000
--- a/libtommath/bn_mp_set_u32.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_U32_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_UNSIGNED(mp_set_u32, uint32_t)
-#endif
diff --git a/libtommath/bn_mp_set_u64.c b/libtommath/bn_mp_set_u64.c
deleted file mode 100644
index 88fab6c..0000000
--- a/libtommath/bn_mp_set_u64.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_U64_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_UNSIGNED(mp_set_u64, uint64_t)
-#endif
diff --git a/libtommath/bn_mp_set_ul.c b/libtommath/bn_mp_set_ul.c
deleted file mode 100644
index adfd85c..0000000
--- a/libtommath/bn_mp_set_ul.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_UL_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_UNSIGNED(mp_set_ul, unsigned long)
-#endif
diff --git a/libtommath/bn_mp_set_ull.c b/libtommath/bn_mp_set_ull.c
deleted file mode 100644
index 8fbc1bd..0000000
--- a/libtommath/bn_mp_set_ull.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SET_ULL_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-MP_SET_UNSIGNED(mp_set_ull, unsigned long long)
-#endif
diff --git a/libtommath/bn_mp_sqrmod.c b/libtommath/bn_mp_sqrmod.c
deleted file mode 100644
index 626ea2c..0000000
--- a/libtommath/bn_mp_sqrmod.c
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SQRMOD_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* c = a * a (mod b) */
-mp_err mp_sqrmod(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   mp_err  err;
-   mp_int  t;
-
-   if ((err = mp_init(&t)) != MP_OKAY) {
-      return err;
-   }
-
-   if ((err = mp_sqr(a, &t)) != MP_OKAY) {
-      goto LBL_ERR;
-   }
-   err = mp_mod(&t, b, c);
-
-LBL_ERR:
-   mp_clear(&t);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_sqrtmod_prime.c b/libtommath/bn_mp_sqrtmod_prime.c
deleted file mode 100644
index a833ed7..0000000
--- a/libtommath/bn_mp_sqrtmod_prime.c
+++ /dev/null
@@ -1,118 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SQRTMOD_PRIME_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* Tonelli-Shanks algorithm
- * https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
- * https://gmplib.org/list-archives/gmp-discuss/2013-April/005300.html
- *
- */
-
-mp_err mp_sqrtmod_prime(const mp_int *n, const mp_int *prime, mp_int *ret)
-{
-   mp_err err;
-   int legendre;
-   mp_int t1, C, Q, S, Z, M, T, R, two;
-   mp_digit i;
-
-   /* first handle the simple cases */
-   if (mp_cmp_d(n, 0uL) == MP_EQ) {
-      mp_zero(ret);
-      return MP_OKAY;
-   }
-   if (mp_cmp_d(prime, 2uL) == MP_EQ)                            return MP_VAL; /* prime must be odd */
-   if ((err = mp_kronecker(n, prime, &legendre)) != MP_OKAY)        return err;
-   if (legendre == -1)                                           return MP_VAL; /* quadratic non-residue mod prime */
-
-   if ((err = mp_init_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, NULL)) != MP_OKAY) {
-      return err;
-   }
-
-   /* SPECIAL CASE: if prime mod 4 == 3
-    * compute directly: err = n^(prime+1)/4 mod prime
-    * Handbook of Applied Cryptography algorithm 3.36
-    */
-   if ((err = mp_mod_d(prime, 4uL, &i)) != MP_OKAY)               goto cleanup;
-   if (i == 3u) {
-      if ((err = mp_add_d(prime, 1uL, &t1)) != MP_OKAY)           goto cleanup;
-      if ((err = mp_div_2(&t1, &t1)) != MP_OKAY)                  goto cleanup;
-      if ((err = mp_div_2(&t1, &t1)) != MP_OKAY)                  goto cleanup;
-      if ((err = mp_exptmod(n, &t1, prime, ret)) != MP_OKAY)      goto cleanup;
-      err = MP_OKAY;
-      goto cleanup;
-   }
-
-   /* NOW: Tonelli-Shanks algorithm */
-
-   /* factor out powers of 2 from prime-1, defining Q and S as: prime-1 = Q*2^S */
-   if ((err = mp_copy(prime, &Q)) != MP_OKAY)                    goto cleanup;
-   if ((err = mp_sub_d(&Q, 1uL, &Q)) != MP_OKAY)                 goto cleanup;
-   /* Q = prime - 1 */
-   mp_zero(&S);
-   /* S = 0 */
-   while (MP_IS_EVEN(&Q)) {
-      if ((err = mp_div_2(&Q, &Q)) != MP_OKAY)                    goto cleanup;
-      /* Q = Q / 2 */
-      if ((err = mp_add_d(&S, 1uL, &S)) != MP_OKAY)               goto cleanup;
-      /* S = S + 1 */
-   }
-
-   /* find a Z such that the Legendre symbol (Z|prime) == -1 */
-   mp_set_u32(&Z, 2u);
-   /* Z = 2 */
-   for (;;) {
-      if ((err = mp_kronecker(&Z, prime, &legendre)) != MP_OKAY)     goto cleanup;
-      if (legendre == -1) break;
-      if ((err = mp_add_d(&Z, 1uL, &Z)) != MP_OKAY)               goto cleanup;
-      /* Z = Z + 1 */
-   }
-
-   if ((err = mp_exptmod(&Z, &Q, prime, &C)) != MP_OKAY)         goto cleanup;
-   /* C = Z ^ Q mod prime */
-   if ((err = mp_add_d(&Q, 1uL, &t1)) != MP_OKAY)                goto cleanup;
-   if ((err = mp_div_2(&t1, &t1)) != MP_OKAY)                    goto cleanup;
-   /* t1 = (Q + 1) / 2 */
-   if ((err = mp_exptmod(n, &t1, prime, &R)) != MP_OKAY)         goto cleanup;
-   /* R = n ^ ((Q + 1) / 2) mod prime */
-   if ((err = mp_exptmod(n, &Q, prime, &T)) != MP_OKAY)          goto cleanup;
-   /* T = n ^ Q mod prime */
-   if ((err = mp_copy(&S, &M)) != MP_OKAY)                       goto cleanup;
-   /* M = S */
-   mp_set_u32(&two, 2u);
-
-   for (;;) {
-      if ((err = mp_copy(&T, &t1)) != MP_OKAY)                    goto cleanup;
-      i = 0;
-      for (;;) {
-         if (mp_cmp_d(&t1, 1uL) == MP_EQ) break;
-         if ((err = mp_exptmod(&t1, &two, prime, &t1)) != MP_OKAY) goto cleanup;
-         i++;
-      }
-      if (i == 0u) {
-         if ((err = mp_copy(&R, ret)) != MP_OKAY)                  goto cleanup;
-         err = MP_OKAY;
-         goto cleanup;
-      }
-      if ((err = mp_sub_d(&M, i, &t1)) != MP_OKAY)                goto cleanup;
-      if ((err = mp_sub_d(&t1, 1uL, &t1)) != MP_OKAY)             goto cleanup;
-      if ((err = mp_exptmod(&two, &t1, prime, &t1)) != MP_OKAY)   goto cleanup;
-      /* t1 = 2 ^ (M - i - 1) */
-      if ((err = mp_exptmod(&C, &t1, prime, &t1)) != MP_OKAY)     goto cleanup;
-      /* t1 = C ^ (2 ^ (M - i - 1)) mod prime */
-      if ((err = mp_sqrmod(&t1, prime, &C)) != MP_OKAY)           goto cleanup;
-      /* C = (t1 * t1) mod prime */
-      if ((err = mp_mulmod(&R, &t1, prime, &R)) != MP_OKAY)       goto cleanup;
-      /* R = (R * t1) mod prime */
-      if ((err = mp_mulmod(&T, &C, prime, &T)) != MP_OKAY)        goto cleanup;
-      /* T = (T * C) mod prime */
-      mp_set(&M, i);
-      /* M = i */
-   }
-
-cleanup:
-   mp_clear_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, NULL);
-   return err;
-}
-
-#endif
diff --git a/libtommath/bn_mp_submod.c b/libtommath/bn_mp_submod.c
deleted file mode 100644
index 5ebd374..0000000
--- a/libtommath/bn_mp_submod.c
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_SUBMOD_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* d = a - b (mod c) */
-mp_err mp_submod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d)
-{
-   mp_err err;
-   mp_int t;
-
-   if ((err = mp_init(&t)) != MP_OKAY) {
-      return err;
-   }
-
-   if ((err = mp_sub(a, b, &t)) != MP_OKAY) {
-      goto LBL_ERR;
-   }
-   err = mp_mod(&t, c, d);
-
-LBL_ERR:
-   mp_clear(&t);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_mp_to_sbin.c b/libtommath/bn_mp_to_sbin.c
deleted file mode 100644
index dbaf53e..0000000
--- a/libtommath/bn_mp_to_sbin.c
+++ /dev/null
@@ -1,22 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_MP_TO_SBIN_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* store in signed [big endian] format */
-mp_err mp_to_sbin(const mp_int *a, unsigned char *buf, size_t maxlen, size_t *written)
-{
-   mp_err err;
-   if (maxlen == 0u) {
-      return MP_BUF;
-   }
-   if ((err = mp_to_ubin(a, buf + 1, maxlen - 1u, written)) != MP_OKAY) {
-      return err;
-   }
-   if (written != NULL) {
-      (*written)++;
-   }
-   buf[0] = (a->sign == MP_ZPOS) ? (unsigned char)0 : (unsigned char)1;
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_prime_tab.c b/libtommath/bn_prime_tab.c
deleted file mode 100644
index a6c07f8..0000000
--- a/libtommath/bn_prime_tab.c
+++ /dev/null
@@ -1,61 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_PRIME_TAB_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-const mp_digit ltm_prime_tab[] = {
-   0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
-   0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
-   0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
-   0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F,
-#ifndef MP_8BIT
-   0x0083,
-   0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
-   0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
-   0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
-   0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
-
-   0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
-   0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
-   0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
-   0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
-   0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
-   0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
-   0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
-   0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
-
-   0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
-   0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
-   0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
-   0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
-   0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
-   0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
-   0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
-   0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
-
-   0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
-   0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
-   0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
-   0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
-   0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
-   0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
-   0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
-   0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
-#endif
-};
-
-#if defined(__GNUC__) && __GNUC__ >= 4
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-const mp_digit *s_mp_prime_tab = ltm_prime_tab;
-#pragma GCC diagnostic pop
-#elif defined(_MSC_VER) && _MSC_VER >= 1500
-#pragma warning(push)
-#pragma warning(disable: 4996)
-const mp_digit *s_mp_prime_tab = ltm_prime_tab;
-#pragma warning(pop)
-#else
-const mp_digit *s_mp_prime_tab = ltm_prime_tab;
-#endif
-
-#endif
diff --git a/libtommath/bn_s_mp_exptmod.c b/libtommath/bn_s_mp_exptmod.c
deleted file mode 100644
index c3bfa95..0000000
--- a/libtommath/bn_s_mp_exptmod.c
+++ /dev/null
@@ -1,198 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_EXPTMOD_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-#ifdef MP_LOW_MEM
-#   define TAB_SIZE 32
-#   define MAX_WINSIZE 5
-#else
-#   define TAB_SIZE 256
-#   define MAX_WINSIZE 0
-#endif
-
-mp_err s_mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
-{
-   mp_int  M[TAB_SIZE], res, mu;
-   mp_digit buf;
-   mp_err   err;
-   int      bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
-   mp_err(*redux)(mp_int *x, const mp_int *m, const mp_int *mu);
-
-   /* find window size */
-   x = mp_count_bits(X);
-   if (x <= 7) {
-      winsize = 2;
-   } else if (x <= 36) {
-      winsize = 3;
-   } else if (x <= 140) {
-      winsize = 4;
-   } else if (x <= 450) {
-      winsize = 5;
-   } else if (x <= 1303) {
-      winsize = 6;
-   } else if (x <= 3529) {
-      winsize = 7;
-   } else {
-      winsize = 8;
-   }
-
-   winsize = MAX_WINSIZE ? MP_MIN(MAX_WINSIZE, winsize) : winsize;
-
-   /* init M array */
-   /* init first cell */
-   if ((err = mp_init(&M[1])) != MP_OKAY) {
-      return err;
-   }
-
-   /* now init the second half of the array */
-   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
-      if ((err = mp_init(&M[x])) != MP_OKAY) {
-         for (y = 1<<(winsize-1); y < x; y++) {
-            mp_clear(&M[y]);
-         }
-         mp_clear(&M[1]);
-         return err;
-      }
-   }
-
-   /* create mu, used for Barrett reduction */
-   if ((err = mp_init(&mu)) != MP_OKAY)                           goto LBL_M;
-
-   if (redmode == 0) {
-      if ((err = mp_reduce_setup(&mu, P)) != MP_OKAY)             goto LBL_MU;
-      redux = mp_reduce;
-   } else {
-      if ((err = mp_reduce_2k_setup_l(P, &mu)) != MP_OKAY)        goto LBL_MU;
-      redux = mp_reduce_2k_l;
-   }
-
-   /* create M table
-    *
-    * The M table contains powers of the base,
-    * e.g. M[x] = G**x mod P
-    *
-    * The first half of the table is not
-    * computed though accept for M[0] and M[1]
-    */
-   if ((err = mp_mod(G, P, &M[1])) != MP_OKAY)                    goto LBL_MU;
-
-   /* compute the value at M[1<<(winsize-1)] by squaring
-    * M[1] (winsize-1) times
-    */
-   if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_MU;
-
-   for (x = 0; x < (winsize - 1); x++) {
-      /* square it */
-      if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)],
-                        &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_MU;
-
-      /* reduce modulo P */
-      if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, &mu)) != MP_OKAY) goto LBL_MU;
-   }
-
-   /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
-    * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
-    */
-   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
-      if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY)     goto LBL_MU;
-      if ((err = redux(&M[x], P, &mu)) != MP_OKAY)                goto LBL_MU;
-   }
-
-   /* setup result */
-   if ((err = mp_init(&res)) != MP_OKAY)                          goto LBL_MU;
-   mp_set(&res, 1uL);
-
-   /* set initial mode and bit cnt */
-   mode   = 0;
-   bitcnt = 1;
-   buf    = 0;
-   digidx = X->used - 1;
-   bitcpy = 0;
-   bitbuf = 0;
-
-   for (;;) {
-      /* grab next digit as required */
-      if (--bitcnt == 0) {
-         /* if digidx == -1 we are out of digits */
-         if (digidx == -1) {
-            break;
-         }
-         /* read next digit and reset the bitcnt */
-         buf    = X->dp[digidx--];
-         bitcnt = (int)MP_DIGIT_BIT;
-      }
-
-      /* grab the next msb from the exponent */
-      y     = (buf >> (mp_digit)(MP_DIGIT_BIT - 1)) & 1uL;
-      buf <<= (mp_digit)1;
-
-      /* if the bit is zero and mode == 0 then we ignore it
-       * These represent the leading zero bits before the first 1 bit
-       * in the exponent.  Technically this opt is not required but it
-       * does lower the # of trivial squaring/reductions used
-       */
-      if ((mode == 0) && (y == 0)) {
-         continue;
-      }
-
-      /* if the bit is zero and mode == 1 then we square */
-      if ((mode == 1) && (y == 0)) {
-         if ((err = mp_sqr(&res, &res)) != MP_OKAY)               goto LBL_RES;
-         if ((err = redux(&res, P, &mu)) != MP_OKAY)              goto LBL_RES;
-         continue;
-      }
-
-      /* else we add it to the window */
-      bitbuf |= (y << (winsize - ++bitcpy));
-      mode    = 2;
-
-      if (bitcpy == winsize) {
-         /* ok window is filled so square as required and multiply  */
-         /* square first */
-         for (x = 0; x < winsize; x++) {
-            if ((err = mp_sqr(&res, &res)) != MP_OKAY)            goto LBL_RES;
-            if ((err = redux(&res, P, &mu)) != MP_OKAY)           goto LBL_RES;
-         }
-
-         /* then multiply */
-         if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY)  goto LBL_RES;
-         if ((err = redux(&res, P, &mu)) != MP_OKAY)             goto LBL_RES;
-
-         /* empty window and reset */
-         bitcpy = 0;
-         bitbuf = 0;
-         mode   = 1;
-      }
-   }
-
-   /* if bits remain then square/multiply */
-   if ((mode == 2) && (bitcpy > 0)) {
-      /* square then multiply if the bit is set */
-      for (x = 0; x < bitcpy; x++) {
-         if ((err = mp_sqr(&res, &res)) != MP_OKAY)               goto LBL_RES;
-         if ((err = redux(&res, P, &mu)) != MP_OKAY)              goto LBL_RES;
-
-         bitbuf <<= 1;
-         if ((bitbuf & (1 << winsize)) != 0) {
-            /* then multiply */
-            if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY)     goto LBL_RES;
-            if ((err = redux(&res, P, &mu)) != MP_OKAY)           goto LBL_RES;
-         }
-      }
-   }
-
-   mp_exch(&res, Y);
-   err = MP_OKAY;
-LBL_RES:
-   mp_clear(&res);
-LBL_MU:
-   mp_clear(&mu);
-LBL_M:
-   mp_clear(&M[1]);
-   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
-      mp_clear(&M[x]);
-   }
-   return err;
-}
-#endif
diff --git a/libtommath/bn_s_mp_exptmod_fast.c b/libtommath/bn_s_mp_exptmod_fast.c
deleted file mode 100644
index 682ded8..0000000
--- a/libtommath/bn_s_mp_exptmod_fast.c
+++ /dev/null
@@ -1,254 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_EXPTMOD_FAST_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
- *
- * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
- * The value of k changes based on the size of the exponent.
- *
- * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
- */
-
-#ifdef MP_LOW_MEM
-#   define TAB_SIZE 32
-#   define MAX_WINSIZE 5
-#else
-#   define TAB_SIZE 256
-#   define MAX_WINSIZE 0
-#endif
-
-mp_err s_mp_exptmod_fast(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
-{
-   mp_int  M[TAB_SIZE], res;
-   mp_digit buf, mp;
-   int     bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
-   mp_err   err;
-
-   /* use a pointer to the reduction algorithm.  This allows us to use
-    * one of many reduction algorithms without modding the guts of
-    * the code with if statements everywhere.
-    */
-   mp_err(*redux)(mp_int *x, const mp_int *n, mp_digit rho);
-
-   /* find window size */
-   x = mp_count_bits(X);
-   if (x <= 7) {
-      winsize = 2;
-   } else if (x <= 36) {
-      winsize = 3;
-   } else if (x <= 140) {
-      winsize = 4;
-   } else if (x <= 450) {
-      winsize = 5;
-   } else if (x <= 1303) {
-      winsize = 6;
-   } else if (x <= 3529) {
-      winsize = 7;
-   } else {
-      winsize = 8;
-   }
-
-   winsize = MAX_WINSIZE ? MP_MIN(MAX_WINSIZE, winsize) : winsize;
-
-   /* init M array */
-   /* init first cell */
-   if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) {
-      return err;
-   }
-
-   /* now init the second half of the array */
-   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
-      if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
-         for (y = 1<<(winsize-1); y < x; y++) {
-            mp_clear(&M[y]);
-         }
-         mp_clear(&M[1]);
-         return err;
-      }
-   }
-
-   /* determine and setup reduction code */
-   if (redmode == 0) {
-      if (MP_HAS(MP_MONTGOMERY_SETUP)) {
-         /* now setup montgomery  */
-         if ((err = mp_montgomery_setup(P, &mp)) != MP_OKAY)      goto LBL_M;
-      } else {
-         err = MP_VAL;
-         goto LBL_M;
-      }
-
-      /* automatically pick the comba one if available (saves quite a few calls/ifs) */
-      if (MP_HAS(S_MP_MONTGOMERY_REDUCE_FAST) &&
-          (((P->used * 2) + 1) < MP_WARRAY) &&
-          (P->used < MP_MAXFAST)) {
-         redux = s_mp_montgomery_reduce_fast;
-      } else if (MP_HAS(MP_MONTGOMERY_REDUCE)) {
-         /* use slower baseline Montgomery method */
-         redux = mp_montgomery_reduce;
-      } else {
-         err = MP_VAL;
-         goto LBL_M;
-      }
-   } else if (redmode == 1) {
-      if (MP_HAS(MP_DR_SETUP) && MP_HAS(MP_DR_REDUCE)) {
-         /* setup DR reduction for moduli of the form B**k - b */
-         mp_dr_setup(P, &mp);
-         redux = mp_dr_reduce;
-      } else {
-         err = MP_VAL;
-         goto LBL_M;
-      }
-   } else if (MP_HAS(MP_REDUCE_2K_SETUP) && MP_HAS(MP_REDUCE_2K)) {
-      /* setup DR reduction for moduli of the form 2**k - b */
-      if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY)          goto LBL_M;
-      redux = mp_reduce_2k;
-   } else {
-      err = MP_VAL;
-      goto LBL_M;
-   }
-
-   /* setup result */
-   if ((err = mp_init_size(&res, P->alloc)) != MP_OKAY)           goto LBL_M;
-
-   /* create M table
-    *
-
-    *
-    * The first half of the table is not computed though accept for M[0] and M[1]
-    */
-
-   if (redmode == 0) {
-      if (MP_HAS(MP_MONTGOMERY_CALC_NORMALIZATION)) {
-         /* now we need R mod m */
-         if ((err = mp_montgomery_calc_normalization(&res, P)) != MP_OKAY) goto LBL_RES;
-
-         /* now set M[1] to G * R mod m */
-         if ((err = mp_mulmod(G, &res, P, &M[1])) != MP_OKAY)     goto LBL_RES;
-      } else {
-         err = MP_VAL;
-         goto LBL_RES;
-      }
-   } else {
-      mp_set(&res, 1uL);
-      if ((err = mp_mod(G, P, &M[1])) != MP_OKAY)                 goto LBL_RES;
-   }
-
-   /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
-   if ((err = mp_copy(&M[1], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
-
-   for (x = 0; x < (winsize - 1); x++) {
-      if ((err = mp_sqr(&M[(size_t)1 << (winsize - 1)], &M[(size_t)1 << (winsize - 1)])) != MP_OKAY) goto LBL_RES;
-      if ((err = redux(&M[(size_t)1 << (winsize - 1)], P, mp)) != MP_OKAY) goto LBL_RES;
-   }
-
-   /* create upper table */
-   for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
-      if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY)     goto LBL_RES;
-      if ((err = redux(&M[x], P, mp)) != MP_OKAY)                 goto LBL_RES;
-   }
-
-   /* set initial mode and bit cnt */
-   mode   = 0;
-   bitcnt = 1;
-   buf    = 0;
-   digidx = X->used - 1;
-   bitcpy = 0;
-   bitbuf = 0;
-
-   for (;;) {
-      /* grab next digit as required */
-      if (--bitcnt == 0) {
-         /* if digidx == -1 we are out of digits so break */
-         if (digidx == -1) {
-            break;
-         }
-         /* read next digit and reset bitcnt */
-         buf    = X->dp[digidx--];
-         bitcnt = (int)MP_DIGIT_BIT;
-      }
-
-      /* grab the next msb from the exponent */
-      y     = (mp_digit)(buf >> (MP_DIGIT_BIT - 1)) & 1uL;
-      buf <<= (mp_digit)1;
-
-      /* if the bit is zero and mode == 0 then we ignore it
-       * These represent the leading zero bits before the first 1 bit
-       * in the exponent.  Technically this opt is not required but it
-       * does lower the # of trivial squaring/reductions used
-       */
-      if ((mode == 0) && (y == 0)) {
-         continue;
-      }
-
-      /* if the bit is zero and mode == 1 then we square */
-      if ((mode == 1) && (y == 0)) {
-         if ((err = mp_sqr(&res, &res)) != MP_OKAY)               goto LBL_RES;
-         if ((err = redux(&res, P, mp)) != MP_OKAY)               goto LBL_RES;
-         continue;
-      }
-
-      /* else we add it to the window */
-      bitbuf |= (y << (winsize - ++bitcpy));
-      mode    = 2;
-
-      if (bitcpy == winsize) {
-         /* ok window is filled so square as required and multiply  */
-         /* square first */
-         for (x = 0; x < winsize; x++) {
-            if ((err = mp_sqr(&res, &res)) != MP_OKAY)            goto LBL_RES;
-            if ((err = redux(&res, P, mp)) != MP_OKAY)            goto LBL_RES;
-         }
-
-         /* then multiply */
-         if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY)   goto LBL_RES;
-         if ((err = redux(&res, P, mp)) != MP_OKAY)               goto LBL_RES;
-
-         /* empty window and reset */
-         bitcpy = 0;
-         bitbuf = 0;
-         mode   = 1;
-      }
-   }
-
-   /* if bits remain then square/multiply */
-   if ((mode == 2) && (bitcpy > 0)) {
-      /* square then multiply if the bit is set */
-      for (x = 0; x < bitcpy; x++) {
-         if ((err = mp_sqr(&res, &res)) != MP_OKAY)               goto LBL_RES;
-         if ((err = redux(&res, P, mp)) != MP_OKAY)               goto LBL_RES;
-
-         /* get next bit of the window */
-         bitbuf <<= 1;
-         if ((bitbuf & (1 << winsize)) != 0) {
-            /* then multiply */
-            if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY)     goto LBL_RES;
-            if ((err = redux(&res, P, mp)) != MP_OKAY)            goto LBL_RES;
-         }
-      }
-   }
-
-   if (redmode == 0) {
-      /* fixup result if Montgomery reduction is used
-       * recall that any value in a Montgomery system is
-       * actually multiplied by R mod n.  So we have
-       * to reduce one more time to cancel out the factor
-       * of R.
-       */
-      if ((err = redux(&res, P, mp)) != MP_OKAY)                  goto LBL_RES;
-   }
-
-   /* swap res with Y */
-   mp_exch(&res, Y);
-   err = MP_OKAY;
-LBL_RES:
-   mp_clear(&res);
-LBL_M:
-   mp_clear(&M[1]);
-   for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
-      mp_clear(&M[x]);
-   }
-   return err;
-}
-#endif
diff --git a/libtommath/bn_s_mp_get_bit.c b/libtommath/bn_s_mp_get_bit.c
deleted file mode 100644
index 28598df..0000000
--- a/libtommath/bn_s_mp_get_bit.c
+++ /dev/null
@@ -1,21 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_GET_BIT_C
-
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* Get bit at position b and return MP_YES if the bit is 1, MP_NO if it is 0 */
-mp_bool s_mp_get_bit(const mp_int *a, unsigned int b)
-{
-   mp_digit bit;
-   int limb = (int)(b / MP_DIGIT_BIT);
-
-   if (limb >= a->used) {
-      return MP_NO;
-   }
-
-   bit = (mp_digit)1 << (b % MP_DIGIT_BIT);
-   return ((a->dp[limb] & bit) != 0u) ? MP_YES : MP_NO;
-}
-
-#endif
diff --git a/libtommath/bn_s_mp_invmod_fast.c b/libtommath/bn_s_mp_invmod_fast.c
deleted file mode 100644
index 677d7ab..0000000
--- a/libtommath/bn_s_mp_invmod_fast.c
+++ /dev/null
@@ -1,118 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_INVMOD_FAST_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* computes the modular inverse via binary extended euclidean algorithm,
- * that is c = 1/a mod b
- *
- * Based on slow invmod except this is optimized for the case where b is
- * odd as per HAC Note 14.64 on pp. 610
- */
-mp_err s_mp_invmod_fast(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   mp_int  x, y, u, v, B, D;
-   mp_sign neg;
-   mp_err  err;
-
-   /* 2. [modified] b must be odd   */
-   if (MP_IS_EVEN(b)) {
-      return MP_VAL;
-   }
-
-   /* init all our temps */
-   if ((err = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
-      return err;
-   }
-
-   /* x == modulus, y == value to invert */
-   if ((err = mp_copy(b, &x)) != MP_OKAY)                         goto LBL_ERR;
-
-   /* we need y = |a| */
-   if ((err = mp_mod(a, b, &y)) != MP_OKAY)                       goto LBL_ERR;
-
-   /* if one of x,y is zero return an error! */
-   if (MP_IS_ZERO(&x) || MP_IS_ZERO(&y)) {
-      err = MP_VAL;
-      goto LBL_ERR;
-   }
-
-   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
-   if ((err = mp_copy(&x, &u)) != MP_OKAY)                        goto LBL_ERR;
-   if ((err = mp_copy(&y, &v)) != MP_OKAY)                        goto LBL_ERR;
-   mp_set(&D, 1uL);
-
-top:
-   /* 4.  while u is even do */
-   while (MP_IS_EVEN(&u)) {
-      /* 4.1 u = u/2 */
-      if ((err = mp_div_2(&u, &u)) != MP_OKAY)                    goto LBL_ERR;
-
-      /* 4.2 if B is odd then */
-      if (MP_IS_ODD(&B)) {
-         if ((err = mp_sub(&B, &x, &B)) != MP_OKAY)               goto LBL_ERR;
-      }
-      /* B = B/2 */
-      if ((err = mp_div_2(&B, &B)) != MP_OKAY)                    goto LBL_ERR;
-   }
-
-   /* 5.  while v is even do */
-   while (MP_IS_EVEN(&v)) {
-      /* 5.1 v = v/2 */
-      if ((err = mp_div_2(&v, &v)) != MP_OKAY)                    goto LBL_ERR;
-
-      /* 5.2 if D is odd then */
-      if (MP_IS_ODD(&D)) {
-         /* D = (D-x)/2 */
-         if ((err = mp_sub(&D, &x, &D)) != MP_OKAY)               goto LBL_ERR;
-      }
-      /* D = D/2 */
-      if ((err = mp_div_2(&D, &D)) != MP_OKAY)                    goto LBL_ERR;
-   }
-
-   /* 6.  if u >= v then */
-   if (mp_cmp(&u, &v) != MP_LT) {
-      /* u = u - v, B = B - D */
-      if ((err = mp_sub(&u, &v, &u)) != MP_OKAY)                  goto LBL_ERR;
-
-      if ((err = mp_sub(&B, &D, &B)) != MP_OKAY)                  goto LBL_ERR;
-   } else {
-      /* v - v - u, D = D - B */
-      if ((err = mp_sub(&v, &u, &v)) != MP_OKAY)                  goto LBL_ERR;
-
-      if ((err = mp_sub(&D, &B, &D)) != MP_OKAY)                  goto LBL_ERR;
-   }
-
-   /* if not zero goto step 4 */
-   if (!MP_IS_ZERO(&u)) {
-      goto top;
-   }
-
-   /* now a = C, b = D, gcd == g*v */
-
-   /* if v != 1 then there is no inverse */
-   if (mp_cmp_d(&v, 1uL) != MP_EQ) {
-      err = MP_VAL;
-      goto LBL_ERR;
-   }
-
-   /* b is now the inverse */
-   neg = a->sign;
-   while (D.sign == MP_NEG) {
-      if ((err = mp_add(&D, b, &D)) != MP_OKAY)                   goto LBL_ERR;
-   }
-
-   /* too big */
-   while (mp_cmp_mag(&D, b) != MP_LT) {
-      if ((err = mp_sub(&D, b, &D)) != MP_OKAY)                   goto LBL_ERR;
-   }
-
-   mp_exch(&D, c);
-   c->sign = neg;
-   err = MP_OKAY;
-
-LBL_ERR:
-   mp_clear_multi(&x, &y, &u, &v, &B, &D, NULL);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_s_mp_invmod_slow.c b/libtommath/bn_s_mp_invmod_slow.c
deleted file mode 100644
index 4c5db33..0000000
--- a/libtommath/bn_s_mp_invmod_slow.c
+++ /dev/null
@@ -1,119 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_INVMOD_SLOW_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* hac 14.61, pp608 */
-mp_err s_mp_invmod_slow(const mp_int *a, const mp_int *b, mp_int *c)
-{
-   mp_int  x, y, u, v, A, B, C, D;
-   mp_err  err;
-
-   /* b cannot be negative */
-   if ((b->sign == MP_NEG) || MP_IS_ZERO(b)) {
-      return MP_VAL;
-   }
-
-   /* init temps */
-   if ((err = mp_init_multi(&x, &y, &u, &v,
-                            &A, &B, &C, &D, NULL)) != MP_OKAY) {
-      return err;
-   }
-
-   /* x = a, y = b */
-   if ((err = mp_mod(a, b, &x)) != MP_OKAY)                       goto LBL_ERR;
-   if ((err = mp_copy(b, &y)) != MP_OKAY)                         goto LBL_ERR;
-
-   /* 2. [modified] if x,y are both even then return an error! */
-   if (MP_IS_EVEN(&x) && MP_IS_EVEN(&y)) {
-      err = MP_VAL;
-      goto LBL_ERR;
-   }
-
-   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
-   if ((err = mp_copy(&x, &u)) != MP_OKAY)                        goto LBL_ERR;
-   if ((err = mp_copy(&y, &v)) != MP_OKAY)                        goto LBL_ERR;
-   mp_set(&A, 1uL);
-   mp_set(&D, 1uL);
-
-top:
-   /* 4.  while u is even do */
-   while (MP_IS_EVEN(&u)) {
-      /* 4.1 u = u/2 */
-      if ((err = mp_div_2(&u, &u)) != MP_OKAY)                    goto LBL_ERR;
-
-      /* 4.2 if A or B is odd then */
-      if (MP_IS_ODD(&A) || MP_IS_ODD(&B)) {
-         /* A = (A+y)/2, B = (B-x)/2 */
-         if ((err = mp_add(&A, &y, &A)) != MP_OKAY)               goto LBL_ERR;
-         if ((err = mp_sub(&B, &x, &B)) != MP_OKAY)               goto LBL_ERR;
-      }
-      /* A = A/2, B = B/2 */
-      if ((err = mp_div_2(&A, &A)) != MP_OKAY)                    goto LBL_ERR;
-      if ((err = mp_div_2(&B, &B)) != MP_OKAY)                    goto LBL_ERR;
-   }
-
-   /* 5.  while v is even do */
-   while (MP_IS_EVEN(&v)) {
-      /* 5.1 v = v/2 */
-      if ((err = mp_div_2(&v, &v)) != MP_OKAY)                    goto LBL_ERR;
-
-      /* 5.2 if C or D is odd then */
-      if (MP_IS_ODD(&C) || MP_IS_ODD(&D)) {
-         /* C = (C+y)/2, D = (D-x)/2 */
-         if ((err = mp_add(&C, &y, &C)) != MP_OKAY)               goto LBL_ERR;
-         if ((err = mp_sub(&D, &x, &D)) != MP_OKAY)               goto LBL_ERR;
-      }
-      /* C = C/2, D = D/2 */
-      if ((err = mp_div_2(&C, &C)) != MP_OKAY)                    goto LBL_ERR;
-      if ((err = mp_div_2(&D, &D)) != MP_OKAY)                    goto LBL_ERR;
-   }
-
-   /* 6.  if u >= v then */
-   if (mp_cmp(&u, &v) != MP_LT) {
-      /* u = u - v, A = A - C, B = B - D */
-      if ((err = mp_sub(&u, &v, &u)) != MP_OKAY)                  goto LBL_ERR;
-
-      if ((err = mp_sub(&A, &C, &A)) != MP_OKAY)                  goto LBL_ERR;
-
-      if ((err = mp_sub(&B, &D, &B)) != MP_OKAY)                  goto LBL_ERR;
-   } else {
-      /* v - v - u, C = C - A, D = D - B */
-      if ((err = mp_sub(&v, &u, &v)) != MP_OKAY)                  goto LBL_ERR;
-
-      if ((err = mp_sub(&C, &A, &C)) != MP_OKAY)                  goto LBL_ERR;
-
-      if ((err = mp_sub(&D, &B, &D)) != MP_OKAY)                  goto LBL_ERR;
-   }
-
-   /* if not zero goto step 4 */
-   if (!MP_IS_ZERO(&u)) {
-      goto top;
-   }
-
-   /* now a = C, b = D, gcd == g*v */
-
-   /* if v != 1 then there is no inverse */
-   if (mp_cmp_d(&v, 1uL) != MP_EQ) {
-      err = MP_VAL;
-      goto LBL_ERR;
-   }
-
-   /* if its too low */
-   while (mp_cmp_d(&C, 0uL) == MP_LT) {
-      if ((err = mp_add(&C, b, &C)) != MP_OKAY)                   goto LBL_ERR;
-   }
-
-   /* too big */
-   while (mp_cmp_mag(&C, b) != MP_LT) {
-      if ((err = mp_sub(&C, b, &C)) != MP_OKAY)                   goto LBL_ERR;
-   }
-
-   /* C is now the inverse */
-   mp_exch(&C, c);
-   err = MP_OKAY;
-LBL_ERR:
-   mp_clear_multi(&x, &y, &u, &v, &A, &B, &C, &D, NULL);
-   return err;
-}
-#endif
diff --git a/libtommath/bn_s_mp_montgomery_reduce_fast.c b/libtommath/bn_s_mp_montgomery_reduce_fast.c
deleted file mode 100644
index 3f0c672..0000000
--- a/libtommath/bn_s_mp_montgomery_reduce_fast.c
+++ /dev/null
@@ -1,159 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_MONTGOMERY_REDUCE_FAST_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* computes xR**-1 == x (mod N) via Montgomery Reduction
- *
- * This is an optimized implementation of montgomery_reduce
- * which uses the comba method to quickly calculate the columns of the
- * reduction.
- *
- * Based on Algorithm 14.32 on pp.601 of HAC.
-*/
-mp_err s_mp_montgomery_reduce_fast(mp_int *x, const mp_int *n, mp_digit rho)
-{
-   int     ix, olduse;
-   mp_err  err;
-   mp_word W[MP_WARRAY];
-
-   if (x->used > MP_WARRAY) {
-      return MP_VAL;
-   }
-
-   /* get old used count */
-   olduse = x->used;
-
-   /* grow a as required */
-   if (x->alloc < (n->used + 1)) {
-      if ((err = mp_grow(x, n->used + 1)) != MP_OKAY) {
-         return err;
-      }
-   }
-
-   /* first we have to get the digits of the input into
-    * an array of double precision words W[...]
-    */
-   {
-      mp_word *_W;
-      mp_digit *tmpx;
-
-      /* alias for the W[] array */
-      _W   = W;
-
-      /* alias for the digits of  x*/
-      tmpx = x->dp;
-
-      /* copy the digits of a into W[0..a->used-1] */
-      for (ix = 0; ix < x->used; ix++) {
-         *_W++ = *tmpx++;
-      }
-
-      /* zero the high words of W[a->used..m->used*2] */
-      if (ix < ((n->used * 2) + 1)) {
-         MP_ZERO_BUFFER(_W, sizeof(mp_word) * (size_t)(((n->used * 2) + 1) - ix));
-      }
-   }
-
-   /* now we proceed to zero successive digits
-    * from the least significant upwards
-    */
-   for (ix = 0; ix < n->used; ix++) {
-      /* mu = ai * m' mod b
-       *
-       * We avoid a double precision multiplication (which isn't required)
-       * by casting the value down to a mp_digit.  Note this requires
-       * that W[ix-1] have  the carry cleared (see after the inner loop)
-       */
-      mp_digit mu;
-      mu = ((W[ix] & MP_MASK) * rho) & MP_MASK;
-
-      /* a = a + mu * m * b**i
-       *
-       * This is computed in place and on the fly.  The multiplication
-       * by b**i is handled by offseting which columns the results
-       * are added to.
-       *
-       * Note the comba method normally doesn't handle carries in the
-       * inner loop In this case we fix the carry from the previous
-       * column since the Montgomery reduction requires digits of the
-       * result (so far) [see above] to work.  This is
-       * handled by fixing up one carry after the inner loop.  The
-       * carry fixups are done in order so after these loops the
-       * first m->used words of W[] have the carries fixed
-       */
-      {
-         int iy;
-         mp_digit *tmpn;
-         mp_word *_W;
-
-         /* alias for the digits of the modulus */
-         tmpn = n->dp;
-
-         /* Alias for the columns set by an offset of ix */
-         _W = W + ix;
-
-         /* inner loop */
-         for (iy = 0; iy < n->used; iy++) {
-            *_W++ += (mp_word)mu * (mp_word)*tmpn++;
-         }
-      }
-
-      /* now fix carry for next digit, W[ix+1] */
-      W[ix + 1] += W[ix] >> (mp_word)MP_DIGIT_BIT;
-   }
-
-   /* now we have to propagate the carries and
-    * shift the words downward [all those least
-    * significant digits we zeroed].
-    */
-   {
-      mp_digit *tmpx;
-      mp_word *_W, *_W1;
-
-      /* nox fix rest of carries */
-
-      /* alias for current word */
-      _W1 = W + ix;
-
-      /* alias for next word, where the carry goes */
-      _W = W + ++ix;
-
-      for (; ix < ((n->used * 2) + 1); ix++) {
-         *_W++ += *_W1++ >> (mp_word)MP_DIGIT_BIT;
-      }
-
-      /* copy out, A = A/b**n
-       *
-       * The result is A/b**n but instead of converting from an
-       * array of mp_word to mp_digit than calling mp_rshd
-       * we just copy them in the right order
-       */
-
-      /* alias for destination word */
-      tmpx = x->dp;
-
-      /* alias for shifted double precision result */
-      _W = W + n->used;
-
-      for (ix = 0; ix < (n->used + 1); ix++) {
-         *tmpx++ = *_W++ & (mp_word)MP_MASK;
-      }
-
-      /* zero oldused digits, if the input a was larger than
-       * m->used+1 we'll have to clear the digits
-       */
-      MP_ZERO_DIGITS(tmpx, olduse - ix);
-   }
-
-   /* set the max used and clamp */
-   x->used = n->used + 1;
-   mp_clamp(x);
-
-   /* if A >= m then A = A - m */
-   if (mp_cmp_mag(x, n) != MP_LT) {
-      return s_mp_sub(x, n, x);
-   }
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_s_mp_mul_high_digs.c b/libtommath/bn_s_mp_mul_high_digs.c
deleted file mode 100644
index c9dd355..0000000
--- a/libtommath/bn_s_mp_mul_high_digs.c
+++ /dev/null
@@ -1,68 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_MUL_HIGH_DIGS_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* multiplies |a| * |b| and does not compute the lower digs digits
- * [meant to get the higher part of the product]
- */
-mp_err s_mp_mul_high_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
-{
-   mp_int   t;
-   int      pa, pb, ix, iy;
-   mp_err   err;
-   mp_digit u;
-   mp_word  r;
-   mp_digit tmpx, *tmpt, *tmpy;
-
-   if (digs < 0) {
-      return MP_VAL;
-   }
-
-   /* can we use the fast multiplier? */
-   if (MP_HAS(S_MP_MUL_HIGH_DIGS_FAST)
-       && ((a->used + b->used + 1) < MP_WARRAY)
-       && (MP_MIN(a->used, b->used) < MP_MAXFAST)) {
-      return s_mp_mul_high_digs_fast(a, b, c, digs);
-   }
-
-   if ((err = mp_init_size(&t, a->used + b->used + 1)) != MP_OKAY) {
-      return err;
-   }
-   t.used = a->used + b->used + 1;
-
-   pa = a->used;
-   pb = b->used;
-   for (ix = 0; ix < pa; ix++) {
-      /* clear the carry */
-      u = 0;
-
-      /* left hand side of A[ix] * B[iy] */
-      tmpx = a->dp[ix];
-
-      /* alias to the address of where the digits will be stored */
-      tmpt = &(t.dp[digs]);
-
-      /* alias for where to read the right hand side from */
-      tmpy = b->dp + (digs - ix);
-
-      for (iy = digs - ix; iy < pb; iy++) {
-         /* calculate the double precision result */
-         r       = (mp_word)*tmpt +
-                   ((mp_word)tmpx * (mp_word)*tmpy++) +
-                   (mp_word)u;
-
-         /* get the lower part */
-         *tmpt++ = (mp_digit)(r & (mp_word)MP_MASK);
-
-         /* carry the carry */
-         u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
-      }
-      *tmpt = u;
-   }
-   mp_clamp(&t);
-   mp_exch(&t, c);
-   mp_clear(&t);
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_s_mp_mul_high_digs_fast.c b/libtommath/bn_s_mp_mul_high_digs_fast.c
deleted file mode 100644
index 0796f72..0000000
--- a/libtommath/bn_s_mp_mul_high_digs_fast.c
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_MUL_HIGH_DIGS_FAST_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* this is a modified version of s_mp_mul_digs_fast that only produces
- * output digits *above* digs.  See the comments for s_mp_mul_digs_fast
- * to see how it works.
- *
- * This is used in the Barrett reduction since for one of the multiplications
- * only the higher digits were needed.  This essentially halves the work.
- *
- * Based on Algorithm 14.12 on pp.595 of HAC.
- */
-mp_err s_mp_mul_high_digs_fast(const mp_int *a, const mp_int *b, mp_int *c, int digs)
-{
-   int     olduse, pa, ix, iz;
-   mp_err   err;
-   mp_digit W[MP_WARRAY];
-   mp_word  _W;
-
-   if (digs < 0) {
-      return MP_VAL;
-   }
-
-   /* grow the destination as required */
-   pa = a->used + b->used;
-   if (c->alloc < pa) {
-      if ((err = mp_grow(c, pa)) != MP_OKAY) {
-         return err;
-      }
-   }
-
-   /* number of output digits to produce */
-   pa = a->used + b->used;
-   _W = 0;
-   for (ix = digs; ix < pa; ix++) {
-      int      tx, ty, iy;
-      mp_digit *tmpx, *tmpy;
-
-      /* get offsets into the two bignums */
-      ty = MP_MIN(b->used-1, ix);
-      tx = ix - ty;
-
-      /* setup temp aliases */
-      tmpx = a->dp + tx;
-      tmpy = b->dp + ty;
-
-      /* this is the number of times the loop will iterrate, essentially its
-         while (tx++ < a->used && ty-- >= 0) { ... }
-       */
-      iy = MP_MIN(a->used-tx, ty+1);
-
-      /* execute loop */
-      for (iz = 0; iz < iy; iz++) {
-         _W += (mp_word)*tmpx++ * (mp_word)*tmpy--;
-      }
-
-      /* store term */
-      W[ix] = (mp_digit)_W & MP_MASK;
-
-      /* make next carry */
-      _W = _W >> (mp_word)MP_DIGIT_BIT;
-   }
-
-   /* setup dest */
-   olduse  = c->used;
-   c->used = pa;
-
-   {
-      mp_digit *tmpc;
-
-      tmpc = c->dp + digs;
-      for (ix = digs; ix < pa; ix++) {
-         /* now extract the previous digit [below the carry] */
-         *tmpc++ = W[ix];
-      }
-
-      /* clear unused digits [that existed in the old copy of c] */
-      MP_ZERO_DIGITS(tmpc, olduse - ix);
-   }
-   mp_clamp(c);
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_s_mp_prime_is_divisible.c b/libtommath/bn_s_mp_prime_is_divisible.c
deleted file mode 100644
index ffd5093..0000000
--- a/libtommath/bn_s_mp_prime_is_divisible.c
+++ /dev/null
@@ -1,35 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_PRIME_IS_DIVISIBLE_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* determines if an integers is divisible by one
- * of the first PRIME_SIZE primes or not
- *
- * sets result to 0 if not, 1 if yes
- */
-mp_err s_mp_prime_is_divisible(const mp_int *a, mp_bool *result)
-{
-   int      ix;
-   mp_err   err;
-   mp_digit res;
-
-   /* default to not */
-   *result = MP_NO;
-
-   for (ix = 0; ix < PRIVATE_MP_PRIME_TAB_SIZE; ix++) {
-      /* what is a mod LBL_prime_tab[ix] */
-      if ((err = mp_mod_d(a, s_mp_prime_tab[ix], &res)) != MP_OKAY) {
-         return err;
-      }
-
-      /* is the residue zero? */
-      if (res == 0u) {
-         *result = MP_YES;
-         return MP_OKAY;
-      }
-   }
-
-   return MP_OKAY;
-}
-#endif
diff --git a/libtommath/bn_s_mp_rand_jenkins.c b/libtommath/bn_s_mp_rand_jenkins.c
deleted file mode 100644
index c64afac..0000000
--- a/libtommath/bn_s_mp_rand_jenkins.c
+++ /dev/null
@@ -1,52 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_RAND_JENKINS_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* Bob Jenkins' http://burtleburtle.net/bob/rand/smallprng.html */
-/* Chosen for speed and a good "mix" */
-typedef struct {
-   uint64_t a;
-   uint64_t b;
-   uint64_t c;
-   uint64_t d;
-} ranctx;
-
-static ranctx jenkins_x;
-
-#define rot(x,k) (((x)<<(k))|((x)>>(64-(k))))
-static uint64_t s_rand_jenkins_val(void)
-{
-   uint64_t e = jenkins_x.a - rot(jenkins_x.b, 7);
-   jenkins_x.a = jenkins_x.b ^ rot(jenkins_x.c, 13);
-   jenkins_x.b = jenkins_x.c + rot(jenkins_x.d, 37);
-   jenkins_x.c = jenkins_x.d + e;
-   jenkins_x.d = e + jenkins_x.a;
-   return jenkins_x.d;
-}
-
-void s_mp_rand_jenkins_init(uint64_t seed)
-{
-   int i;
-   jenkins_x.a = 0xf1ea5eedULL;
-   jenkins_x.b = jenkins_x.c = jenkins_x.d = seed;
-   for (i = 0; i < 20; ++i) {
-      (void)s_rand_jenkins_val();
-   }
-}
-
-mp_err s_mp_rand_jenkins(void *p, size_t n)
-{
-   char *q = (char *)p;
-   while (n > 0u) {
-      int i;
-      uint64_t x = s_rand_jenkins_val();
-      for (i = 0; (i < 8) && (n > 0u); ++i, --n) {
-         *q++ = (char)(x & 0xFFuLL);
-         x >>= 8;
-      }
-   }
-   return MP_OKAY;
-}
-
-#endif
diff --git a/libtommath/bn_s_mp_rand_platform.c b/libtommath/bn_s_mp_rand_platform.c
deleted file mode 100644
index 27339bf..0000000
--- a/libtommath/bn_s_mp_rand_platform.c
+++ /dev/null
@@ -1,148 +0,0 @@
-#include "tommath_private.h"
-#ifdef BN_S_MP_RAND_PLATFORM_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* First the OS-specific special cases
- * - *BSD
- * - Windows
- */
-#if defined(__FreeBSD__) || defined(__OpenBSD__) || defined(__NetBSD__) || defined(__DragonFly__)
-#define BN_S_READ_ARC4RANDOM_C
-static mp_err s_read_arc4random(void *p, size_t n)
-{
-   arc4random_buf(p, n);
-   return MP_OKAY;
-}
-#endif
-
-#if defined(_WIN32) || defined(_WIN32_WCE)
-#define BN_S_READ_WINCSP_C
-
-#ifndef _WIN32_WINNT
-#define _WIN32_WINNT 0x0400
-#endif
-#ifdef _WIN32_WCE
-#define UNDER_CE
-#define ARM
-#endif
-
-#define WIN32_LEAN_AND_MEAN
-#include <windows.h>
-#include <wincrypt.h>
-
-static mp_err s_read_wincsp(void *p, size_t n)
-{
-   static HCRYPTPROV hProv = 0;
-   if (hProv == 0) {
-      HCRYPTPROV h = 0;
-      if (!CryptAcquireContext(&h, NULL, MS_DEF_PROV, PROV_RSA_FULL,
-                               (CRYPT_VERIFYCONTEXT | CRYPT_MACHINE_KEYSET)) &&
-          !CryptAcquireContext(&h, NULL, MS_DEF_PROV, PROV_RSA_FULL,
-                               CRYPT_VERIFYCONTEXT | CRYPT_MACHINE_KEYSET | CRYPT_NEWKEYSET)) {
-         return MP_ERR;
-      }
-      hProv = h;
-   }
-   return CryptGenRandom(hProv, (DWORD)n, (BYTE *)p) == TRUE ? MP_OKAY : MP_ERR;
-}
-#endif /* WIN32 */
-
-#if !defined(BN_S_READ_WINCSP_C) && defined(__linux__) && defined(__GLIBC_PREREQ)
-#if __GLIBC_PREREQ(2, 25)
-#define BN_S_READ_GETRANDOM_C
-#include <sys/random.h>
-#include <errno.h>
-
-static mp_err s_read_getrandom(void *p, size_t n)
-{
-   char *q = (char *)p;
-   while (n > 0u) {
-      ssize_t ret = getrandom(q, n, 0);
-      if (ret < 0) {
-         if (errno == EINTR) {
-            continue;
-         }
-         return MP_ERR;
-      }
-      q += ret;
-      n -= (size_t)ret;
-   }
-   return MP_OKAY;
-}
-#endif
-#endif
-
-/* We assume all platforms besides windows provide "/dev/urandom".
- * In case yours doesn't, define MP_NO_DEV_URANDOM at compile-time.
- */
-#if !defined(BN_S_READ_WINCSP_C) && !defined(MP_NO_DEV_URANDOM)
-#define BN_S_READ_URANDOM_C
-#ifndef MP_DEV_URANDOM
-#define MP_DEV_URANDOM "/dev/urandom"
-#endif
-#include <fcntl.h>
-#include <errno.h>
-#include <unistd.h>
-
-static mp_err s_read_urandom(void *p, size_t n)
-{
-   int fd;
-   char *q = (char *)p;
-
-   do {
-      fd = open(MP_DEV_URANDOM, O_RDONLY);
-   } while ((fd == -1) && (errno == EINTR));
-   if (fd == -1) return MP_ERR;
-
-   while (n > 0u) {
-      ssize_t ret = read(fd, p, n);
-      if (ret < 0) {
-         if (errno == EINTR) {
-            continue;
-         }
-         close(fd);
-         return MP_ERR;
-      }
-      q += ret;
-      n -= (size_t)ret;
-   }
-
-   close(fd);
-   return MP_OKAY;
-}
-#endif
-
-#if defined(MP_PRNG_ENABLE_LTM_RNG)
-#define BN_S_READ_LTM_RNG
-unsigned long (*ltm_rng)(unsigned char *out, unsigned long outlen, void (*callback)(void));
-void (*ltm_rng_callback)(void);
-
-static mp_err s_read_ltm_rng(void *p, size_t n)
-{
-   unsigned long res;
-   if (ltm_rng == NULL) return MP_ERR;
-   res = ltm_rng(p, n, ltm_rng_callback);
-   if (res != n) return MP_ERR;
-   return MP_OKAY;
-}
-#endif
-
-mp_err s_read_arc4random(void *p, size_t n);
-mp_err s_read_wincsp(void *p, size_t n);
-mp_err s_read_getrandom(void *p, size_t n);
-mp_err s_read_urandom(void *p, size_t n);
-mp_err s_read_ltm_rng(void *p, size_t n);
-
-mp_err s_mp_rand_platform(void *p, size_t n)
-{
-   mp_err err = MP_ERR;
-   if ((err != MP_OKAY) && MP_HAS(S_READ_ARC4RANDOM)) err = s_read_arc4random(p, n);
-   if ((err != MP_OKAY) && MP_HAS(S_READ_WINCSP))     err = s_read_wincsp(p, n);
-   if ((err != MP_OKAY) && MP_HAS(S_READ_GETRANDOM))  err = s_read_getrandom(p, n);
-   if ((err != MP_OKAY) && MP_HAS(S_READ_URANDOM))    err = s_read_urandom(p, n);
-   if ((err != MP_OKAY) && MP_HAS(S_READ_LTM_RNG))    err = s_read_ltm_rng(p, n);
-   return err;
-}
-
-#endif
diff --git a/libtommath/testme.sh b/libtommath/testme.sh
deleted file mode 100755
index 40fa32d..0000000
--- a/libtommath/testme.sh
+++ /dev/null
@@ -1,394 +0,0 @@
-#!/bin/bash
-#
-# return values of this script are:
-#   0  success
-# 128  a test failed
-#  >0  the number of timed-out tests
-# 255  parsing of parameters failed
-
-set -e
-
-if [ -f /proc/cpuinfo ]
-then
-  MAKE_JOBS=$(( ($(cat /proc/cpuinfo | grep -E '^processor[[:space:]]*:' | tail -n -1 | cut -d':' -f2) + 1) * 2 + 1 ))
-else
-  MAKE_JOBS=8
-fi
-
-ret=0
-TEST_CFLAGS=""
-
-_help()
-{
-  echo "Usage options for $(basename $0) [--with-cc=arg [other options]]"
-  echo
-  echo "Executing this script without any parameter will only run the default"
-  echo "configuration that has automatically been determined for the"
-  echo "architecture you're running."
-  echo
-  echo "    --with-cc=*             The compiler(s) to use for the tests"
-  echo "                            This is an option that will be iterated."
-  echo
-  echo "    --test-vs-mtest=*       Run test vs. mtest for '*' operations."
-  echo "                            Only the first of each options will be"
-  echo "                            taken into account."
-  echo
-  echo "To be able to specify options a compiler has to be given with"
-  echo "the option --with-cc=compilername"
-  echo "All other options will be tested with all MP_xBIT configurations."
-  echo
-  echo "    --with-{m64,m32,mx32}   The architecture(s) to build and test"
-  echo "                            for, e.g. --with-mx32."
-  echo "                            This is an option that will be iterated,"
-  echo "                            multiple selections are possible."
-  echo "                            The mx32 architecture is not supported"
-  echo "                            by clang and will not be executed."
-  echo
-  echo "    --cflags=*              Give an option to the compiler,"
-  echo "                            e.g. --cflags=-g"
-  echo "                            This is an option that will always be"
-  echo "                            passed as parameter to CC."
-  echo
-  echo "    --make-option=*         Give an option to make,"
-  echo "                            e.g. --make-option=\"-f makefile.shared\""
-  echo "                            This is an option that will always be"
-  echo "                            passed as parameter to make."
-  echo
-  echo "    --with-low-mp           Also build&run tests with -DMP_{8,16,32}BIT."
-  echo
-  echo "    --mtest-real-rand       Use real random data when running mtest."
-  echo
-  echo "    --with-valgrind"
-  echo "    --with-valgrind=*       Run in valgrind (slow!)."
-  echo
-  echo "    --with-travis-valgrind  Run with valgrind on Travis on specific branches."
-  echo
-  echo "    --valgrind-options      Additional Valgrind options"
-  echo "                            Some of the options like e.g.:"
-  echo "                            --track-origins=yes add a lot of extra"
-  echo "                            runtime and may trigger the 30 minutes"
-  echo "                            timeout."
-  echo
-  echo "Godmode:"
-  echo
-  echo "    --all                   Choose all architectures and gcc and clang"
-  echo "                            as compilers but does not run valgrind."
-  echo
-  echo "    --format                Runs the various source-code formatters"
-  echo "                            and generators and checks if the sources"
-  echo "                            are clean."
-  echo
-  echo "    -h"
-  echo "    --help                  This message"
-  echo
-  echo "    -v"
-  echo "    --version               Prints the version. It is just the number"
-  echo "                            of git commits to this file, no deeper"
-  echo "                            meaning attached"
-  exit 0
-}
-
-_die()
-{
-  echo "error $2 while $1"
-  if [ "$2" != "124" ]
-  then
-    exit 128
-  else
-    echo "assuming timeout while running test - continue"
-    local _tail=""
-    which tail >/dev/null && _tail="tail -n 1 test_${suffix}.log" && \
-    echo "last line of test_"${suffix}".log was:" && $_tail && echo ""
-    ret=$(( $ret + 1 ))
-  fi
-}
-
-_make()
-{
-  echo -ne " Compile $1 $2"
-  suffix=$(echo ${1}${2}  | tr ' ' '_')
-  CC="$1" CFLAGS="$2 $TEST_CFLAGS" make -j$MAKE_JOBS $3 $MAKE_OPTIONS > /dev/null 2>gcc_errors_${suffix}.log
-  errcnt=$(wc -l < gcc_errors_${suffix}.log)
-  if [[ ${errcnt} -gt 1 ]]; then
-    echo " failed"
-    cat gcc_errors_${suffix}.log
-    exit 128
-  fi
-}
-
-
-_runtest()
-{
-  make clean > /dev/null
-  local _timeout=""
-  which timeout >/dev/null && _timeout="timeout --foreground 90"
-  if [[ "$MAKE_OPTIONS" =~ "tune" ]]
-  then
-    # "make tune" will run "tune_it.sh" automatically, hence "autotune", but it cannot
-    # get switched off without some effort, so we just let it run twice for testing purposes
-    echo -e "\rRun autotune $1 $2"
-    _make "$1" "$2" ""
-    $_timeout $TUNE_CMD > test_${suffix}.log || _die "running autotune" $?
-  else
-    _make "$1" "$2" "test"
-    echo -e "\rRun test $1 $2"
-    $_timeout ./test > test_${suffix}.log || _die "running tests" $?
-  fi
-}
-
-# This is not much more of a C&P of _runtest with a different timeout
-# and the additional valgrind call.
-# TODO: merge
-_runvalgrind()
-{
-  make clean > /dev/null
-  local _timeout=""
-  # 30 minutes? Yes. Had it at 20 minutes and the Valgrind run needed over 25 minutes.
-  # A bit too close for comfort.
-  which timeout >/dev/null && _timeout="timeout --foreground 1800"
-echo "MAKE_OPTIONS = \"$MAKE_OPTIONS\""
-  if [[ "$MAKE_OPTIONS" =~ "tune"  ]]
-  then
-echo "autotune branch"
-    _make "$1" "$2" ""
-    # The shell used for /bin/sh is DASH 0.5.7-4ubuntu1 on the author's machine which fails valgrind, so
-    # we just run on instance of etc/tune with the same options as in etc/tune_it.sh
-    echo -e "\rRun etc/tune $1 $2 once inside valgrind"
-    $_timeout $VALGRIND_BIN $VALGRIND_OPTS $TUNE_CMD > test_${suffix}.log || _die "running etc/tune" $?
-  else
-    _make "$1" "$2" "test"
-    echo -e "\rRun test $1 $2 inside valgrind"
-    $_timeout $VALGRIND_BIN $VALGRIND_OPTS ./test > test_${suffix}.log || _die "running tests" $?
-  fi
-}
-
-
-_banner()
-{
-  echo "uname="$(uname -a)
-  [[ "$#" != "0" ]] && (echo $1=$($1 -dumpversion)) || true
-}
-
-_exit()
-{
-  if [ "$ret" == "0" ]
-  then
-    echo "Tests successful"
-  else
-    echo "$ret tests timed out"
-  fi
-
-  exit $ret
-}
-
-ARCHFLAGS=""
-COMPILERS=""
-CFLAGS=""
-WITH_LOW_MP=""
-TEST_VS_MTEST=""
-MTEST_RAND=""
-# timed with an AMD A8-6600K
-# 25 minutes
-#VALGRIND_OPTS=" --track-origins=yes --leak-check=full --show-leak-kinds=all --error-exitcode=1 "
-# 9 minutes (14 minutes with --test-vs-mtest=333333 --mtest-real-rand)
-VALGRIND_OPTS=" --leak-check=full --show-leak-kinds=all --error-exitcode=1 "
-#VALGRIND_OPTS=""
-VALGRIND_BIN=""
-CHECK_FORMAT=""
-TUNE_CMD="./etc/tune -t -r 10 -L 3"
-
-alive_pid=0
-
-function kill_alive() {
-  disown $alive_pid || true
-  kill $alive_pid 2>/dev/null
-}
-
-function start_alive_printing() {
-  [ "$alive_pid" == "0" ] || return 0;
-  for i in `seq 1 10` ; do sleep 300 && echo "Tests still in Progress..."; done &
-  alive_pid=$!
-  trap kill_alive EXIT
-}
-
-while [ $# -gt 0 ];
-do
-  case $1 in
-    "--with-m64" | "--with-m32" | "--with-mx32")
-      ARCHFLAGS="$ARCHFLAGS ${1:6}"
-    ;;
-    --with-cc=*)
-      COMPILERS="$COMPILERS ${1#*=}"
-    ;;
-    --cflags=*)
-      CFLAGS="$CFLAGS ${1#*=}"
-    ;;
-    --valgrind-options=*)
-      VALGRIND_OPTS="$VALGRIND_OPTS ${1#*=}"
-    ;;
-    --with-valgrind*)
-      if [[ ${1#*d} != "" ]]
-      then
-        VALGRIND_BIN="${1#*=}"
-      else
-        VALGRIND_BIN="valgrind"
-      fi
-      start_alive_printing
-    ;;
-    --with-travis-valgrind*)
-      if [[ ("$TRAVIS_BRANCH" == "develop" && "$TRAVIS_PULL_REQUEST" == "false") || "$TRAVIS_BRANCH" == *"valgrind"* || "$TRAVIS_COMMIT_MESSAGE" == *"valgrind"* ]]
-      then
-        if [[ ${1#*d} != "" ]]
-        then
-          VALGRIND_BIN="${1#*=}"
-        else
-          VALGRIND_BIN="valgrind"
-        fi
-        start_alive_printing
-      fi
-    ;;
-    --make-option=*)
-      MAKE_OPTIONS="$MAKE_OPTIONS ${1#*=}"
-    ;;
-    --with-low-mp)
-      WITH_LOW_MP="1"
-    ;;
-    --test-vs-mtest=*)
-      TEST_VS_MTEST="${1#*=}"
-      if ! [ "$TEST_VS_MTEST" -eq "$TEST_VS_MTEST" ] 2> /dev/null
-      then
-         echo "--test-vs-mtest Parameter has to be int"
-         exit 255
-      fi
-      start_alive_printing
-    ;;
-    --mtest-real-rand)
-      MTEST_RAND="-DLTM_MTEST_REAL_RAND"
-    ;;
-    --format)
-      CHECK_FORMAT="1"
-    ;;
-    --all)
-      COMPILERS="gcc clang"
-      ARCHFLAGS="-m64 -m32 -mx32"
-    ;;
-    --help | -h)
-      _help
-    ;;
-    --version | -v)
-      echo $(git rev-list HEAD --count -- testme.sh) || echo "Unknown. Please run in original libtommath git repository."
-      exit 0
-    ;;
-    *)
-      echo "Ignoring option ${1}"
-    ;;
-  esac
-  shift
-done
-
-function _check_git() {
-  git update-index --refresh >/dev/null || true
-  git diff-index --quiet HEAD -- . || ( echo "FAILURE: $*" && exit 1 )
-}
-
-if [[ "$CHECK_FORMAT" == "1" ]]
-then
-  make astyle
-  _check_git "make astyle"
-  perl helper.pl --update-files
-  _check_git "helper.pl --update-files"
-  perl helper.pl --check-all
-  _check_git "helper.pl --check-all"
-  exit $?
-fi
-
-[[ "$VALGRIND_BIN" == "" ]] && VALGRIND_OPTS=""
-
-# default to CC environment variable if no compiler is defined but some other options
-if [[ "$COMPILERS" == "" ]] && [[ "$ARCHFLAGS$MAKE_OPTIONS$CFLAGS" != "" ]]
-then
-   COMPILERS="$CC"
-# default to CC environment variable and run only default config if no option is given
-elif [[ "$COMPILERS" == "" ]]
-then
-  _banner "$CC"
-  if [[ "$VALGRIND_BIN" != "" ]]
-  then
-    _runvalgrind "$CC" ""
-  else
-    _runtest "$CC" ""
-  fi
-  _exit
-fi
-
-
-archflags=( $ARCHFLAGS )
-compilers=( $COMPILERS )
-
-# choosing a compiler without specifying an architecture will use the default architecture
-if [ "${#archflags[@]}" == "0" ]
-then
-  archflags[0]=" "
-fi
-
-_banner
-
-if [[ "$TEST_VS_MTEST" != "" ]]
-then
-   make clean > /dev/null
-   _make "${compilers[0]} ${archflags[0]}" "$CFLAGS" "mtest_opponent"
-   echo
-   _make "gcc" "$MTEST_RAND" "mtest"
-   echo
-   echo "Run test vs. mtest for $TEST_VS_MTEST iterations"
-   _timeout=""
-   which timeout >/dev/null && _timeout="timeout --foreground 1800"
-   $_timeout ./mtest/mtest $TEST_VS_MTEST | $VALGRIND_BIN $VALGRIND_OPTS  ./mtest_opponent > valgrind_test.log 2> test_vs_mtest_err.log
-   retval=$?
-   head -n 5 valgrind_test.log
-   tail -n 2 valgrind_test.log
-   exit $retval
-fi
-
-for i in "${compilers[@]}"
-do
-  if [ -z "$(which $i)" ]
-  then
-    echo "Skipped compiler $i, file not found"
-    continue
-  fi
-  compiler_version=$(echo "$i="$($i -dumpversion))
-  if [ "$compiler_version" == "clang=4.2.1" ]
-  then
-    # one of my versions of clang complains about some stuff in stdio.h and stdarg.h ...
-    TEST_CFLAGS="-Wno-typedef-redefinition"
-  else
-    TEST_CFLAGS=""
-  fi
-  echo $compiler_version
-
-  for a in "${archflags[@]}"
-  do
-    if [[ $(expr "$i" : "clang") -ne 0 && "$a" == "-mx32" ]]
-    then
-      echo "clang -mx32 tests skipped"
-      continue
-    fi
-    if [[ "$VALGRIND_BIN" != "" ]]
-    then
-      _runvalgrind "$i $a" "$CFLAGS"
-      [ "$WITH_LOW_MP" != "1" ] && continue
-      _runvalgrind "$i $a" "-DMP_8BIT $CFLAGS"
-      _runvalgrind "$i $a" "-DMP_16BIT $CFLAGS"
-      _runvalgrind "$i $a" "-DMP_32BIT $CFLAGS"
-    else
-      _runtest "$i $a" "$CFLAGS"
-      [ "$WITH_LOW_MP" != "1" ] && continue
-      _runtest "$i $a" "-DMP_8BIT $CFLAGS"
-      _runtest "$i $a" "-DMP_16BIT $CFLAGS"
-      _runtest "$i $a" "-DMP_32BIT $CFLAGS"
-    fi
-  done
-done
-
-_exit
diff --git a/unix/Makefile.in b/unix/Makefile.in
index ce07a6a..2ebadbb 100644
--- a/unix/Makefile.in
+++ b/unix/Makefile.in
@@ -489,13 +489,8 @@ STUB_SRCS = \
 	$(GENERIC_DIR)/tclOOStubLib.c
 
 TOMMATH_SRCS = \
-	$(TOMMATH_DIR)/bn_cutoffs.c \
-	$(TOMMATH_DIR)/bn_deprecated.c \
-	$(TOMMATH_DIR)/bn_mp_2expt.c \
-	$(TOMMATH_DIR)/bn_mp_abs.c \
 	$(TOMMATH_DIR)/bn_mp_add.c \
 	$(TOMMATH_DIR)/bn_mp_add_d.c \
-	$(TOMMATH_DIR)/bn_mp_addmod.c \
 	$(TOMMATH_DIR)/bn_mp_and.c \
 	$(TOMMATH_DIR)/bn_mp_clamp.c \
 	$(TOMMATH_DIR)/bn_mp_clear.c \
@@ -504,136 +499,55 @@ TOMMATH_SRCS = \
 	$(TOMMATH_DIR)/bn_mp_cmp_d.c \
 	$(TOMMATH_DIR)/bn_mp_cmp_mag.c \
 	$(TOMMATH_DIR)/bn_mp_cnt_lsb.c \
-	$(TOMMATH_DIR)/bn_mp_complement.c \
 	$(TOMMATH_DIR)/bn_mp_copy.c \
 	$(TOMMATH_DIR)/bn_mp_count_bits.c \
-	$(TOMMATH_DIR)/bn_mp_decr.c \
 	$(TOMMATH_DIR)/bn_mp_div.c \
 	$(TOMMATH_DIR)/bn_mp_div_2.c \
 	$(TOMMATH_DIR)/bn_mp_div_2d.c \
 	$(TOMMATH_DIR)/bn_mp_div_3.c \
 	$(TOMMATH_DIR)/bn_mp_div_d.c \
-	$(TOMMATH_DIR)/bn_mp_dr_is_modulus.c \
-	$(TOMMATH_DIR)/bn_mp_dr_reduce.c \
-	$(TOMMATH_DIR)/bn_mp_dr_setup.c \
-	$(TOMMATH_DIR)/bn_mp_error_to_string.c \
 	$(TOMMATH_DIR)/bn_mp_exch.c \
 	$(TOMMATH_DIR)/bn_mp_expt_u32.c \
-	$(TOMMATH_DIR)/bn_mp_exptmod.c \
-	$(TOMMATH_DIR)/bn_mp_exteuclid.c \
-	$(TOMMATH_DIR)/bn_mp_fread.c \
-	$(TOMMATH_DIR)/bn_mp_from_sbin.c \
-	$(TOMMATH_DIR)/bn_mp_from_ubin.c \
-	$(TOMMATH_DIR)/bn_mp_fwrite.c \
-	$(TOMMATH_DIR)/bn_mp_gcd.c \
-	$(TOMMATH_DIR)/bn_mp_get_double.c \
-	$(TOMMATH_DIR)/bn_mp_get_i32.c \
-	$(TOMMATH_DIR)/bn_mp_get_i64.c \
-	$(TOMMATH_DIR)/bn_mp_get_l.c \
-	$(TOMMATH_DIR)/bn_mp_get_mag_u32.c \
-	$(TOMMATH_DIR)/bn_mp_get_mag_u64.c \
-	$(TOMMATH_DIR)/bn_mp_get_mag_ul.c \
 	$(TOMMATH_DIR)/bn_mp_grow.c \
-	$(TOMMATH_DIR)/bn_mp_incr.c \
 	$(TOMMATH_DIR)/bn_mp_init.c \
 	$(TOMMATH_DIR)/bn_mp_init_copy.c \
-	$(TOMMATH_DIR)/bn_mp_init_i32.c \
-	$(TOMMATH_DIR)/bn_mp_init_i64.c \
-	$(TOMMATH_DIR)/bn_mp_init_l.c \
 	$(TOMMATH_DIR)/bn_mp_init_multi.c \
 	$(TOMMATH_DIR)/bn_mp_init_set.c \
 	$(TOMMATH_DIR)/bn_mp_init_size.c \
-	$(TOMMATH_DIR)/bn_mp_init_u32.c \
-	$(TOMMATH_DIR)/bn_mp_init_u64.c \
-	$(TOMMATH_DIR)/bn_mp_init_ul.c \
-	$(TOMMATH_DIR)/bn_mp_invmod.c \
-	$(TOMMATH_DIR)/bn_mp_is_square.c \
-	$(TOMMATH_DIR)/bn_mp_iseven.c \
-	$(TOMMATH_DIR)/bn_mp_isodd.c \
-	$(TOMMATH_DIR)/bn_mp_kronecker.c \
-	$(TOMMATH_DIR)/bn_mp_lcm.c \
-	$(TOMMATH_DIR)/bn_mp_log_u32.c \
 	$(TOMMATH_DIR)/bn_mp_lshd.c \
 	$(TOMMATH_DIR)/bn_mp_mod.c \
 	$(TOMMATH_DIR)/bn_mp_mod_2d.c \
-	$(TOMMATH_DIR)/bn_mp_mod_d.c \
-	$(TOMMATH_DIR)/bn_mp_montgomery_calc_normalization.c \
-	$(TOMMATH_DIR)/bn_mp_montgomery_reduce.c \
-	$(TOMMATH_DIR)/bn_mp_montgomery_setup.c \
 	$(TOMMATH_DIR)/bn_mp_mul.c \
 	$(TOMMATH_DIR)/bn_mp_mul_2.c \
 	$(TOMMATH_DIR)/bn_mp_mul_2d.c \
 	$(TOMMATH_DIR)/bn_mp_mul_d.c \
-	$(TOMMATH_DIR)/bn_mp_mulmod.c \
 	$(TOMMATH_DIR)/bn_mp_neg.c \
 	$(TOMMATH_DIR)/bn_mp_or.c \
 	$(TOMMATH_DIR)/bn_mp_pack.c \
 	$(TOMMATH_DIR)/bn_mp_pack_count.c \
-	$(TOMMATH_DIR)/bn_mp_prime_fermat.c \
-	$(TOMMATH_DIR)/bn_mp_prime_frobenius_underwood.c \
-	$(TOMMATH_DIR)/bn_mp_prime_is_prime.c \
-	$(TOMMATH_DIR)/bn_mp_prime_miller_rabin.c \
-	$(TOMMATH_DIR)/bn_mp_prime_next_prime.c \
-	$(TOMMATH_DIR)/bn_mp_prime_rabin_miller_trials.c \
-	$(TOMMATH_DIR)/bn_mp_prime_rand.c \
-	$(TOMMATH_DIR)/bn_mp_prime_strong_lucas_selfridge.c \
 	$(TOMMATH_DIR)/bn_mp_radix_size.c \
 	$(TOMMATH_DIR)/bn_mp_radix_smap.c \
-	$(TOMMATH_DIR)/bn_mp_rand.c \
 	$(TOMMATH_DIR)/bn_mp_read_radix.c \
-	$(TOMMATH_DIR)/bn_mp_reduce.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_2k.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_2k_l.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_2k_setup.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_2k_setup_l.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_is_2k.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_is_2k_l.c \
-	$(TOMMATH_DIR)/bn_mp_reduce_setup.c \
-	$(TOMMATH_DIR)/bn_mp_root_u32.c \
 	$(TOMMATH_DIR)/bn_mp_rshd.c \
-	$(TOMMATH_DIR)/bn_mp_sbin_size.c \
 	$(TOMMATH_DIR)/bn_mp_set.c \
-	$(TOMMATH_DIR)/bn_mp_set_double.c \
-	$(TOMMATH_DIR)/bn_mp_set_i32.c \
-	$(TOMMATH_DIR)/bn_mp_set_i64.c \
-	$(TOMMATH_DIR)/bn_mp_set_l.c \
-	$(TOMMATH_DIR)/bn_mp_set_u32.c \
-	$(TOMMATH_DIR)/bn_mp_set_u64.c \
-	$(TOMMATH_DIR)/bn_mp_set_ul.c \
 	$(TOMMATH_DIR)/bn_mp_shrink.c \
 	$(TOMMATH_DIR)/bn_mp_signed_rsh.c \
 	$(TOMMATH_DIR)/bn_mp_sqr.c \
-	$(TOMMATH_DIR)/bn_mp_sqrmod.c \
 	$(TOMMATH_DIR)/bn_mp_sqrt.c \
-	$(TOMMATH_DIR)/bn_mp_sqrtmod_prime.c \
 	$(TOMMATH_DIR)/bn_mp_sub.c \
 	$(TOMMATH_DIR)/bn_mp_sub_d.c \
-	$(TOMMATH_DIR)/bn_mp_submod.c \
 	$(TOMMATH_DIR)/bn_mp_to_radix.c \
-	$(TOMMATH_DIR)/bn_mp_to_sbin.c \
 	$(TOMMATH_DIR)/bn_mp_to_ubin.c \
 	$(TOMMATH_DIR)/bn_mp_ubin_size.c \
 	$(TOMMATH_DIR)/bn_mp_unpack.c \
 	$(TOMMATH_DIR)/bn_mp_xor.c \
 	$(TOMMATH_DIR)/bn_mp_zero.c \
-	$(TOMMATH_DIR)/bn_prime_tab.c \
 	$(TOMMATH_DIR)/bn_s_mp_add.c \
 	$(TOMMATH_DIR)/bn_s_mp_balance_mul.c \
-	$(TOMMATH_DIR)/bn_s_mp_exptmod.c \
-	$(TOMMATH_DIR)/bn_s_mp_exptmod_fast.c \
-	$(TOMMATH_DIR)/bn_s_mp_get_bit.c \
-	$(TOMMATH_DIR)/bn_s_mp_invmod_fast.c \
-	$(TOMMATH_DIR)/bn_s_mp_invmod_slow.c \
 	$(TOMMATH_DIR)/bn_s_mp_karatsuba_mul.c \
 	$(TOMMATH_DIR)/bn_s_mp_karatsuba_sqr.c \
-	$(TOMMATH_DIR)/bn_s_mp_montgomery_reduce_fast.c \
 	$(TOMMATH_DIR)/bn_s_mp_mul_digs.c \
 	$(TOMMATH_DIR)/bn_s_mp_mul_digs_fast.c \
-	$(TOMMATH_DIR)/bn_s_mp_mul_high_digs.c \
-	$(TOMMATH_DIR)/bn_s_mp_mul_high_digs_fast.c \
-	$(TOMMATH_DIR)/bn_s_mp_prime_is_divisible.c \
-	$(TOMMATH_DIR)/bn_s_mp_rand_jenkins.c \
-	$(TOMMATH_DIR)/bn_s_mp_rand_platform.c \
 	$(TOMMATH_DIR)/bn_s_mp_reverse.c \
 	$(TOMMATH_DIR)/bn_s_mp_sqr.c \
 	$(TOMMATH_DIR)/bn_s_mp_sqr_fast.c \