Skip to content

Commit 563af5f

Browse files
sjaeckelsjaeckel
committed
also refactor 'check result' part into mp_log()
This code was duplicated as well. Signed-off-by: Steffen Jaeckel <s@jaeckel.eu>
1 parent 6bc7d15 commit 563af5f

File tree

1 file changed

+75
-129
lines changed

1 file changed

+75
-129
lines changed

‎mp_log.c‎

Lines changed: 75 additions & 129 deletions
Original file line numberDiff line numberDiff line change
@@ -7,106 +7,28 @@
77

88
static mp_err s_approx_log_d(const mp_int *a, const mp_int *b, int *lb)
99
{
10-
mp_int bn, La, Lb;
11-
int n;
12-
mp_word la_word, lb_word;
13-
10+
mp_word La, Lb;
1411
mp_err err;
15-
mp_ord cmp;
16-
17-
if ((err = mp_init_multi(&La, &Lb, &bn, NULL)) != MP_OKAY) {
18-
return err;
19-
}
2012

2113
/* Approximation of the individual logarithms with low precision */
22-
if ((err = s_mp_fp_log_d(a, &la_word)) != MP_OKAY) goto LTM_ERR;
23-
if ((err = s_mp_fp_log_d(b, &lb_word)) != MP_OKAY) goto LTM_ERR;
14+
if ((err = s_mp_fp_log_d(a, &La)) != MP_OKAY) goto LTM_ERR;
15+
if ((err = s_mp_fp_log_d(b, &Lb)) != MP_OKAY) goto LTM_ERR;
2416

2517
/* Approximation of log_b(a) with low precision. */
26-
n = (int)(((la_word - (lb_word + 1)/2) / lb_word) + 1);
18+
*lb = (int)(((La - (Lb + 1)/2) / Lb) + 1);
2719
/* TODO: just floor it instead? Multiplication is cheaper than division. */
28-
/* n = (int)(la_word / lb_word); */
29-
30-
/* Check result. Result is wrong by 2(two) at most. */
31-
if ((err = mp_expt_n(b, n, &bn)) != MP_OKAY) {
32-
/* If the approximation overshot we can give it another try */
33-
if (err == MP_OVF) {
34-
n--;
35-
/* But only one */
36-
if ((err = mp_expt_n(b, n, &bn)) != MP_OKAY) goto LTM_ERR;
37-
} else {
38-
goto LTM_ERR;
39-
}
40-
}
41-
42-
cmp = mp_cmp(&bn, a);
43-
44-
/* The rare case of a perfect power makes a perfect shortcut, too. */
45-
if (cmp == MP_EQ) {
46-
*lb = n;
47-
goto LTM_OUT;
48-
}
49-
50-
/* We have to make at least one multiplication because it could still be a perfect power. */
51-
if (cmp == MP_LT) {
52-
do {
53-
/* Full big-integer operations are to be avoided if possible */
54-
if (b->used == 1) {
55-
if ((err = mp_mul_d(&bn, b->dp[0], &bn)) != MP_OKAY) {
56-
if (err == MP_OVF) {
57-
goto LTM_OUT;
58-
}
59-
goto LTM_ERR;
60-
}
61-
} else {
62-
if ((err = mp_mul(&bn, b, &bn)) != MP_OKAY) {
63-
if (err == MP_OVF) {
64-
goto LTM_OUT;
65-
}
66-
goto LTM_ERR;
67-
}
68-
}
69-
n++;
70-
} while ((cmp = mp_cmp(&bn, a)) == MP_LT);
71-
/* Overshot, take it back. */
72-
if (cmp == MP_GT) {
73-
n--;
74-
}
75-
goto LTM_OUT;
76-
}
77-
78-
/* But it can overestimate, too, for example if "a" is closely below some "b^k" */
79-
if (cmp == MP_GT) {
80-
do {
81-
if (b->used == 1) {
82-
/* These are cheaper exact divisions, but that function is not available in LTM */
83-
if ((err = mp_div_d(&bn, b->dp[0], &bn, NULL)) != MP_OKAY) goto LTM_ERR;
84-
85-
} else {
86-
if ((err = mp_div(&bn, b, &bn, NULL)) != MP_OKAY) goto LTM_ERR;
87-
}
88-
n--;
89-
} while ((cmp = mp_cmp(&bn, a)) == MP_GT);
90-
}
91-
92-
LTM_OUT:
93-
*lb = n;
20+
/* *lb = (int)(la_word / lb_word); */
9421
err = MP_OKAY;
9522
LTM_ERR:
96-
mp_clear_multi(&La, &Lb, &bn, NULL);
9723
return err;
9824
}
9925

