diff --git a/bn_mp_factors_add.c b/bn_mp_factors_add.c
new file mode 100644
index 000000000..2dd5c5f82
--- /dev/null
+++ b/bn_mp_factors_add.c
@@ -0,0 +1,30 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_ADD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Add one element (mp_int) to a factor-array (an array of mp_int's) */
+int mp_factors_add(const mp_int *a, mp_factors *f)
+{
+   int new_size = f->alloc;
+   mp_int *tmp;
+   mp_int t;
+   if (f->alloc <= f->length) {
+      new_size += LTM_TRIAL_GROWTH;
+      tmp = (mp_int *) XREALLOC(f->factors, sizeof(*tmp) * f->alloc,
+                                sizeof(*tmp) * (size_t)(new_size));
+      if (tmp == NULL) {
+         return MP_MEM;
+      }
+      f->factors = tmp;
+      f->alloc = new_size;
+   }
+   mp_init(&t);
+   mp_copy(a,&t);
+   f->factors[f->length] = t;
+   f->length++;
+   return MP_OKAY;
+}
+
+#endif
+
diff --git a/bn_mp_factors_clear.c b/bn_mp_factors_clear.c
new file mode 100644
index 000000000..c68c00a89
--- /dev/null
+++ b/bn_mp_factors_clear.c
@@ -0,0 +1,22 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_CLEAR_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Clear and free factor-array (an array of mp_int's) */
+void mp_factors_clear(mp_factors *f)
+{
+   int i;
+   if (f->factors != NULL) {
+      for (i = 0; i < f->length; i++) {
+         mp_clear(&f->factors[i]);
+      }
+      XFREE(f->factors, sizeof(mp_int) * f->alloc);
+      f->factors = NULL;
+      f->alloc = 0;
+      f->length = 0;
+   }
+}
+
+#endif
+
diff --git a/bn_mp_factors_compress.c b/bn_mp_factors_compress.c
new file mode 100644
index 000000000..c0f99ca45
--- /dev/null
+++ b/bn_mp_factors_compress.c
@@ -0,0 +1,68 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_COMPRESS_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+
+/* Change (sorted) factors array into format prime,prime-count,...,prime,prime-count */
+int mp_factors_compress(mp_factors *f, int sorted, mp_factors *g)
+{
+   int i, e = MP_OKAY;
+   mp_int count;
+   mp_factors in, out, t;
+
+   if ((e = mp_init(&count)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+   if ((e = mp_factors_init(&out)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+   if ((e = mp_factors_init(&in)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+   for (i = 0; i < f->length; i++) {
+      if ((e = mp_factors_add(&(f->factors[i]), &in)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+   }
+   if (sorted == MP_NO) {
+      mp_factors_sort(&in);
+   }
+
+   for (i = 1; i < in.length; i++) {
+      if ((e = mp_incr(&count)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (mp_cmp(&(in.factors[i-1]), &(in.factors[i])) != MP_EQ) {
+         if ((e = mp_factors_add(&(in.factors[i-1]), &out)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+         if ((e = mp_factors_add(&count, &out)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+         mp_zero(&count);
+      }
+   }
+   /* and the last one, too */
+   if ((e = mp_incr(&count)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+   if ((e = mp_factors_add(&(in.factors[i-1]), &out)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+   if ((e = mp_factors_add(&count, &out)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+   t = *g;
+   *g = out;
+   out = t;
+
+LTM_ERR:
+   mp_factors_clear(&in);
+   mp_clear(&count);
+   mp_factors_clear(&out);
+   return e;
+}
+#endif
diff --git a/bn_mp_factors_init.c b/bn_mp_factors_init.c
new file mode 100644
index 000000000..f180d586d
--- /dev/null
+++ b/bn_mp_factors_init.c
@@ -0,0 +1,19 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_INIT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Init factor-array (an array of mp_int's) */
+int mp_factors_init(mp_factors *f)
+{
+   f->factors = (mp_int *) XCALLOC((size_t)(LTM_TRIAL_GROWTH), sizeof(mp_int));
+   if ((f->factors) == NULL) {
+      return MP_MEM;
+   }
+   f->alloc = LTM_TRIAL_GROWTH;
+   f->length = 0;
+   return MP_OKAY;
+}
+
+#endif
+
diff --git a/bn_mp_factors_print.c b/bn_mp_factors_print.c
new file mode 100644
index 000000000..ec711e05f
--- /dev/null
+++ b/bn_mp_factors_print.c
@@ -0,0 +1,32 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_PRINT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Print the elememts of a factor-array (an array of mp_int's) */
+int mp_factors_print(const mp_factors *f, int base, char delimiter, FILE *stream)
+{
+   int i, e = MP_OKAY;
+   int d;
+   if (delimiter != 0) {
+      d = (int)delimiter;
+   } else {
+      d = ',';
+   }
+   if (f->factors != NULL) {
+      for (i = 0; i < f->length; i++) {
+         if ((e = mp_fwrite(&(f->factors[i]),base, stream)) != MP_OKAY) {
+            return e;
+         }
+         if (i < ((f->length)-1)) {
+            if (fputc(d, stream) == EOF) {
+               return MP_VAL;
+            }
+         }
+      }
+   }
+   return e;
+}
+
+#endif
+
diff --git a/bn_mp_factors_product.c b/bn_mp_factors_product.c
new file mode 100644
index 000000000..8cfb4ee3a
--- /dev/null
+++ b/bn_mp_factors_product.c
@@ -0,0 +1,63 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_PRODUCT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Compute the product of the elements of "factors" */
+static int s_product_lowlevel(const mp_factors *factors, int start, int n, mp_int *result)
+{
+   int i, first_half, second_half;
+   mp_int tmp;
+   int e = MP_OKAY;
+
+   if (n == 0) {
+      mp_set(result, 1uL);
+      return MP_OKAY;
+   }
+   /* Do the rest linearily. Faster for primorials at least,  but YMMV */
+   if (n <= 64) {
+      mp_set_long(result,1uL);
+      for (i = 0; i < n; i++) {
+         if ((e = mp_mul(result, &(factors->factors[i+start]), result)) != MP_OKAY) {
+            return e;
+         }
+      }
+      return MP_OKAY;
+   }
+
+   first_half = n / 2;
+   second_half = n - first_half;
+   if ((e = s_product_lowlevel(factors, start, second_half, result)) != MP_OKAY) {
+      return e;
+   }
+   if ((e = mp_init(&tmp)) != MP_OKAY) {
+      return e;
+   }
+   if ((e = s_product_lowlevel(factors, start+second_half, first_half, &tmp)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+   if ((e = mp_mul(result, &tmp, result)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+LTM_ERR:
+   mp_clear(&tmp);
+   return e;
+}
+
+int mp_factors_product(const mp_factors *factors, mp_int *p)
+{
+   int e, r;
+   r = factors->length;
+   if (r == 0) {
+      return MP_VAL;
+   }
+   if (r == 1) {
+      return mp_copy(&(factors->factors[0]), p);
+   }
+   if ((e = s_product_lowlevel(factors, 0, r, p)) != MP_OKAY) {
+      return e;
+   }
+   return MP_OKAY;
+}
+#endif
diff --git a/bn_mp_factors_sort.c b/bn_mp_factors_sort.c
new file mode 100644
index 000000000..fb4a6511b
--- /dev/null
+++ b/bn_mp_factors_sort.c
@@ -0,0 +1,27 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_SORT_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Sort factor-array (an array of mp_int's) */
+/* TODO: Sorts inline. Work on copy instead? */
+void mp_factors_sort(mp_factors *f)
+{
+   int i, idx;
+   mp_int tmp;
+
+   /*
+      It will be almost always be already sorted and even if it is unsorted
+      it will just show a bushy tail. So insertion sort it is.
+   */
+   for (i = 1 ; i < f->length; i++) {
+      idx = i;
+      while ((idx > 0) && (mp_cmp(&(f->factors[idx-1]),&(f->factors[idx])) == MP_GT)) {
+         mp_exch(&(f->factors[idx]), &tmp);
+         mp_exch(&(f->factors[idx-1]), &(f->factors[idx]));
+         mp_exch(&tmp, &(f->factors[idx-1]));
+         idx--;
+      }
+   }
+}
+#endif
diff --git a/bn_mp_factors_zero.c b/bn_mp_factors_zero.c
new file mode 100644
index 000000000..784e0ebd5
--- /dev/null
+++ b/bn_mp_factors_zero.c
@@ -0,0 +1,33 @@
+#include "tommath_private.h"
+#ifdef BN_MP_FACTORS_ZERO_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* Zero the elements of a factor-array (an array of mp_int's) and realloc to default
+   size in that order */
+int mp_factors_zero(mp_factors *f)
+{
+   int i, e = MP_OKAY;
+   mp_int *tmp;
+
+   if (f->factors != NULL) {
+      for (i = 0; i < f->length; i++) {
+         mp_clear(&(f->factors[i]));
+      }
+   }
+
+   tmp = (mp_int *) XREALLOC(f->factors, sizeof(*tmp) * f->alloc,
+                             sizeof(*tmp) * LTM_TRIAL_GROWTH);
+   if (tmp == NULL) {
+      return MP_MEM;
+   }
+
+   f->factors = tmp;
+   f->alloc = LTM_TRIAL_GROWTH;
+   f->length = 0;
+
+   return e;
+}
+
+#endif
+
diff --git a/callgraph.txt b/callgraph.txt
index eeeea1383..42e672814 100644
--- a/callgraph.txt
+++ b/callgraph.txt
@@ -2496,6 +2496,226 @@ BN_MP_EXTEUCLID_C
 |   +--->BN_MP_CLEAR_C
 
 
+BN_MP_FACTORS_ADD_C
++--->BN_MP_INIT_C
++--->BN_MP_COPY_C
+|   +--->BN_MP_GROW_C
+
+
+BN_MP_FACTORS_CLEAR_C
++--->BN_MP_CLEAR_C
+
+
+BN_MP_FACTORS_COMPRESS_C
++--->BN_MP_INIT_C
++--->BN_MP_FACTORS_INIT_C
++--->BN_MP_FACTORS_ADD_C
+|   +--->BN_MP_COPY_C
+|   |   +--->BN_MP_GROW_C
++--->BN_MP_FACTORS_SORT_C
+|   +--->BN_MP_CMP_C
+|   |   +--->BN_MP_CMP_MAG_C
+|   +--->BN_MP_EXCH_C
++--->BN_MP_INCR_C
+|   +--->BN_MP_SET_C
+|   |   +--->BN_MP_ZERO_C
+|   +--->BN_MP_DECR_C
+|   |   +--->BN_MP_ZERO_C
+|   |   +--->BN_MP_SUB_D_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_ADD_D_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   +--->BN_MP_ADD_D_C
+|   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_SUB_D_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_CLAMP_C
++--->BN_MP_CMP_C
+|   +--->BN_MP_CMP_MAG_C
++--->BN_MP_ZERO_C
++--->BN_MP_FACTORS_CLEAR_C
+|   +--->BN_MP_CLEAR_C
++--->BN_MP_CLEAR_C
+
+
+BN_MP_FACTORS_INIT_C
+
+
+BN_MP_FACTORS_PRINT_C
++--->BN_MP_FWRITE_C
+|   +--->BN_MP_RADIX_SIZE_C
+|   |   +--->BN_MP_COUNT_BITS_C
+|   |   +--->BN_MP_INIT_COPY_C
+|   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_COPY_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_DIV_D_C
+|   |   |   +--->BN_MP_COPY_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_DIV_2D_C
+|   |   |   |   +--->BN_MP_ZERO_C
+|   |   |   |   +--->BN_MP_MOD_2D_C
+|   |   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   |   +--->BN_MP_RSHD_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_DIV_3_C
+|   |   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   |   +--->BN_MP_EXCH_C
+|   |   |   |   +--->BN_MP_CLEAR_C
+|   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_EXCH_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_CLEAR_C
+|   +--->BN_MP_TORADIX_C
+|   |   +--->BN_MP_INIT_COPY_C
+|   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_COPY_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_DIV_D_C
+|   |   |   +--->BN_MP_COPY_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_DIV_2D_C
+|   |   |   |   +--->BN_MP_ZERO_C
+|   |   |   |   +--->BN_MP_MOD_2D_C
+|   |   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   |   +--->BN_MP_RSHD_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_DIV_3_C
+|   |   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   |   +--->BN_MP_EXCH_C
+|   |   |   |   +--->BN_MP_CLEAR_C
+|   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_EXCH_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_CLEAR_C
+
+
+BN_MP_FACTORS_PRODUCT_C
++--->BN_MP_SET_C
+|   +--->BN_MP_ZERO_C
++--->BN_MP_SET_LONG_C
++--->BN_MP_MUL_C
+|   +--->BN_FAST_S_MP_MUL_DIGS_C
+|   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_CLAMP_C
+|   +--->BN_MP_BALANCE_MUL_C
+|   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_INIT_C
+|   |   +--->BN_MP_INIT_MULTI_C
+|   |   |   +--->BN_MP_INIT_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_LSHD_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_RSHD_C
+|   |   |   |   +--->BN_MP_ZERO_C
+|   |   +--->BN_MP_ADD_C
+|   |   |   +--->BN_S_MP_ADD_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_CMP_MAG_C
+|   |   |   +--->BN_S_MP_SUB_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_EXCH_C
+|   |   +--->BN_MP_CLEAR_MULTI_C
+|   +--->BN_MP_TOOM_MUL_C
+|   |   +--->BN_MP_INIT_MULTI_C
+|   |   |   +--->BN_MP_INIT_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_MOD_2D_C
+|   |   |   +--->BN_MP_ZERO_C
+|   |   |   +--->BN_MP_COPY_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_COPY_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_RSHD_C
+|   |   |   +--->BN_MP_ZERO_C
+|   |   +--->BN_MP_MUL_2_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_ADD_C
+|   |   |   +--->BN_S_MP_ADD_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_CMP_MAG_C
+|   |   |   +--->BN_S_MP_SUB_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_SUB_C
+|   |   |   +--->BN_S_MP_ADD_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_CMP_MAG_C
+|   |   |   +--->BN_S_MP_SUB_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_DIV_2_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_MUL_2D_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_LSHD_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_MUL_D_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_DIV_3_C
+|   |   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   |   +--->BN_MP_INIT_C
+|   |   |   +--->BN_MP_CLAMP_C
+|   |   |   +--->BN_MP_EXCH_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   |   +--->BN_MP_LSHD_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_CLEAR_MULTI_C
+|   |   |   +--->BN_MP_CLEAR_C
+|   +--->BN_MP_KARATSUBA_MUL_C
+|   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_INIT_C
+|   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_S_MP_ADD_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_ADD_C
+|   |   |   +--->BN_MP_CMP_MAG_C
+|   |   |   +--->BN_S_MP_SUB_C
+|   |   |   |   +--->BN_MP_GROW_C
+|   |   +--->BN_S_MP_SUB_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   +--->BN_MP_LSHD_C
+|   |   |   +--->BN_MP_GROW_C
+|   |   |   +--->BN_MP_RSHD_C
+|   |   |   |   +--->BN_MP_ZERO_C
+|   |   +--->BN_MP_CLEAR_C
+|   +--->BN_S_MP_MUL_DIGS_C
+|   |   +--->BN_MP_INIT_SIZE_C
+|   |   |   +--->BN_MP_INIT_C
+|   |   +--->BN_MP_CLAMP_C
+|   |   +--->BN_MP_EXCH_C
+|   |   +--->BN_MP_CLEAR_C
++--->BN_MP_INIT_C
++--->BN_MP_CLEAR_C
++--->BN_MP_COPY_C
+|   +--->BN_MP_GROW_C
+
+
+BN_MP_FACTORS_SORT_C
++--->BN_MP_CMP_C
+|   +--->BN_MP_CMP_MAG_C
++--->BN_MP_EXCH_C
+
+
+BN_MP_FACTORS_ZERO_C
++--->BN_MP_CLEAR_C
+
+
 BN_MP_FREAD_C
 +--->BN_MP_ZERO_C
 +--->BN_MP_MUL_D_C
diff --git a/demo/test.c b/demo/test.c
index 3c61ff512..609030b7c 100644
--- a/demo/test.c
+++ b/demo/test.c
@@ -1643,6 +1643,172 @@ static int test_mp_balance_mul(void)
    return EXIT_FAILURE;
 }
 
+static int test_mp_factors_sort(void)
+{
+   int e, i, j, k;
+   mp_factors factors;
+   mp_int t, primorial;
+   const char *p_29 = "6469693230";
+   const mp_digit tests[9][10] = {
+      {2, 3, 5, 7, 11, 13, 17, 19, 23, 29},
+      {29, 23, 19, 17, 13, 11, 7, 5, 3, 2},
+      {19, 2, 29, 3, 23, 17, 5, 11, 13, 7},
+      {23, 2, 11, 5, 7, 19, 17, 3, 13, 29},
+      {3, 2, 5, 7, 11, 13, 17, 19, 23, 29},
+      {2, 3, 5, 7, 11, 13, 17, 19, 29, 23},
+      {2, 2, 2, 2,  2,  2,  2,  2,  2, 23},
+      {23, 2, 2, 2,  2,  2,  2,  2,  2, 2},
+      {2, 2, 2, 2,  2,  23,  2,  2,  2, 2},
+   };
+   unsigned long p_2_23 = 11776uL;
+
+   if ((e = mp_factors_init(&factors)) != MP_OKAY) {
+      return e;
+   }
+   if ((e = mp_init_multi(&t, &primorial, NULL)) != MP_OKAY) {
+      mp_factors_clear(&factors);
+      return e;
+   }
+   if ((e = mp_read_radix(&primorial, p_29, 10)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+   for (i = 0; i < 6; i++) {
+      mp_factors_zero(&factors);
+      for (j = 0; j < 10; j++) {
+         mp_set(&t,tests[i][j]);
+         if ((e = mp_factors_add(&t, &factors)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+      }
+      mp_factors_sort(&factors);
+      if ((e = mp_factors_product(&factors, &t)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (mp_cmp(&t, &primorial) != MP_EQ) {
+         goto LTM_ERR;
+      }
+      for (k = 0; k < (factors.length - 1); k++) {
+         if (mp_cmp(&(factors.factors[k]), &(factors.factors[k+1])) == MP_GT) {
+            goto LTM_ERR;
+         }
+      }
+   }
+
+   if ((e = mp_set_long(&primorial, p_2_23)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+   for (i = 6; i < 9; i++) {
+      mp_factors_zero(&factors);
+      for (j = 0; j < 10; j++) {
+         mp_set(&t,tests[i][j]);
+         if ((e = mp_factors_add(&t, &factors)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+      }
+      mp_factors_sort(&factors);
+      if ((e = mp_factors_product(&factors, &t)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (mp_cmp(&t, &primorial) != MP_EQ) {
+         goto LTM_ERR;
+      }
+      for (k = 0; k < (factors.length - 1); k++) {
+         if (mp_cmp(&(factors.factors[k]), &(factors.factors[k+1])) == MP_GT) {
+            goto LTM_ERR;
+         }
+      }
+   }
+
+   mp_clear_multi(&t, &primorial, NULL);
+   mp_factors_clear(&factors);
+   return EXIT_SUCCESS;
+LTM_ERR:
+   mp_clear_multi(&t, &primorial, NULL);
+   mp_factors_clear(&factors);
+   return EXIT_FAILURE;
+}
+
+static int test_mp_factors_product(void)
+{
+   int e;
+   int i, j;
+   mp_factors factors;
+   mp_int t, primorial;
+   const char *p_29 = "6469693230";
+#ifndef MP_8BIT
+   const char *p_1000 =
+      "lWyCCT0aF3uOiLv+kFrAb/EkIbfXTJx+MMPOxI1Dl8aLBCvFRKN/CUTUHwyDdk1sdC8rdjKRFGTd3wE46/iQibo1eoA6l7QW3H2PpQPNb46iK1m2/2Rmhr+Cmxu8eB+KHZoXa7H8pEXdIfufN0hWvfw7o2Sph35IoxakyRoz0nCZiPXLlt4tVsMNHljt+hK67sycluoiFhWR+iCJhtFqX1I21PJok/eNSkr2QM";
+#endif
+   const mp_digit tests_a[3][10] = {
+      {2, 3, 5, 7, 11, 13, 17, 19, 23, 29},
+      {29, 23, 19, 17, 13, 11, 7, 5, 3, 2},
+      {29, 23, 19, 17,  0, 11, 7, 5, 3, 2}
+   };
+
+   if ((e = mp_factors_init(&factors)) != MP_OKAY) {
+      return e;
+   }
+   if ((e = mp_init_multi(&t, &primorial, NULL)) != MP_OKAY) {
+      mp_factors_clear(&factors);
+      return e;
+   }
+
+   /* zero length */
+   if ((e = mp_factors_product(&factors, &t)) == MP_OKAY) {
+      goto LTM_ERR;
+   }
+   /* zero in array */
+   mp_factors_zero(&factors);
+   for (j = 0; j < 10; j++) {
+      mp_set(&t,tests_a[2][j]);
+      mp_factors_add(&t, &factors);
+   }
+   mp_factors_product(&factors, &t);
+   if (mp_cmp_d(&t, 0uL) != MP_EQ) {
+      goto LTM_ERR;
+   }
+
+   mp_read_radix(&primorial, p_29, 10);
+
+   for (i = 0; i < 2; i++) {
+      mp_factors_zero(&factors);
+      for (j = 0; j < 10; j++) {
+         mp_set(&t,tests_a[i][j]);
+         mp_factors_add(&t, &factors);
+      }
+      mp_factors_product(&factors, &t);
+      if (mp_cmp(&t, &primorial) != MP_EQ) {
+         puts("AAA");
+         goto LTM_ERR;
+      }
+   }
+
+   /* not enough primes < 127 to trigger the binary-splitting */
+#ifndef MP_8BIT
+   mp_factors_zero(&factors);
+   for (i = 0; i < 168; i++) {
+      mp_set(&t, ltm_prime_tab[i]);
+      mp_factors_add(&t, &factors);
+   }
+   mp_factors_product(&factors, &t);
+   mp_read_radix(&primorial, p_1000, 64);
+   if (mp_cmp(&t, &primorial) != MP_EQ) {
+      goto LTM_ERR;
+   }
+#endif
+
+   mp_clear_multi(&t, &primorial, NULL);
+   mp_factors_clear(&factors);
+   return EXIT_SUCCESS;
+LTM_ERR:
+   mp_clear_multi(&t, &primorial, NULL);
+   mp_factors_clear(&factors);
+   return EXIT_FAILURE;
+}
+
+
 int unit_tests(void)
 {
    static const struct {
@@ -1678,7 +1844,9 @@ int unit_tests(void)
       T(mp_tc_xor),
       T(mp_incr),
       T(mp_decr),
-      T(mp_balance_mul)
+      T(mp_balance_mul),
+      T(mp_factors_sort),
+      T(mp_factors_product)
 #undef T
    };
    unsigned long i;
diff --git a/doc/bn.tex b/doc/bn.tex
index c83b2e399..8b86ab2b0 100644
--- a/doc/bn.tex
+++ b/doc/bn.tex
@@ -2116,6 +2116,57 @@ \subsection{To ASCII}
 int mp_fwrite(const mp_int *a, int radix, FILE *stream);
 \end{alltt}
 
+\chapter{Factorizing}
+\subsection{Managing Lists of Factors}
+The data structure used for the list of factors is a simple array, packed into a \texttt{struct} to surround it with some useful information like the number of elements (\texttt{length}) and the amount of memory allocated (\texttt{alloc}), also counted as a number of elements.
+\begin{alltt}
+typedef struct  {
+   int length, alloc;
+   mp_int *factors;
+} mp_factors;
+\end{alltt}
+The macro \texttt{LTM\_TRIAL\_GROWTH} holds the amount of memory, measured in numbers of elements, by which \texttt{realloc} will increase the memory available to the list \texttt{factors}. It is defined in \texttt{tommath.h} and the default value is $64$.
+
+A couple of functions are available to make the use of it a bit more comfortable.
+
+\begin{description}
+\item
+\index{mp\_factors\_init}
+\verb!int mp_factors_init(mp_factors *f);!\\
+Initialize the factor list by allocating a certain amount of memory and setting \texttt{length = 0} and \texttt{alloc} to the amount of memory pre-allocated. The exact amount is defined at compile time by the macro \texttt{LTM\_TRIAL\_GROWTH} in \texttt{tommath.h}.
+\item
+\index{mp\_factors\_clear}
+\verb!void mp_factors_clear(mp_factors *f);!\\
+This function free's all memory used.
+\item
+\index{mp\_factors\_zero}
+\verb!int mp_factors_zero(mp_factors *f);!\\
+Remove the elements of the factor list and allocate (fresh) memory of default size in that order.
+\item
+\index{mp\_factors\_add}
+\verb!int mp_factors_add(const mp_int *a, mp_factors *f);!\\
+Add a factor of type \texttt{mp\_int} to the list.
+\item
+\index{mp\_factors\_sort}
+\verb!int mp_factors_sort(mp_factors *f);!\\
+The factors in the list are not necessarily in increasing order. This functions changes that. It does it with the insert-sort algorithm, a good choice for the task it has been written for (the list is most likely already ordered) but not for many other tasks involving large lists in random order.
+\item
+\index{mp\_factors\_print}
+\texttt{int mp\_factors\_print(mp\_factors *f, int base, char delimiter,\\
+\hphantom{int mp\_factors\_print(} FILE *stream);}\\
+Prints the element of the list in base \texttt{base} to \texttt{stream} with the delimiter \texttt{delimiter}. The default delimiter is a comma (ASCII \texttt{0x2c}).
+\item
+\index{mp\_factors\_product}
+\verb!int mp_factors_product(const mp_factors *factors, mp_int *p);!\\
+Multiplies all elements of the list. It does not recognize sparse lists, every zero in the list gets multiplied, too. It does multiply the list with a binary-splitting algorithm which assumes a highly or better fully sorted list to work optimally.
+\item
+\index{mp\_factors\_compress}
+\verb!int mp_factors_compress(const mp_factors *f, int sorted, mp_factors *g)!\\
+Change (sorted) factors array into format ``prime, prime-count, \textellipsis, prime, prime-count''. The type of prime-count is therefore a bigint, too, which makes it useful for low \texttt{MP\_xBIT}.
+\end{description}
+
+
+
 
 \subsection{From ASCII}
 \index{mp\_read\_radix}
diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj
index 69bce8769..6a1da0897 100644
--- a/libtommath_VS2008.vcproj
+++ b/libtommath_VS2008.vcproj
@@ -468,6 +468,38 @@
 			RelativePath="bn_mp_exteuclid.c"
 			>
 		</File>
+		<File
+			RelativePath="bn_mp_factors_add.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_clear.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_compress.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_init.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_print.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_product.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_sort.c"
+			>
+		</File>
+		<File
+			RelativePath="bn_mp_factors_zero.c"
+			>
+		</File>
 		<File
 			RelativePath="bn_mp_fread.c"
 			>
diff --git a/makefile b/makefile
index 4ca43ecad..fdf1aadc9 100644
--- a/makefile
+++ b/makefile
@@ -32,14 +32,16 @@ bn_mp_addmod.o bn_mp_and.o bn_mp_balance_mul.o bn_mp_clamp.o bn_mp_clear.o bn_mp
 bn_mp_cmp_d.o bn_mp_cmp_mag.o bn_mp_cnt_lsb.o bn_mp_complement.o bn_mp_copy.o bn_mp_count_bits.o \
 bn_mp_decr.o bn_mp_div.o bn_mp_div_2.o bn_mp_div_2d.o bn_mp_div_3.o bn_mp_div_d.o bn_mp_dr_is_modulus.o \
 bn_mp_dr_reduce.o bn_mp_dr_setup.o bn_mp_exch.o bn_mp_export.o bn_mp_expt_d.o bn_mp_expt_d_ex.o \
-bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o \
-bn_mp_get_bit.o bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o \
-bn_mp_import.o bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o \
-bn_mp_iseven.o bn_mp_isodd.o bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o \
-bn_mp_lcm.o bn_mp_lshd.o bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o \
-bn_mp_montgomery_reduce.o bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o \
-bn_mp_mulmod.o bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
+bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_factors_add.o bn_mp_factors_clear.o \
+bn_mp_factors_compress.o bn_mp_factors_init.o bn_mp_factors_print.o bn_mp_factors_product.o \
+bn_mp_factors_sort.o bn_mp_factors_zero.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o bn_mp_get_bit.o \
+bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o bn_mp_import.o \
+bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o \
+bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o \
+bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \
+bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \
+bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \
+bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
 bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \
 bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \
 bn_mp_prime_random_ex.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \
diff --git a/makefile.mingw b/makefile.mingw
index bb2552d95..cf23c6313 100644
--- a/makefile.mingw
+++ b/makefile.mingw
@@ -35,14 +35,16 @@ bn_mp_addmod.o bn_mp_and.o bn_mp_balance_mul.o bn_mp_clamp.o bn_mp_clear.o bn_mp
 bn_mp_cmp_d.o bn_mp_cmp_mag.o bn_mp_cnt_lsb.o bn_mp_complement.o bn_mp_copy.o bn_mp_count_bits.o \
 bn_mp_decr.o bn_mp_div.o bn_mp_div_2.o bn_mp_div_2d.o bn_mp_div_3.o bn_mp_div_d.o bn_mp_dr_is_modulus.o \
 bn_mp_dr_reduce.o bn_mp_dr_setup.o bn_mp_exch.o bn_mp_export.o bn_mp_expt_d.o bn_mp_expt_d_ex.o \
-bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o \
-bn_mp_get_bit.o bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o \
-bn_mp_import.o bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o \
-bn_mp_iseven.o bn_mp_isodd.o bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o \
-bn_mp_lcm.o bn_mp_lshd.o bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o \
-bn_mp_montgomery_reduce.o bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o \
-bn_mp_mulmod.o bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
+bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_factors_add.o bn_mp_factors_clear.o \
+bn_mp_factors_compress.o bn_mp_factors_init.o bn_mp_factors_print.o bn_mp_factors_product.o \
+bn_mp_factors_sort.o bn_mp_factors_zero.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o bn_mp_get_bit.o \
+bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o bn_mp_import.o \
+bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o \
+bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o \
+bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \
+bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \
+bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \
+bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
 bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \
 bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \
 bn_mp_prime_random_ex.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \
diff --git a/makefile.msvc b/makefile.msvc
index 43f1c6d44..a7bb967d0 100644
--- a/makefile.msvc
+++ b/makefile.msvc
@@ -27,14 +27,16 @@ bn_mp_addmod.obj bn_mp_and.obj bn_mp_balance_mul.obj bn_mp_clamp.obj bn_mp_clear
 bn_mp_cmp_d.obj bn_mp_cmp_mag.obj bn_mp_cnt_lsb.obj bn_mp_complement.obj bn_mp_copy.obj bn_mp_count_bits.obj \
 bn_mp_decr.obj bn_mp_div.obj bn_mp_div_2.obj bn_mp_div_2d.obj bn_mp_div_3.obj bn_mp_div_d.obj bn_mp_dr_is_modulus.obj \
 bn_mp_dr_reduce.obj bn_mp_dr_setup.obj bn_mp_exch.obj bn_mp_export.obj bn_mp_expt_d.obj bn_mp_expt_d_ex.obj \
-bn_mp_exptmod.obj bn_mp_exptmod_fast.obj bn_mp_exteuclid.obj bn_mp_fread.obj bn_mp_fwrite.obj bn_mp_gcd.obj \
-bn_mp_get_bit.obj bn_mp_get_double.obj bn_mp_get_int.obj bn_mp_get_long.obj bn_mp_get_long_long.obj bn_mp_grow.obj \
-bn_mp_import.obj bn_mp_incr.obj bn_mp_init.obj bn_mp_init_copy.obj bn_mp_init_multi.obj bn_mp_init_set.obj \
-bn_mp_init_set_int.obj bn_mp_init_size.obj bn_mp_invmod.obj bn_mp_invmod_slow.obj bn_mp_is_square.obj \
-bn_mp_iseven.obj bn_mp_isodd.obj bn_mp_jacobi.obj bn_mp_karatsuba_mul.obj bn_mp_karatsuba_sqr.obj bn_mp_kronecker.obj \
-bn_mp_lcm.obj bn_mp_lshd.obj bn_mp_mod.obj bn_mp_mod_2d.obj bn_mp_mod_d.obj bn_mp_montgomery_calc_normalization.obj \
-bn_mp_montgomery_reduce.obj bn_mp_montgomery_setup.obj bn_mp_mul.obj bn_mp_mul_2.obj bn_mp_mul_2d.obj bn_mp_mul_d.obj \
-bn_mp_mulmod.obj bn_mp_n_root.obj bn_mp_n_root_ex.obj bn_mp_neg.obj bn_mp_or.obj bn_mp_prime_fermat.obj \
+bn_mp_exptmod.obj bn_mp_exptmod_fast.obj bn_mp_exteuclid.obj bn_mp_factors_add.obj bn_mp_factors_clear.obj \
+bn_mp_factors_compress.obj bn_mp_factors_init.obj bn_mp_factors_print.obj bn_mp_factors_product.obj \
+bn_mp_factors_sort.obj bn_mp_factors_zero.obj bn_mp_fread.obj bn_mp_fwrite.obj bn_mp_gcd.obj bn_mp_get_bit.obj \
+bn_mp_get_double.obj bn_mp_get_int.obj bn_mp_get_long.obj bn_mp_get_long_long.obj bn_mp_grow.obj bn_mp_import.obj \
+bn_mp_incr.obj bn_mp_init.obj bn_mp_init_copy.obj bn_mp_init_multi.obj bn_mp_init_set.obj bn_mp_init_set_int.obj \
+bn_mp_init_size.obj bn_mp_invmod.obj bn_mp_invmod_slow.obj bn_mp_is_square.obj bn_mp_iseven.obj bn_mp_isodd.obj \
+bn_mp_jacobi.obj bn_mp_karatsuba_mul.obj bn_mp_karatsuba_sqr.obj bn_mp_kronecker.obj bn_mp_lcm.obj bn_mp_lshd.obj \
+bn_mp_mod.obj bn_mp_mod_2d.obj bn_mp_mod_d.obj bn_mp_montgomery_calc_normalization.obj bn_mp_montgomery_reduce.obj \
+bn_mp_montgomery_setup.obj bn_mp_mul.obj bn_mp_mul_2.obj bn_mp_mul_2d.obj bn_mp_mul_d.obj bn_mp_mulmod.obj \
+bn_mp_n_root.obj bn_mp_n_root_ex.obj bn_mp_neg.obj bn_mp_or.obj bn_mp_prime_fermat.obj \
 bn_mp_prime_frobenius_underwood.obj bn_mp_prime_is_divisible.obj bn_mp_prime_is_prime.obj \
 bn_mp_prime_miller_rabin.obj bn_mp_prime_next_prime.obj bn_mp_prime_rabin_miller_trials.obj \
 bn_mp_prime_random_ex.obj bn_mp_prime_strong_lucas_selfridge.obj bn_mp_radix_size.obj bn_mp_radix_smap.obj \
diff --git a/makefile.shared b/makefile.shared
index 012cd5bff..533aefa61 100644
--- a/makefile.shared
+++ b/makefile.shared
@@ -29,14 +29,16 @@ bn_mp_addmod.o bn_mp_and.o bn_mp_balance_mul.o bn_mp_clamp.o bn_mp_clear.o bn_mp
 bn_mp_cmp_d.o bn_mp_cmp_mag.o bn_mp_cnt_lsb.o bn_mp_complement.o bn_mp_copy.o bn_mp_count_bits.o \
 bn_mp_decr.o bn_mp_div.o bn_mp_div_2.o bn_mp_div_2d.o bn_mp_div_3.o bn_mp_div_d.o bn_mp_dr_is_modulus.o \
 bn_mp_dr_reduce.o bn_mp_dr_setup.o bn_mp_exch.o bn_mp_export.o bn_mp_expt_d.o bn_mp_expt_d_ex.o \
-bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o \
-bn_mp_get_bit.o bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o \
-bn_mp_import.o bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o \
-bn_mp_iseven.o bn_mp_isodd.o bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o \
-bn_mp_lcm.o bn_mp_lshd.o bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o \
-bn_mp_montgomery_reduce.o bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o \
-bn_mp_mulmod.o bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
+bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_factors_add.o bn_mp_factors_clear.o \
+bn_mp_factors_compress.o bn_mp_factors_init.o bn_mp_factors_print.o bn_mp_factors_product.o \
+bn_mp_factors_sort.o bn_mp_factors_zero.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o bn_mp_get_bit.o \
+bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o bn_mp_import.o \
+bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o \
+bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o \
+bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \
+bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \
+bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \
+bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
 bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \
 bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \
 bn_mp_prime_random_ex.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \
diff --git a/makefile.unix b/makefile.unix
index 856e18cee..bb428eec2 100644
--- a/makefile.unix
+++ b/makefile.unix
@@ -36,14 +36,16 @@ bn_mp_addmod.o bn_mp_and.o bn_mp_balance_mul.o bn_mp_clamp.o bn_mp_clear.o bn_mp
 bn_mp_cmp_d.o bn_mp_cmp_mag.o bn_mp_cnt_lsb.o bn_mp_complement.o bn_mp_copy.o bn_mp_count_bits.o \
 bn_mp_decr.o bn_mp_div.o bn_mp_div_2.o bn_mp_div_2d.o bn_mp_div_3.o bn_mp_div_d.o bn_mp_dr_is_modulus.o \
 bn_mp_dr_reduce.o bn_mp_dr_setup.o bn_mp_exch.o bn_mp_export.o bn_mp_expt_d.o bn_mp_expt_d_ex.o \
-bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o \
-bn_mp_get_bit.o bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o \
-bn_mp_import.o bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o \
-bn_mp_init_set_int.o bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o \
-bn_mp_iseven.o bn_mp_isodd.o bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o \
-bn_mp_lcm.o bn_mp_lshd.o bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o \
-bn_mp_montgomery_reduce.o bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o \
-bn_mp_mulmod.o bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
+bn_mp_exptmod.o bn_mp_exptmod_fast.o bn_mp_exteuclid.o bn_mp_factors_add.o bn_mp_factors_clear.o \
+bn_mp_factors_compress.o bn_mp_factors_init.o bn_mp_factors_print.o bn_mp_factors_product.o \
+bn_mp_factors_sort.o bn_mp_factors_zero.o bn_mp_fread.o bn_mp_fwrite.o bn_mp_gcd.o bn_mp_get_bit.o \
+bn_mp_get_double.o bn_mp_get_int.o bn_mp_get_long.o bn_mp_get_long_long.o bn_mp_grow.o bn_mp_import.o \
+bn_mp_incr.o bn_mp_init.o bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o \
+bn_mp_init_size.o bn_mp_invmod.o bn_mp_invmod_slow.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o \
+bn_mp_jacobi.o bn_mp_karatsuba_mul.o bn_mp_karatsuba_sqr.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \
+bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \
+bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \
+bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \
 bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \
 bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \
 bn_mp_prime_random_ex.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \
diff --git a/tommath.h b/tommath.h
index 41c9ec0b2..b4e01856d 100644
--- a/tommath.h
+++ b/tommath.h
@@ -567,6 +567,48 @@ int mp_prime_next_prime(mp_int *a, int t, int bbs_style);
  */
 int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat);
 
+/*
+   Simple data structure to hold some mp_int's.
+   Is used in LTM to hold the factors from factoring, about
+   100 or so max (large primorials of powers notwithstanding),
+   always with linear r/w access.
+ */
+typedef struct  {
+   int length, alloc;
+   mp_int *factors;
+} mp_factors;
+
+#ifndef LTM_TRIAL_GROWTH
+#define LTM_TRIAL_GROWTH 64
+#endif
+
+/* Init factor-array (an array of mp_int's) */
+int mp_factors_init(mp_factors *f);
+
+/* Clear and free factor-array (an array of mp_int's) */
+void mp_factors_clear(mp_factors *f);
+
+/* Zero the elements of a factor-array (an array of mp_int's) and realloc to def
+   size in that order */
+int mp_factors_zero(mp_factors *f);
+
+/* Add one element (mp_int) to a factor-array (an array of mp_int's) */
+int mp_factors_add(const mp_int *a, mp_factors *f);
+
+/* Sort factor-array (an array of mp_int's) */
+void mp_factors_sort(mp_factors *f);
+
+/* Factor integer. Good for up to about 35-40 bit large factors */
+int mp_factor(const mp_int *z, mp_factors *factors);
+
+/* Compute the product of the elements of "factors" */
+int mp_factors_product(const mp_factors *factors, mp_int *p);
+
+/* Change (sorted) factors array into format prime,prime-count,...,prime,prime-count */
+int mp_factors_compress(mp_factors *f, int sorted, mp_factors *g);
+
+
+
 /* ---> radix conversion <--- */
 int mp_count_bits(const mp_int *a);
 
@@ -588,6 +630,8 @@ int mp_radix_size(const mp_int *a, int radix, int *size);
 #ifndef LTM_NO_FILE
 int mp_fread(mp_int *a, int radix, FILE *stream);
 int mp_fwrite(const mp_int *a, int radix, FILE *stream);
+/* Print the elememts of a factor-array (an array of mp_int's) */
+int mp_factors_print(const mp_factors *f, int base, char delimiter, FILE *stream);
 #endif
 
 #define mp_read_raw(mp, str, len) mp_read_signed_bin((mp), (str), (len))
diff --git a/tommath_class.h b/tommath_class.h
index 56a6cf239..df8fec78b 100644
--- a/tommath_class.h
+++ b/tommath_class.h
@@ -49,6 +49,14 @@
 #   define BN_MP_EXPTMOD_C
 #   define BN_MP_EXPTMOD_FAST_C
 #   define BN_MP_EXTEUCLID_C
+#   define BN_MP_FACTORS_ADD_C
+#   define BN_MP_FACTORS_CLEAR_C
+#   define BN_MP_FACTORS_COMPRESS_C
+#   define BN_MP_FACTORS_INIT_C
+#   define BN_MP_FACTORS_PRINT_C
+#   define BN_MP_FACTORS_PRODUCT_C
+#   define BN_MP_FACTORS_SORT_C
+#   define BN_MP_FACTORS_ZERO_C
 #   define BN_MP_FREAD_C
 #   define BN_MP_FWRITE_C
 #   define BN_MP_GCD_C
@@ -421,6 +429,58 @@
 #   define BN_MP_CLEAR_MULTI_C
 #endif
 
+#if defined(BN_MP_FACTORS_ADD_C)
+#   define BN_MP_FACTORS_C
+#   define BN_MP_INIT_C
+#   define BN_MP_COPY_C
+#endif
+
+#if defined(BN_MP_FACTORS_CLEAR_C)
+#   define BN_MP_CLEAR_C
+#endif
+
+#if defined(BN_MP_FACTORS_COMPRESS_C)
+#   define BN_MP_INIT_C
+#   define BN_MP_FACTORS_INIT_C
+#   define BN_MP_FACTORS_ADD_C
+#   define BN_MP_FACTORS_SORT_C
+#   define BN_MP_INCR_C
+#   define BN_MP_CMP_C
+#   define BN_MP_ZERO_C
+#   define BN_MP_FACTORS_CLEAR_C
+#   define BN_MP_CLEAR_C
+#endif
+
+#if defined(BN_MP_FACTORS_INIT_C)
+#   define BN_MP_FACTORS_C
+#endif
+
+#if defined(BN_MP_FACTORS_PRINT_C)
+#   define BN_MP_FACTORS_C
+#   define BN_MP_FWRITE_C
+#endif
+
+#if defined(BN_MP_FACTORS_PRODUCT_C)
+#   define BN_MP_FACTORS_C
+#   define BN_MP_SET_C
+#   define BN_MP_SET_LONG_C
+#   define BN_MP_MUL_C
+#   define BN_MP_INIT_C
+#   define BN_MP_CLEAR_C
+#   define BN_MP_COPY_C
+#endif
+
+#if defined(BN_MP_FACTORS_SORT_C)
+#   define BN_MP_FACTORS_C
+#   define BN_MP_CMP_C
+#   define BN_MP_EXCH_C
+#endif
+
+#if defined(BN_MP_FACTORS_ZERO_C)
+#   define BN_MP_FACTORS_C
+#   define BN_MP_CLEAR_C
+#endif
+
 #if defined(BN_MP_FREAD_C)
 #   define BN_MP_ZERO_C
 #   define BN_MP_S_RMAP_REVERSE_SZ_C