Skip to content

Commit 2a5d8ed

Browse files
author
Christoph Zurnieden
committed
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)
1 parent ae40a87 commit 2a5d8ed

13 files changed

Lines changed: 381 additions & 157 deletions

‎demo/test.c‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1930,6 +1930,26 @@ static int test_mp_root_n(void)
19301930
EXPECT(mp_cmp(&r, &c) == MP_EQ);
19311931
}
19321932
}
1933+
1934+
/* 0^(1/x) = 0 with x != 0 is allowed, test */
1935+
mp_set(&a, 0);
1936+
DO(mp_root_n(&a, 2, &c));
1937+
EXPECT(mp_cmp_d(&c, 0) == MP_EQ);
1938+
1939+
/* Not allowed: division by zero */
1940+
mp_set(&a, 2);
1941+
EXPECT(mp_root_n(&a, 0, &c) == MP_VAL);
1942+
1943+
/* root^base == input with small input and base */
1944+
mp_set(&a, 4);
1945+
DO(mp_root_n(&a, 2, &c));
1946+
EXPECT(mp_cmp_d(&c, 2) == MP_EQ);
1947+
1948+
/* (root^base)^(1/(base + 1)) with small root */
1949+
DO(mp_2expt(&a, 48));
1950+
DO(mp_root_n(&a, 49, &c));
1951+
EXPECT(mp_cmp_d(&c, 1) == MP_EQ);
1952+
19331953
mp_clear_multi(&a, &c, &r, NULL);
19341954
return EXIT_SUCCESS;
19351955
LBL_ERR:

‎libtommath_VS2008.vcproj‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -836,6 +836,10 @@
836836
RelativePath="s_mp_exptmod_fast.c"
837837
>
838838
</File>
839+
<File
840+
RelativePath="s_mp_fp_exp2.c"
841+
>
842+
</File>
839843
<File
840844
RelativePath="s_mp_fp_log.c"
841845
>
@@ -916,6 +920,10 @@
916920
RelativePath="s_mp_rand_source.c"
917921
>
918922
</File>
923+
<File
924+
RelativePath="s_mp_root_n.c"
925+
>
926+
</File>
919927
<File
920928
RelativePath="s_mp_sqr.c"
921929
>

‎makefile‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -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
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 \
4646
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 \
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 \
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 s_mp_sqr.o
5252
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 \
5353
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
5454

55+
5556
#END_INS
5657

5758
$(LIBNAME): $(OBJECTS)

‎makefile.mingw‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -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
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 \
4848
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 \
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 \
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 s_mp_sqr.o \
5454
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 \
5555
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
5656

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

‎makefile.msvc‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,11 +42,11 @@ mp_reduce_setup.obj mp_root_n.obj mp_rshd.obj mp_sbin_size.obj mp_set.obj mp_set
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 \
4444
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 \
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 \
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 s_mp_sqr.obj \
5050
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 \
5151
s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj
5252

‎makefile.shared‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -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
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 \
4343
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 \
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 \
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 s_mp_sqr.o \
4949
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 \
5050
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
5151

52+
5253
#END_INS
5354

5455
objs: $(OBJECTS)

‎makefile.unix‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -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
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 \
4949
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 \
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 \
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 s_mp_sqr.o \
5555
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 \
5656
s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o
5757

5858

59+
5960
HEADERS_PUB=tommath.h
6061
HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)
6162

‎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

