diff --git a/demo/test.c b/demo/test.c
index 2fa6e08d..6ec6fdd3 100644
--- a/demo/test.c
+++ b/demo/test.c
@@ -1930,6 +1930,26 @@ static int test_mp_root_n(void)
          EXPECT(mp_cmp(&r, &c) == MP_EQ);
       }
    }
+
+   /* 0^(1/x) = 0 with x != 0 is allowed, test */
+   mp_set(&a, 0);
+   DO(mp_root_n(&a, 2, &c));
+   EXPECT(mp_cmp_d(&c, 0) == MP_EQ);
+
+   /* Not allowed: division by zero */
+   mp_set(&a, 2);
+   EXPECT(mp_root_n(&a, 0, &c) == MP_VAL);
+
+   /* root^base == input with small input and base */
+   mp_set(&a, 4);
+   DO(mp_root_n(&a, 2, &c));
+   EXPECT(mp_cmp_d(&c, 2) == MP_EQ);
+
+   /* (root^base)^(1/(base + 1)) with small root */
+   DO(mp_2expt(&a, 48));
+   DO(mp_root_n(&a, 49, &c));
+   EXPECT(mp_cmp_d(&c, 1) == MP_EQ);
+
    mp_clear_multi(&a, &c, &r, NULL);
    return EXIT_SUCCESS;
 LBL_ERR:
diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj
index 080b7e0d..97a99d60 100644
--- a/libtommath_VS2008.vcproj
+++ b/libtommath_VS2008.vcproj
@@ -836,6 +836,10 @@
 			RelativePath="s_mp_exptmod_fast.c"
 			>
 		</File>
+		<File
+			RelativePath="s_mp_fp_exp2.c"
+			>
+		</File>
 		<File
 			RelativePath="s_mp_fp_log.c"
 			>
@@ -916,6 +920,10 @@
 			RelativePath="s_mp_rand_source.c"
 			>
 		</File>
+		<File
+			RelativePath="s_mp_root_n.c"
+			>
+		</File>
 		<File
 			RelativePath="s_mp_sqr.c"
 			>
diff --git a/makefile b/makefile
index 6daa522b..8e4619cc 100644
--- a/makefile
+++ b/makefile
@@ -44,13 +44,14 @@ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o
 mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \
 mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \
 mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
-s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \
-s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_exp2.o \
+s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
 s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \
 s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \
-s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \
-s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \
-s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_root_n.o \
+s_mp_sqr.o s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o \
+s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+
 
 #END_INS
 
diff --git a/makefile.mingw b/makefile.mingw
index 10f2bd59..b78e12b6 100644
--- a/makefile.mingw
+++ b/makefile.mingw
@@ -46,13 +46,14 @@ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o
 mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \
 mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \
 mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
-s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \
-s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_exp2.o \
+s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
 s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \
 s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \
-s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \
-s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \
-s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_root_n.o \
+s_mp_sqr.o s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o \
+s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+
 
 HEADERS_PUB=tommath.h
 HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)
diff --git a/makefile.msvc b/makefile.msvc
index f160f2cc..4398633d 100644
--- a/makefile.msvc
+++ b/makefile.msvc
@@ -42,13 +42,13 @@ mp_reduce_setup.obj mp_root_n.obj mp_rshd.obj mp_sbin_size.obj mp_set.obj mp_set
 mp_set_l.obj mp_set_u32.obj mp_set_u64.obj mp_set_ul.obj mp_shrink.obj mp_signed_rsh.obj mp_sqrmod.obj mp_sqrt.obj \
 mp_sqrtmod_prime.obj mp_sub.obj mp_sub_d.obj mp_submod.obj mp_to_radix.obj mp_to_sbin.obj mp_to_ubin.obj mp_ubin_size.obj \
 mp_unpack.obj mp_warray_free.obj mp_xor.obj mp_zero.obj s_mp_add.obj s_mp_copy_digs.obj s_mp_div_3.obj \
-s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj s_mp_exptmod.obj s_mp_exptmod_fast.obj s_mp_fp_log.obj \
-s_mp_fp_log_d.obj s_mp_get_bit.obj s_mp_invmod.obj s_mp_invmod_odd.obj s_mp_log_2expt.obj \
+s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj s_mp_exptmod.obj s_mp_exptmod_fast.obj s_mp_fp_exp2.obj \
+s_mp_fp_log.obj s_mp_fp_log_d.obj s_mp_get_bit.obj s_mp_invmod.obj s_mp_invmod_odd.obj s_mp_log_2expt.obj \
 s_mp_montgomery_reduce_comba.obj s_mp_mul.obj s_mp_mul_balance.obj s_mp_mul_comba.obj s_mp_mul_high.obj \
 s_mp_mul_high_comba.obj s_mp_mul_karatsuba.obj s_mp_mul_toom.obj s_mp_prime_is_divisible.obj s_mp_prime_tab.obj \
