From 2a5d8edbdd5621c68c5455eae8cc1bfbc66f1d39 Mon Sep 17 00:00:00 2001 From: Christoph Zurnieden Date: Fri, 24 Apr 2026 02:07:06 +0200 Subject: [PATCH 1/2] Refinement of and bugfixes in mp_root_n: - Implementation of a fixed point exp2 function to get better start-values - fixed treatment of perfect powers - fixed 0^(1/x) == 0 with x != 0 - added test for b == 0 (division by zero) --- demo/test.c | 20 +++++ libtommath_VS2008.vcproj | 8 ++ makefile | 7 +- makefile.mingw | 7 +- makefile.msvc | 6 +- makefile.shared | 7 +- makefile.unix | 7 +- mp_root_n.c | 129 +----------------------------- s_mp_fp_exp2.c | 116 +++++++++++++++++++++++++++ s_mp_root_n.c | 167 +++++++++++++++++++++++++++++++++++++++ sources.cmake | 2 + tommath_class.h | 60 ++++++++++---- tommath_private.h | 2 + 13 files changed, 381 insertions(+), 157 deletions(-) create mode 100644 s_mp_fp_exp2.c create mode 100644 s_mp_root_n.c diff --git a/demo/test.c b/demo/test.c index 2fa6e08db..6ec6fdd37 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 080b7e0de..97a99d609 100644 --- a/libtommath_VS2008.vcproj +++ b/libtommath_VS2008.vcproj @@ -836,6 +836,10 @@ RelativePath="s_mp_exptmod_fast.c" > + + @@ -916,6 +920,10 @@ RelativePath="s_mp_rand_source.c" > + + diff --git a/makefile b/makefile index 6daa522be..c4c28e491 100644 --- a/makefile +++ b/makefile @@ -44,14 +44,15 @@ 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_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 $(LIBNAME): $(OBJECTS) diff --git a/makefile.mingw b/makefile.mingw index 10f2bd59d..dc4f53f2a 100644 --- a/makefile.mingw +++ b/makefile.mingw @@ -46,14 +46,15 @@ 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_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 f160f2ccd..944ddbb3c 100644 --- a/makefile.msvc +++ b/makefile.msvc @@ -42,11 +42,11 @@ 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_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 diff --git a/makefile.shared b/makefile.shared index 20bd8e03a..91ca53025 100644 --- a/makefile.shared +++ b/makefile.shared @@ -41,14 +41,15 @@ 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_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 objs: $(OBJECTS) diff --git a/makefile.unix b/makefile.unix index 06dadddb4..a0b350eb6 100644 --- a/makefile.unix +++ b/makefile.unix @@ -47,15 +47,16 @@ 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_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/mp_root_n.c b/mp_root_n.c index d904df883..7440709be 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 000000000..3480081c5 --- /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 000000000..c339d9aa6 --- /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 ba3012efb..50c1acc9a 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 40ddfd4b4..4531f57af 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 @@ -1249,10 +1253,36 @@ #if defined(S_MP_RAND_PLATFORM_C) #endif +<<<<<<< HEAD #if defined(S_MP_RAND_SOURCE_C) # define S_MP_RAND_PLATFORM_C #endif +||||||| parent of 02af6fb8 (Refinement of and bugfixes in mp_root_n:) +======= +#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 + +>>>>>>> 02af6fb8 (Refinement of and bugfixes in mp_root_n:) #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 be620dbc9..400935ebf 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 From a34e7713aa685d2c34f22a3ce01fd4e3da74608f Mon Sep 17 00:00:00 2001 From: Christoph Zurnieden Date: Sat, 25 Apr 2026 22:35:28 +0200 Subject: [PATCH 2/2] Fixed rebase --- makefile | 6 +++--- makefile.mingw | 6 +++--- makefile.msvc | 6 +++--- makefile.shared | 6 +++--- makefile.unix | 6 +++--- tommath_class.h | 4 ---- 6 files changed, 15 insertions(+), 19 deletions(-) diff --git a/makefile b/makefile index c4c28e491..8e4619cc7 100644 --- a/makefile +++ b/makefile @@ -48,9 +48,9 @@ s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_expt 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_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 +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 dc4f53f2a..b78e12b69 100644 --- a/makefile.mingw +++ b/makefile.mingw @@ -50,9 +50,9 @@ s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_expt 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_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 +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/makefile.msvc b/makefile.msvc index 944ddbb3c..4398633dc 100644 --- a/makefile.msvc +++ b/makefile.msvc @@ -46,9 +46,9 @@ s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj s_mp_exptmod.obj s 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_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 +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 91ca53025..0f7909dc6 100644 --- a/makefile.shared +++ b/makefile.shared @@ -45,9 +45,9 @@ s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_expt 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_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 +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 a0b350eb6..aade2d9d6 100644 --- a/makefile.unix +++ b/makefile.unix @@ -51,9 +51,9 @@ s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_expt 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_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 +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 diff --git a/tommath_class.h b/tommath_class.h index 4531f57af..4c3875792 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -1253,13 +1253,10 @@ #if defined(S_MP_RAND_PLATFORM_C) #endif -<<<<<<< HEAD #if defined(S_MP_RAND_SOURCE_C) # define S_MP_RAND_PLATFORM_C #endif -||||||| parent of 02af6fb8 (Refinement of and bugfixes in mp_root_n:) -======= #if defined(S_MP_ROOT_N_C) # define MP_ADD_C # define MP_ADD_D_C @@ -1282,7 +1279,6 @@ # define S_MP_FP_LOG_C #endif ->>>>>>> 02af6fb8 (Refinement of and bugfixes in mp_root_n:) #if defined(S_MP_SQR_C) # define MP_CLAMP_C # define MP_CLEAR_C