Skip to content

Commit a1caf2c

Browse files
author
Christoph Zurnieden
committed
Refinement of and bugfixes in mp_root_n:
 - Changed all fixed point functions to use uint32_t instead of mp_word. - Implementation of a fixed point exp2m1 function to get better start-values for the nth-root function - fixed treatment of perfect powers - fixed 0^(1/x) == 0 with x != 0 - additional test for b == 0 (division by zero)
1 parent ae40a87 commit a1caf2c

18 files changed

Lines changed: 528 additions & 264 deletions

‎demo/test.c‎

Lines changed: 22 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1454,12 +1454,6 @@ static int test_mp_log_n(void)
14541454
int base, lb, size, i;
14551455
const int max_base = MP_MIN(INT_MAX, MP_DIGIT_MAX);
14561456

1457-
if (MP_HAS(S_MP_WORD_TOO_SMALL)) {
1458-
fprintf(stderr, "Testing mp_log_n with restricted size of mp_word.\n");
1459-
} else {
1460-
fprintf(stderr, "Testing mp_log_n with normal size of mp_word.\n");
1461-
}
1462-
14631457
DOR(mp_init(&a));
14641458

14651459
/*
@@ -1930,6 +1924,26 @@ static int test_mp_root_n(void)
19301924
EXPECT(mp_cmp(&r, &c) == MP_EQ);
19311925
}
19321926
}
1927+
1928+
/* 0^(1/x) = 0 with x != 0 is allowed, test */
1929+
mp_set(&a, 0);
1930+
DO(mp_root_n(&a, 2, &c));
1931+
EXPECT(mp_cmp_d(&c, 0) == MP_EQ);
1932+
1933+
/* Not allowed: division by zero */
1934+
mp_set(&a, 2);
1935+
EXPECT(mp_root_n(&a, 0, &c) == MP_VAL);
1936+
1937+
/* root^base == input with small input and base */
1938+
mp_set(&a, 4);
1939+
DO(mp_root_n(&a, 2, &c));
1940+
EXPECT(mp_cmp_d(&c, 2) == MP_EQ);
1941+
1942+
/* (root^base)^(1/(base + 1)) with small root */
1943+
DO(mp_2expt(&a, 48));
1944+
DO(mp_root_n(&a, 49, &c));
1945+
EXPECT(mp_cmp_d(&c, 1) == MP_EQ);
1946+
19331947
mp_clear_multi(&a, &c, &r, NULL);
19341948
return EXIT_SUCCESS;
19351949
LBL_ERR:
@@ -2647,8 +2661,9 @@ static int unit_tests(int argc, char **argv)
26472661
if (j == argc) continue;
26482662
}
26492663

2650-
if (test[i].fn)
2664+
if (test[i].fn) {
26512665
j = test[i].fn();
2666+
}
26522667
} else if (MP_HAS(MULTI_THREADED)) {
26532668
EXPECT(thread_join(&test_threads[i], &res) == 0);
26542669
j = res->ret;

‎libtommath_VS2008.vcproj‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -804,6 +804,10 @@
804804
RelativePath="mp_zero.c"
805805
>
806806
</File>
807+
<File
808+
RelativePath="s_mp_32_umul32.c"
809+
>
810+
</File>
807811
<File
808812
RelativePath="s_mp_add.c"
809813
>
@@ -836,6 +840,10 @@
836840
RelativePath="s_mp_exptmod_fast.c"
837841
>
838842
</File>
843+
<File
844+
RelativePath="s_mp_fp_exp2.c"
845+
>
846+
</File>
839847
<File
840848
RelativePath="s_mp_fp_log.c"
841849
>
@@ -916,6 +924,14 @@
916924
RelativePath="s_mp_rand_source.c"
917925
>
918926
</File>
927+
<File
928+
RelativePath="s_mp_root_n.c"
929+
>
930+
</File>
931+
<File
932+
RelativePath="s_mp_set_word.c"
933+
>
934+
</File>
919935
<File
920936
RelativePath="s_mp_sqr.c"
921937
>

‎makefile‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -43,14 +43,15 @@ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o m
4343
mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \
4444
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 \
4545
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 \
46-
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 \
47-
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 \
48-
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 \
46+
mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_32_umul32.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
47+
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 \
48+
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 \
4949
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 \
5050
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 \
51-
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 \
52-
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 \
53-
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
51+
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 \
52+
s_mp_set_word.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 \
53+
s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
54+
5455

5556
#END_INS
5657

‎makefile.mingw‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -45,14 +45,15 @@ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o m
4545
mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \
4646
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 \
4747
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 \
48-
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 \
49-
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 \
50-
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 \
48+
mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_32_umul32.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
49+
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 \
50+
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 \
5151
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 \
5252
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 \
53-
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 \
54-
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 \
55-
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
53+
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 \
54+
s_mp_set_word.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 \
55+
s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
56+
5657

5758
HEADERS_PUB=tommath.h
5859
HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)