-s_mp_radix_map.obj s_mp_radix_size_overestimate.obj s_mp_rand_platform.obj s_mp_rand_source.obj s_mp_sqr.obj \
-s_mp_sqr_comba.obj s_mp_sqr_karatsuba.obj s_mp_sqr_toom.obj s_mp_sub.obj s_mp_warray.obj s_mp_warray_get.obj \
-s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj
+s_mp_radix_map.obj s_mp_radix_size_overestimate.obj s_mp_rand_platform.obj s_mp_rand_source.obj s_mp_root_n.obj \
+s_mp_sqr.obj s_mp_sqr_comba.obj s_mp_sqr_karatsuba.obj s_mp_sqr_toom.obj s_mp_sub.obj s_mp_warray.obj \
+s_mp_warray_get.obj s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj
 
 HEADERS_PUB=tommath.h
 HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)
diff --git a/makefile.shared b/makefile.shared
index 20bd8e03..0f7909dc 100644
--- a/makefile.shared
+++ b/makefile.shared
@@ -41,13 +41,14 @@ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o
 mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \
 mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \
 mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
-s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \
-s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_exp2.o \
+s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
 s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \
 s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \
-s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \
-s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \
-s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_root_n.o \
+s_mp_sqr.o s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o \
+s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+
 
 #END_INS
 
diff --git a/makefile.unix b/makefile.unix
index 06dadddb..aade2d9d 100644
--- a/makefile.unix
+++ b/makefile.unix
@@ -47,13 +47,14 @@ mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o
 mp_set_l.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_shrink.o mp_signed_rsh.o mp_sqrmod.o mp_sqrt.o \
 mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o mp_to_ubin.o mp_ubin_size.o \
 mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
-s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \
-s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_exp2.o \
+s_mp_fp_log.o s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \
 s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \
 s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \
-s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \
-s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \
-s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_root_n.o \
+s_mp_sqr.o s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o \
+s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
+
 
 
 HEADERS_PUB=tommath.h