10026
static mp_err s_approx_log(const mp_int *a, const mp_int *b, int *lb)
10127
{
102-
mp_int bn, La, Lb, t;
103-
104-
int n;
105-
28+
mp_int La, Lb, t;
10629
mp_err err;
107-
mp_ord cmp;
10830

109-
if ((err = mp_init_multi(&La, &Lb, &bn, &t, NULL)) != MP_OKAY) {
31+
if ((err = mp_init_multi(&La, &Lb, &t, NULL)) != MP_OKAY) {
11032
return err;
11133
}
11234

@@ -119,25 +41,86 @@ static mp_err s_approx_log(const mp_int *a, const mp_int *b, int *lb)
11941
if ((err = mp_div(&t, &Lb, &t, NULL)) != MP_OKAY) goto LTM_ERR;
12042
if ((err = mp_add_d(&t, 1u, &t)) != MP_OKAY) goto LTM_ERR;
12143

122-
n = mp_get_i32(&t);
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+
51+
52+
mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
53+
{
54+
mp_int bn;
55+
int n, fla, flb;
56+
mp_err err;
57+
mp_ord cmp;
58+
59+
if (mp_isneg(a) || mp_iszero(a) || (mp_cmp_d(b, 2u) == MP_LT)) {
60+
return MP_VAL;
61+
}
62+
63+
if (MP_IS_POWER_OF_TWO(b)) {
64+
*lb = MP_LOG2_EXPT(a, b);
65+
return MP_OKAY;
66+
}
67+
68+
/* floor(log_2(x)) for cut-off */
69+
fla = mp_count_bits(a) - 1;
70+
flb = mp_count_bits(b) - 1;
71+
72+
cmp = mp_cmp(a, b);
73+
74+
/* "a < b -> 0" and "(b == a) -> 1" */
75+
if ((cmp == MP_LT) || (cmp == MP_EQ)) {
76+
*lb = (cmp == MP_EQ);
77+
return MP_OKAY;
78+
}
79+
80+
/* "a < b^2 -> 1" (bit-count is sufficient, doesn't need to be exact) */
81+
if (((2 * flb)-1) > fla) {
82+
*lb = 1;
83+
return MP_OKAY;
84+
}
85+
86+
87+
if (MP_HAS(S_MP_LOG_D_POSSIBLE)) {
88+
err = s_approx_log_d(a, b, &n);
89+
} else {
90+
err = s_approx_log(a, b, &n);
91+
}
92+
if (err != MP_OKAY) {
93+
return err;
94+
}
12395

96+
if ((err = mp_init(&bn)) != MP_OKAY) {
97+
return err;
98+
}
99+
100+
/* Check result. Result is wrong by 2(two) at most. */
124101
if ((err = mp_expt_n(b, n, &bn)) != MP_OKAY) {
102+
/* If the approximation overshot we can give it another try */
125103
if (err == MP_OVF) {
126104
n--;
105+
/* But only one */
127106
if ((err = mp_expt_n(b, n, &bn)) != MP_OKAY) goto LTM_ERR;
128107
} else {
129108
goto LTM_ERR;
130109
}
131110
}
111+
132112
cmp = mp_cmp(&bn, a);
133113

114+
/* The rare case of a perfect power makes a perfect shortcut, too. */
134115
if (cmp == MP_EQ) {
135116
*lb = n;
136117
goto LTM_OUT;
137118
}
138119

120+
/* We have to make at least one multiplication because it could still be a perfect power. */
139121
if (cmp == MP_LT) {
140122
do {
123+
/* Full big-integer operations are to be avoided if possible */
141124
if (b->used == 1) {
142125
if ((err = mp_mul_d(&bn, b->dp[0], &bn)) != MP_OKAY) {
143126
if (err == MP_OVF) {
@@ -155,19 +138,22 @@ static mp_err s_approx_log(const mp_int *a, const mp_int *b, int *lb)
155138
}
156139
n++;
157140
} while ((cmp = mp_cmp(&bn, a)) == MP_LT);
141+
/* Overshot, take it back. */
158142
if (cmp == MP_GT) {
159143
n--;
160144
}
161145
goto LTM_OUT;
162146
}
163147

148+
/* But it can overestimate, too, for example if "a" is closely below some "b^k" */
164149
if (cmp == MP_GT) {
165150
do {
166151
if (b->used == 1) {
167-
if ((err = mp_div_d(&bn, b->dp[0], &bn, NULL)) != MP_OKAY) goto LTM_ERR;
152+
/* These are cheaper exact divisions, but that function is not available in LTM */
153+
if ((err = mp_div_d(&bn, b->dp[0], &bn, NULL)) != MP_OKAY) goto LTM_ERR;
168154

169155
} else {
170-
if ((err = mp_div(&bn, b, &bn, NULL)) != MP_OKAY) goto LTM_ERR;
156+
if ((err = mp_div(&bn, b, &bn, NULL)) != MP_OKAY) goto LTM_ERR;
171157
}
172158
n--;
173159
} while ((cmp = mp_cmp(&bn, a)) == MP_GT);
@@ -177,48 +163,8 @@ static mp_err s_approx_log(const mp_int *a, const mp_int *b, int *lb)
177163
*lb = n;
178164
err = MP_OKAY;
179165
LTM_ERR:
180-
mp_clear_multi(&La, &Lb, &bn, &t, NULL);
166+
mp_clear(&bn);
181167
return err;
182168
}
183169

184-
185-
mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
186-
{
187-
int fla, flb;
188-
mp_ord cmp;
189-
190-
if (mp_isneg(a) || mp_iszero(a) || (mp_cmp_d(b, 2u) == MP_LT)) {
191-
return MP_VAL;
192-
}
193-
194-
if (MP_IS_POWER_OF_TWO(b)) {
195-
*lb = MP_LOG2_EXPT(a, b);
196-
return MP_OKAY;
197-
}
198-
199-
/* floor(log_2(x)) for cut-off */
200-
fla = mp_count_bits(a) - 1;
201-
flb = mp_count_bits(b) - 1;
202-
203-
cmp = mp_cmp(a, b);
204-
205-
/* "a < b -> 0" and "(b == a) -> 1" */
206-
if ((cmp == MP_LT) || (cmp == MP_EQ)) {
207-
*lb = (cmp == MP_EQ);
208-
return MP_OKAY;
209-
}
210-
211-
/* "a < b^2 -> 1" (bit-count is sufficient, doesn't need to be exact) */
212-
if (((2 * flb)-1) > fla) {
213-
*lb = 1;
214-
return MP_OKAY;
215-
}
216-
217-
218-
if (MP_HAS(S_MP_LOG_D_POSSIBLE)) {
219-
return s_approx_log_d(a, b, lb);
220-
}
221-
return s_approx_log(a, b, lb);
222-
}
223-
224170
#endif

0 commit comments

Comments
 (0)