‎makefile.msvc‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -41,14 +41,14 @@ mp_reduce_2k_l.obj mp_reduce_2k_setup.obj mp_reduce_2k_setup_l.obj mp_reduce_is_
4141
mp_reduce_setup.obj mp_root_n.obj mp_rshd.obj mp_sbin_size.obj mp_set.obj mp_set_double.obj mp_set_i32.obj mp_set_i64.obj \
4242
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 \
4343
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 \
44-
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 \
45-
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 \
46-
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 \
44+
mp_unpack.obj mp_warray_free.obj mp_xor.obj mp_zero.obj s_mp_32_umul32.obj s_mp_add.obj s_mp_copy_digs.obj s_mp_div_3.obj \
45+
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 \
46+
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 \
4747
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 \
4848
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 \
49-
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 \
50-
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 \
51-
s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj
49+
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 \
50+
s_mp_set_word.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 \
51+
s_mp_warray_get.obj s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj
5252

5353
HEADERS_PUB=tommath.h
5454
HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)

‎makefile.shared‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,14 +40,15 @@ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o m
4040
mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \
4141
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 \
4242
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 \
43-
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 \
44-
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 \
45-
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 \
43+
mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_32_umul32.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
44+
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 \
45+
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 \
4646
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 \
4747
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 \
48-
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 \
49-
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 \
50-
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
48+
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 \
49+
s_mp_set_word.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 \
50+
s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
51+
5152

5253
#END_INS
5354

‎makefile.unix‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -46,14 +46,15 @@ mp_reduce_2k_l.o mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o m
4646
mp_reduce_setup.o mp_root_n.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o \
4747
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 \
4848
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 \
49-
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 \
50-
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 \
51-
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 \
49+
mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_32_umul32.o s_mp_add.o s_mp_copy_digs.o s_mp_div_3.o \
50+
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 \
51+
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 \
5252
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 \
5353
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 \
54-
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 \
55-
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 \
56-
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
54+
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 \
55+
s_mp_set_word.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 \
56+
s_mp_warray_get.o s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
57+
5758

5859

5960
HEADERS_PUB=tommath.h

‎mp_log.c‎

Lines changed: 2 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -23,31 +23,6 @@ static mp_err s_approx_log_d(const mp_int *a, const mp_int *b, int *lb)
2323
return err;
2424
}
2525

26-
static mp_err s_approx_log(const mp_int *a, const mp_int *b, int *lb)
27-
{
28-
mp_int La, Lb, t;
29-
mp_err err;
30-
31-
if ((err = mp_init_multi(&La, &Lb, &t, NULL)) != MP_OKAY) {
32-
return err;
33-
}
34-
35-
if ((err = s_mp_fp_log(a, &La)) != MP_OKAY) goto LTM_ERR;
36-
if ((err = s_mp_fp_log(b, &Lb)) != MP_OKAY) goto LTM_ERR;
37-
38-
if ((err = mp_add_d(&Lb, 1u, &t)) != MP_OKAY) goto LTM_ERR;
39-
if ((err = mp_div_2(&t, &t)) != MP_OKAY) goto LTM_ERR;
40-
if ((err = mp_sub(&La, &t, &t)) != MP_OKAY) goto LTM_ERR;
41-
if ((err = mp_div(&t, &Lb, &t, NULL)) != MP_OKAY) goto LTM_ERR;
42-
if ((err = mp_add_d(&t, 1u, &t)) != MP_OKAY) goto LTM_ERR;
43-
44-
*lb = mp_get_i32(&t);
45-
err = MP_OKAY;
46-
LTM_ERR:
47-
mp_clear_multi(&t, &Lb, &La, NULL);
48-
return err;
49-
}
50-
5126
mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
5227
{
5328
mp_int bn;
@@ -82,11 +57,8 @@ mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
8257
return MP_OKAY;
8358
}
8459

85-
if (MP_HAS(S_MP_WORD_TOO_SMALL)) {
86-
err = s_approx_log(a, b, &n);
87-
} else {
88-
err = s_approx_log_d(a, b, &n);
89-
}
60+
err = s_approx_log_d(a, b, &n);
61+
9062
if (err != MP_OKAY) {
9163
return err;
9264
}

‎mp_root_n.c‎