diff --git a/mp_root_n.c b/mp_root_n.c
index d904df88..7440709b 100644
--- a/mp_root_n.c
+++ b/mp_root_n.c
@@ -6,136 +6,11 @@
 /* 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_n(const mp_int *a, int b, mp_int *c)
 {
-   mp_int t1, t2, t3, a_;
-   int    ilog2;
-   mp_err err;
-
-   if (b < 0 || (unsigned)b > (unsigned)MP_DIGIT_MAX) {
-      return MP_VAL;
-   }
-
-   /* input must be positive if b is even */
-   if (((b & 1) == 0) && mp_isneg(a)) {
-      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 > 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 < b) {
-      mp_set(c, 1uL);
-      c->sign = a->sign;
-      err = MP_OKAY;
-      goto LBL_ERR;
-   }
-   ilog2 =  ilog2 / 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_n(&t1, b - 1, &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, (mp_digit)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 (;;) {
-      mp_ord cmp;
-      if ((err = mp_expt_n(&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_n(&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;
-
-LBL_ERR:
-   mp_clear_multi(&t1, &t2, &t3, NULL);
-   return err;
+   return s_mp_root_n(a, b, c, NULL);
 }
 
 #endif
diff --git a/s_mp_fp_exp2.c b/s_mp_fp_exp2.c
new file mode 100644
index 00000000..3480081c
--- /dev/null
+++ b/s_mp_fp_exp2.c
@@ -0,0 +1,116 @@
+#include "tommath_private.h"
+#ifdef S_MP_FP_EXP2_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+static MP_SET_UNSIGNED(mp_set_word, mp_word)
+
+static mp_err s_mp_fp_exp2_fract(const mp_int *a, mp_int *c)
+{
+   mp_err err = MP_OKAY;
+   mp_int u, tmp;
+   /*
+      Generated with Sam Hocevar's LolRemez at https://github.com/samhocevar/lolremez
+
+      lolremez --double -d 10 -r "0:1" "2**(x)"
+         and
+      lolremez --double -d 5 -r "0:1" "2**(x)"
+
+   */
+#if (MP_DIGIT_BIT == 15)
+   mp_word coefficients[] = {
+      0x3e,
+      0x125,
+      0x726,
+      0x1ebc,
+      0x58b9,
+      0x7fff
+   };
+#elif ((MP_DIGIT_BIT == 28) || (MP_DIGIT_BIT == 31))
+   mp_word coefficients[] = {
+      0x15,
+      0xca,
+      0xb2c,
+      0x7fe0,
+      0x50c2e,
+      0x2bb0fc,
+      0x13b2ab7,
+      0x71ac235,
+      0x1ebfbdff,
+      0x58b90bfb,
+      0x80000000
+   };
+#elif (MP_DIGIT_BIT == 60)
+   mp_word coefficients[] = {
+      0x157d3d6ee0,
+      0xca97869248,
+      0xb2c1b72cee4,
+      0x7fe03bc5603e,
+      0x50c2e781f2028,
+      0x2bb0fc474348da,
+      0x13b2ab7bece61d7,
+      0x71ac235a89a9728,
+      0x1ebfbdff845aca73,
+      0x58b90bfbe8dd8d73,
+      0x8000000000000acf
+   };
+#endif
+   int coeff_size = (int)(sizeof(coefficients)/sizeof(coefficients[0]));
+   int i;
+
+   if ((err = mp_init_multi(&u, &tmp, NULL)) != MP_OKAY) {
+      return err;
+   }
+   mp_set_word(&tmp, coefficients[0]);
+   for (i = 1; i < coeff_size; i++) {
+      /* u = ((u * x) >> 63) + coefficients[i] */
+      if ((err = mp_mul(&u, a, &u)) != MP_OKAY)                                      goto LBL_ERR;
+      if ((err = mp_div_2d(&u, MP_PRECISION_FIXED_LOG, &u, NULL)) != MP_OKAY)        goto LBL_ERR;
+      mp_set_word(&tmp, coefficients[i]);
+      if ((err = mp_add(&u, &tmp, &u)) != MP_OKAY)                                   goto LBL_ERR;
+   }
+
+   if (c != NULL) {
+      mp_exch(&u, c);
+   }
+
+LBL_ERR:
+   mp_clear_multi(&u, &tmp, NULL);
+   return err;
+}
+
+/* input a is the output from s_mp_fp_log(), a fixed point log base 2 */
+mp_err s_mp_fp_exp2(const mp_int *a, mp_int *c)
+{
+   mp_err err = MP_OKAY;
+   mp_int tmp, int_part, fract_part;
+
+   if ((err = mp_init_multi(&tmp, &int_part, &fract_part, NULL)) != MP_OKAY) {
+      return err;
+   }
+   /* Cut fore and aft to get the integer part */
+   if ((err = mp_div_2d(a, MP_PRECISION_FIXED_LOG, &int_part, NULL)) != MP_OKAY)          goto LBL_ERR;
+   if ((err = mp_mul_2d(&int_part, MP_PRECISION_FIXED_LOG, &int_part)) != MP_OKAY)        goto LBL_ERR;
+   /* Fractional part */
+   if ((err = mp_sub(a, &int_part, &fract_part)) != MP_OKAY)                              goto LBL_ERR;
+   /* 2^{fractional part} */
+   if ((err = s_mp_fp_exp2_fract(&fract_part, &fract_part)) != MP_OKAY)                   goto LBL_ERR;
+   /* A = 1 << normalized integer part */
+   mp_set(&tmp,1);
+   if ((err = mp_div_2d(&int_part, MP_PRECISION_FIXED_LOG, &int_part, NULL)) != MP_OKAY)  goto LBL_ERR;
+   if ((err = mp_mul_2d(&tmp, (int)mp_get_l(&int_part), &int_part)) != MP_OKAY)           goto LBL_ERR;
+   /* B = A * fractional part */
+   if ((err = mp_mul(&int_part, &fract_part, &tmp)) != MP_OKAY)                           goto LBL_ERR;
+   /* normalize B */
+   if ((err = mp_div_2d(&tmp, MP_PRECISION_FIXED_LOG, &tmp, NULL)) != MP_OKAY)            goto LBL_ERR;
+   /* C = 2^a (Roughly. Very, very roughly) */
+   if (c != NULL) {
+      mp_exch(&tmp, c);
+   }
+
+LBL_ERR:
+   mp_clear_multi(&tmp, &int_part, &fract_part, NULL);
+   return err;
+}
+
+#endif
diff --git a/s_mp_root_n.c b/s_mp_root_n.c
new file mode 100644
index 00000000..c339d9aa
--- /dev/null
+++ b/s_mp_root_n.c
@@ -0,0 +1,167 @@
+#include "tommath_private.h"
+#ifdef S_MP_ROOT_N_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 s_mp_root_n(const mp_int *a, int b, mp_int *c, bool *is_perfect_power)
+{
+   mp_int t1, t2, t3, a_;
+   int    ilog2;
+   mp_err err = MP_OKAY;
+
+   if (is_perfect_power != NULL) {
+      *is_perfect_power = false;
+   }
+
+   if (b == 0) {
+      mp_set(c, 0);
+      return MP_VAL;
+   }
+
+   /* 0^(1/x) = 0 with x != 0 is allowed */
+   if (mp_iszero(a)) {
+      mp_set(c, 0);
+      return MP_OKAY;
+   }
+
+   if (b < 0 || (unsigned)b > (unsigned)MP_DIGIT_MAX) {
+      return MP_VAL;
+   }
+
+   /* input must be positive if b is even */
+   if (((b & 1) == 0) && mp_isneg(a)) {
+      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 > 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 < b) {
+      mp_set(c, 1uL);
+      c->sign = a->sign;
+      err = MP_OKAY;
+      goto LBL_ERR;
+   }
+
+   ilog2 =  ilog2 / b;
+   if (ilog2 == 0) {
+      mp_set(c, 1uL);
+      c->sign = a->sign;
+      err = MP_OKAY;
+      goto LBL_ERR;
+   }
+   /* log2(a), fixed point */
+   if ((err = s_mp_fp_log(&a_, &t2)) != MP_OKAY)                              goto LBL_ERR;
+
+   /* This algorithm is sensitive to start-values that are too low. */
+   mp_set(&t1, 1);
+   if ((err = mp_mul_2d(&t1,MP_PRECISION_FIXED_LOG - 1, &t1)) != MP_OKAY)     goto LBL_ERR;
+   if ((err = mp_add(&t2, &t1, &t2)) != MP_OKAY)                              goto LBL_ERR;
+
+   /* k = log2(a) / b, fixed point */
+   mp_set_i64(&t1, b);
+   if ((err = mp_mul_2d(&t1,MP_PRECISION_FIXED_LOG, &t1)) != MP_OKAY)         goto LBL_ERR;
+   if ((err = mp_mul_2d(&t2,MP_PRECISION_FIXED_LOG, &t2)) != MP_OKAY)         goto LBL_ERR;
+   if ((err = mp_div(&t2, &t1, &t2,NULL)) != MP_OKAY)                         goto LBL_ERR;
+
+   /* 2^k, the start-value */
+   if ((err = s_mp_fp_exp2(&t2, &t2)) != MP_OKAY)                             goto LBL_ERR;
+   /* Add 1 (one) unconditionally for the small roots.  */
+   if ((err = mp_incr(&t2)) != 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_n(&t1, b - 1, &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, (mp_digit)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 root we found is smaller than the actual root */
+   for (;;) {
+      mp_ord cmp;
+      if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY)                          goto LBL_ERR;
+      cmp = mp_cmp(&t2, &a_);
+      if (cmp == MP_EQ) {
+         /* It is a perfect power, flag it */
+         if (is_perfect_power != NULL) {
+            *is_perfect_power = true;
+         }
+         goto LBL_SET;
+      }
+      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_n(&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;
+      }
+   }
+LBL_SET:
+   /* set the result */
+   mp_exch(&t1, c);
+
+   /* set the sign of the result */
+   c->sign = a->sign;
+
+LBL_ERR:
+   mp_clear_multi(&t1, &t2, &t3, NULL);
+   return err;
+}
+
+#endif
diff --git a/sources.cmake b/sources.cmake
index ba3012ef..50c1acc9 100644
--- a/sources.cmake
+++ b/sources.cmake
@@ -133,6 +133,7 @@ s_mp_div_school.c
 s_mp_div_small.c
 s_mp_exptmod.c
 s_mp_exptmod_fast.c
+s_mp_fp_exp2.c
 s_mp_fp_log.c
 s_mp_fp_log_d.c
 s_mp_get_bit.c
@@ -153,6 +154,7 @@ s_mp_radix_map.c
 s_mp_radix_size_overestimate.c
 s_mp_rand_platform.c
 s_mp_rand_source.c
+s_mp_root_n.c
 s_mp_sqr.c
 s_mp_sqr_comba.c
 s_mp_sqr_karatsuba.c
diff --git a/tommath_class.h b/tommath_class.h
index 40ddfd4b..4c387579 100644
--- a/tommath_class.h
+++ b/tommath_class.h
@@ -142,6 +142,7 @@
 #   define S_MP_DIV_SMALL_C
 #   define S_MP_EXPTMOD_C
 #   define S_MP_EXPTMOD_FAST_C
+#   define S_MP_FP_EXP2_C
 #   define S_MP_FP_LOG_C
 #   define S_MP_FP_LOG_D_C
 #   define S_MP_GET_BIT_C
@@ -162,6 +163,7 @@
 #   define S_MP_RADIX_SIZE_OVERESTIMATE_C
 #   define S_MP_RAND_PLATFORM_C
 #   define S_MP_RAND_SOURCE_C
+#   define S_MP_ROOT_N_C
 #   define S_MP_SQR_C
 #   define S_MP_SQR_COMBA_C
 #   define S_MP_SQR_KARATSUBA_C
@@ -813,21 +815,7 @@
 #endif
 
 #if defined(MP_ROOT_N_C)
-#   define MP_2EXPT_C
-#   define MP_ADD_D_C
-#   define MP_CLEAR_MULTI_C
-#   define MP_CMP_C
-#   define MP_COPY_C
-#   define MP_COUNT_BITS_C
-#   define MP_DIV_C
-#   define MP_EXCH_C
-#   define MP_EXPT_N_C
-#   define MP_INIT_MULTI_C
-#   define MP_MUL_C
-#   define MP_MUL_D_C
-#   define MP_SET_C
-#   define MP_SUB_C
-#   define MP_SUB_D_C
+#   define S_MP_ROOT_N_C
 #endif
 
 #if defined(MP_RSHD_C)
@@ -1078,6 +1066,22 @@
 #   define S_MP_MONTGOMERY_REDUCE_COMBA_C
 #endif
 
+#if defined(S_MP_FP_EXP2_C)
+#   define MP_ADD_C
+#   define MP_CLEAR_MULTI_C
+#   define MP_DIV_2D_C
+#   define MP_EXCH_C
+#   define MP_GET_L_C
+#   define MP_INIT_MULTI_C
+#   define MP_MUL_2D_C
+#   define MP_MUL_C
+#   define MP_SET_C
+#   define MP_SET_WORD_C
+#   define MP_SUB_C
+#   define S_MP_FP_EXP2_FRACT_C
+#   define S_MP_ZERO_DIGS_C
+#endif
+
 #if defined(S_MP_FP_LOG_C)
 #   define MP_2EXPT_C
 #   define MP_ADD_C
@@ -1253,6 +1257,28 @@
 #   define S_MP_RAND_PLATFORM_C
 #endif
 
+#if defined(S_MP_ROOT_N_C)
+#   define MP_ADD_C
+#   define MP_ADD_D_C
+#   define MP_CLEAR_MULTI_C
+#   define MP_CMP_C
+#   define MP_COPY_C
+#   define MP_COUNT_BITS_C
+#   define MP_DIV_C
+#   define MP_EXCH_C
+#   define MP_EXPT_N_C
+#   define MP_INIT_MULTI_C
+#   define MP_MUL_2D_C
+#   define MP_MUL_C
+#   define MP_MUL_D_C
+#   define MP_SET_C
+#   define MP_SET_I64_C
+#   define MP_SUB_C
+#   define MP_SUB_D_C
+#   define S_MP_FP_EXP2_C
+#   define S_MP_FP_LOG_C
+#endif
+
 #if defined(S_MP_SQR_C)
 #   define MP_CLAMP_C
 #   define MP_CLEAR_C
diff --git a/tommath_private.h b/tommath_private.h
index be620dbc..400935eb 100644
--- a/tommath_private.h
+++ b/tommath_private.h
@@ -233,6 +233,8 @@ MP_PRIVATE mp_err s_mp_radix_size_overestimate(const mp_int *a, const int radix,
 #define MP_UPPER_LIMIT_FIXED_LOG ( (int) ( (sizeof(mp_word) * CHAR_BIT) - 1))
 MP_PRIVATE mp_err s_mp_fp_log(const mp_int *a, mp_int *c) MP_WUR;
 MP_PRIVATE mp_err s_mp_fp_log_d(const mp_int *a, mp_word *c) MP_WUR;
+MP_PRIVATE mp_err s_mp_fp_exp2(const mp_int *a, mp_int *c) MP_WUR;
+MP_PRIVATE mp_err s_mp_root_n(const mp_int *a, int b, mp_int *c, bool *is_perfect_power) MP_WUR;
 
 #ifdef MP_SMALL_STACK_SIZE