‎s_mp_fp_exp2.c‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
#include "tommath_private.h"
2+
#ifdef S_MP_FP_EXP2_C
3+
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
4+
/* SPDX-License-Identifier: Unlicense */
5+
6+
static MP_SET_UNSIGNED(mp_set_word, mp_word)
7+
8+
static mp_err s_mp_fp_exp2_fract(const mp_int *a, mp_int *c)
9+
{
10+
mp_err err = MP_OKAY;
11+
mp_int u, tmp;
12+
/*
13+
Generated with Sam Hocevar's LolRemez at https://github.com/samhocevar/lolremez
14+
15+
lolremez --double -d 10 -r "0:1" "2**(x)"
16+
and
17+
lolremez --double -d 5 -r "0:1" "2**(x)"
18+
19+
*/
20+
#if (MP_DIGIT_BIT == 15)
21+
mp_word coefficients[] = {
22+
0x3e,
23+
0x125,
24+
0x726,
25+
0x1ebc,
26+
0x58b9,
27+
0x7fff
28+
};
29+
#elif ((MP_DIGIT_BIT == 28) || (MP_DIGIT_BIT == 31))
30+
mp_word coefficients[] = {
31+
0x15,
32+
0xca,
33+
0xb2c,
34+
0x7fe0,
35+
0x50c2e,
36+
0x2bb0fc,
37+
0x13b2ab7,
38+
0x71ac235,
39+
0x1ebfbdff,
40+
0x58b90bfb,
41+
0x80000000
42+
};
43+
#elif (MP_DIGIT_BIT == 60)
44+
mp_word coefficients[] = {
45+
0x157d3d6ee0,
46+
0xca97869248,
47+
0xb2c1b72cee4,
48+
0x7fe03bc5603e,
49+
0x50c2e781f2028,
50+
0x2bb0fc474348da,
51+
0x13b2ab7bece61d7,
52+
0x71ac235a89a9728,
53+
0x1ebfbdff845aca73,
54+
0x58b90bfbe8dd8d73,
55+
0x8000000000000acf
56+
};
57+
#endif
58+
int coeff_size = (int)(sizeof(coefficients)/sizeof(coefficients[0]));
59+
int i;
60+
61+
if ((err = mp_init_multi(&u, &tmp, NULL)) != MP_OKAY) {
62+
return err;
63+
}
64+
mp_set_word(&tmp, coefficients[0]);
65+
for (i = 1; i < coeff_size; i++) {
66+
/* u = ((u * x) >> 63) + coefficients[i] */
67+
if ((err = mp_mul(&u, a, &u)) != MP_OKAY) goto LBL_ERR;
68+
if ((err = mp_div_2d(&u, MP_PRECISION_FIXED_LOG, &u, NULL)) != MP_OKAY) goto LBL_ERR;
69+
mp_set_word(&tmp, coefficients[i]);
70+
if ((err = mp_add(&u, &tmp, &u)) != MP_OKAY) goto LBL_ERR;
71+
}
72+
73+
if (c != NULL) {
74+
mp_exch(&u, c);
75+
}
76+
77+
LBL_ERR:
78+
mp_clear_multi(&u, &tmp, NULL);
79+
return err;
80+
}
81+
82+
/* input a is the output from s_mp_fp_log(), a fixed point log base 2 */
83+
mp_err s_mp_fp_exp2(const mp_int *a, mp_int *c)
84+
{
85+
mp_err err = MP_OKAY;
86+
mp_int tmp, int_part, fract_part;
87+
88+
if ((err = mp_init_multi(&tmp, &int_part, &fract_part, NULL)) != MP_OKAY) {
89+
return err;
90+
}
91+
/* Cut fore and aft to get the integer part */
92+
if ((err = mp_div_2d(a, MP_PRECISION_FIXED_LOG, &int_part, NULL)) != MP_OKAY) goto LBL_ERR;
93+
if ((err = mp_mul_2d(&int_part, MP_PRECISION_FIXED_LOG, &int_part)) != MP_OKAY) goto LBL_ERR;
94+
/* Fractional part */
95+
if ((err = mp_sub(a, &int_part, &fract_part)) != MP_OKAY) goto LBL_ERR;
96+
/* 2^{fractional part} */
97+
if ((err = s_mp_fp_exp2_fract(&fract_part, &fract_part)) != MP_OKAY) goto LBL_ERR;
98+
/* A = 1 << normalized integer part */
99+
mp_set(&tmp,1);
100+
if ((err = mp_div_2d(&int_part, MP_PRECISION_FIXED_LOG, &int_part, NULL)) != MP_OKAY) goto LBL_ERR;
101+
if ((err = mp_mul_2d(&tmp, (int)mp_get_l(&int_part), &int_part)) != MP_OKAY) goto LBL_ERR;
102+
/* B = A * fractional part */
103+
if ((err = mp_mul(&int_part, &fract_part, &tmp)) != MP_OKAY) goto LBL_ERR;
104+
/* normalize B */
105+
if ((err = mp_div_2d(&tmp, MP_PRECISION_FIXED_LOG, &tmp, NULL)) != MP_OKAY) goto LBL_ERR;
106+
/* C = 2^a (Roughly. Very, very roughly) */
107+
if (c != NULL) {
108+
mp_exch(&tmp, c);
109+
}
110+
111+
LBL_ERR:
112+
mp_clear_multi(&tmp, &int_part, &fract_part, NULL);
113+
return err;
114+
}
115+
116+
#endif

0 commit comments

Comments
 (0)