Lines changed: 2 additions & 127 deletions
Original file line numberDiff line numberDiff line change
@@ -6,136 +6,11 @@
66
/* find the n'th root of an integer
77
*
88
* Result found such that (c)**b <= a and (c+1)**b > a
9-
*
10-
* This algorithm uses Newton's approximation
11-
* x[i+1] = x[i] - f(x[i])/f'(x[i])
12-
* which will find the root in log(N) time where
13-
* each step involves a fair bit.
149
*/
10+
1511
mp_err mp_root_n(const mp_int *a, int b, mp_int *c)
1612
{
17-
mp_int t1, t2, t3, a_;
18-
int ilog2;
19-
mp_err err;
20-
21-
if (b < 0 || (unsigned)b > (unsigned)MP_DIGIT_MAX) {
22-
return MP_VAL;
23-
}
24-
25-
/* input must be positive if b is even */
26-
if (((b & 1) == 0) && mp_isneg(a)) {
27-
return MP_VAL;
28-
}
29-
30-
if ((err = mp_init_multi(&t1, &t2, &t3, NULL)) != MP_OKAY) {
31-
return err;
32-
}
33-
34-
/* if a is negative fudge the sign but keep track */
35-
a_ = *a;
36-
a_.sign = MP_ZPOS;
37-
38-
/* Compute seed: 2^(log_2(n)/b + 2)*/
39-
ilog2 = mp_count_bits(a);
40-
41-
/*
42-
If "b" is larger than INT_MAX it is also larger than
43-
log_2(n) because the bit-length of the "n" is measured
44-
with an int and hence the root is always < 2 (two).
45-
*/
46-
if (b > INT_MAX/2) {
47-
mp_set(c, 1uL);
48-
c->sign = a->sign;
49-
err = MP_OKAY;
50-
goto LBL_ERR;
51-
}
52-
53-
/* "b" is smaller than INT_MAX, we can cast safely */
54-
if (ilog2 < b) {
55-
mp_set(c, 1uL);
56-
c->sign = a->sign;
57-
err = MP_OKAY;
58-
goto LBL_ERR;
59-
}
60-
ilog2 = ilog2 / b;
61-
if (ilog2 == 0) {
62-
mp_set(c, 1uL);
63-
c->sign = a->sign;
64-
err = MP_OKAY;
65-
goto LBL_ERR;
66-
}
67-
/* Start value must be larger than root */
68-
ilog2 += 2;
69-
if ((err = mp_2expt(&t2,ilog2)) != MP_OKAY) goto LBL_ERR;
70-
do {
71-
/* t1 = t2 */
72-
if ((err = mp_copy(&t2, &t1)) != MP_OKAY) goto LBL_ERR;
73-
74-
/* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
75-
76-
/* t3 = t1**(b-1) */
77-
if ((err = mp_expt_n(&t1, b - 1, &t3)) != MP_OKAY) goto LBL_ERR;
78-
79-
/* numerator */
80-
/* t2 = t1**b */
81-
if ((err = mp_mul(&t3, &t1, &t2)) != MP_OKAY) goto LBL_ERR;
82-
83-
/* t2 = t1**b - a */
84-
if ((err = mp_sub(&t2, &a_, &t2)) != MP_OKAY) goto LBL_ERR;
85-
86-
/* denominator */
87-
/* t3 = t1**(b-1) * b */
88-
if ((err = mp_mul_d(&t3, (mp_digit)b, &t3)) != MP_OKAY) goto LBL_ERR;
89-
90-
/* t3 = (t1**b - a)/(b * t1**(b-1)) */
91-
if ((err = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY) goto LBL_ERR;
92-
93-
if ((err = mp_sub(&t1, &t3, &t2)) != MP_OKAY) goto LBL_ERR;
94-
95-
/*
96-
Number of rounds is at most log_2(root). If it is more it
97-
got stuck, so break out of the loop and do the rest manually.
98-
*/
99-
if (ilog2-- == 0) {
100-
break;
101-
}
102-
} while (mp_cmp(&t1, &t2) != MP_EQ);
103-
104-
/* result can be off by a few so check */
105-
/* Loop beneath can overshoot by one if found root is smaller than actual root */
106-
for (;;) {
107-
mp_ord cmp;
108-
if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR;
109-
cmp = mp_cmp(&t2, &a_);
110-
if (cmp == MP_EQ) {
111-
err = MP_OKAY;
112-
goto LBL_ERR;
113-
}
114-
if (cmp == MP_LT) {
115-
if ((err = mp_add_d(&t1, 1uL, &t1)) != MP_OKAY) goto LBL_ERR;
116-
} else {
117-
break;
118-
}
119-
}
120-
/* correct overshoot from above or from recurrence */
121-
for (;;) {
122-
if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR;
123-
if (mp_cmp(&t2, &a_) == MP_GT) {
124-
if ((err = mp_sub_d(&t1, 1uL, &t1)) != MP_OKAY) goto LBL_ERR;
125-
} else {
126-
break;
127-
}
128-
}
129-
130-
/* set the result */
131-
mp_exch(&t1, c);
132-
133-
/* set the sign of the result */
134-
c->sign = a->sign;
135-
136-
LBL_ERR:
137-
mp_clear_multi(&t1, &t2, &t3, NULL);
138-
return err;
13+
return s_mp_root_n(a, b, c, NULL);
13914
}
14015

14116
#endif

0 commit comments

Comments
 (0)