33/* LibTomMath, multiple-precision integer library -- Tom St Denis */
44/* SPDX-License-Identifier: Unlicense */
55
6-
76#define MP_LOG2_EXPT (a ,b ) ((mp_count_bits((a)) - 1) / mp_cnt_lsb((b)))
87
9- /*
10- This functions rely on the size of mp_word being larger than INT_MAX and in case
11- there is a really weird architecture we try to check for it. Not a 100% reliable
12- test but it has a safe fallback.
13- */
14- #if ( (UINT_MAX == UINT32_MAX ) && ( MP_WORD_SIZE > 4 ) ) \
15- || \
16- ( (UINT_MAX == UINT16_MAX ) && ( MP_WORD_SIZE > 2 ) )
17-
18- mp_err mp_log (const mp_int * a , const mp_int * b , int * lb )
8+ static mp_err s_approx_log_d (const mp_int * a , const mp_int * b , int * lb )
199{
2010 mp_int bn , La , Lb ;
21- int n , fla , flb ;
11+ int n ;
2212 mp_word la_word , lb_word ;
2313
2414 mp_err err ;
2515 mp_ord cmp ;
2616
27- if (mp_isneg (a ) || mp_iszero (a ) || (mp_cmp_d (b , 2u ) == MP_LT )) {
28- return MP_VAL ;
29- }
30-
31- if (MP_HAS (S_MP_LOG_2EXPT ) && MP_IS_POWER_OF_TWO (b )) {
32- * lb = MP_LOG2_EXPT (a , b );
33- return MP_OKAY ;
34- }
35-
36- /* floor(log_2(x)) for cut-off */
37- fla = mp_count_bits (a ) - 1 ;
38- flb = mp_count_bits (b ) - 1 ;
39-
40- cmp = mp_cmp (a , b );
41-
42- /* "a < b -> 0" and "(b == a) -> 1" */
43- if ((cmp == MP_LT ) || (cmp == MP_EQ )) {
44- * lb = (cmp == MP_EQ );
45- return MP_OKAY ;
46- }
47-
48- /* "a < b^2 -> 1" (bit-count is sufficient, doesn't need to be exact) */
49- if (((2 * flb )- 1 ) > fla ) {
50- * lb = 1 ;
51- return MP_OKAY ;
52- }
53-
5417 if ((err = mp_init_multi (& La , & Lb , & bn , NULL )) != MP_OKAY ) {
5518 return err ;
5619 }
5720
5821 /* Approximation of the individual logarithms with low precision */
59- if ((err = s_mp_fp_log (a , & la_word )) != MP_OKAY ) goto LTM_ERR ;
60- if ((err = s_mp_fp_log (b , & lb_word )) != MP_OKAY ) goto LTM_ERR ;
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 ;
6124
6225 /* Approximation of log_b(a) with low precision. */
6326 n = (int )(((la_word - (lb_word + 1 )/2 ) / lb_word ) + 1 );
@@ -117,10 +80,10 @@ mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
11780 do {
11881 if (b -> used == 1 ) {
11982 /* These are cheaper exact divisions, but that function is not available in LTM */
120- if ((err = mp_div_d (& bn , b -> dp [0 ], & bn , NULL )) != MP_OKAY ) goto LTM_ERR ;
83+ if ((err = mp_div_d (& bn , b -> dp [0 ], & bn , NULL )) != MP_OKAY ) goto LTM_ERR ;
12184
12285 } else {
123- if ((err = mp_div (& bn , b , & bn , NULL )) != MP_OKAY ) goto LTM_ERR ;
86+ if ((err = mp_div (& bn , b , & bn , NULL )) != MP_OKAY ) goto LTM_ERR ;
12487 }
12588 n -- ;
12689 } while ((cmp = mp_cmp (& bn , a )) == MP_GT );
@@ -134,61 +97,34 @@ mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
13497 return err ;
13598}
13699
137-
138- #else
139- /* Bigint version. A bit slower but not _that_ much */
140- mp_err mp_log (const mp_int * a , const mp_int * b , int * lb )
100+ static mp_err s_approx_log (const mp_int * a , const mp_int * b , int * lb )
141101{
142102 mp_int bn , La , Lb , t ;
143103
144- int n , fla , flb ;
104+ int n ;
145105
146106 mp_err err ;
147107 mp_ord cmp ;
148108
149- if (mp_isneg (a ) || mp_iszero (a ) || (mp_cmp_d (b , 2u ) == MP_LT )) {
150- return MP_VAL ;
151- }
152-
153- if (MP_HAS (S_MP_LOG_2EXPT ) && MP_IS_POWER_OF_TWO (b )) {
154- * lb = MP_LOG2_EXPT (a , b );
155- return MP_OKAY ;
156- }
157-
158- fla = mp_count_bits (a ) - 1 ;
159- flb = mp_count_bits (b ) - 1 ;
160-
161- cmp = mp_cmp (a , b );
162-
163- if ((cmp == MP_LT ) || (cmp == MP_EQ )) {
164- * lb = (cmp == MP_EQ );
165- return MP_OKAY ;
166- }
167-
168- if ((2 * flb ) > fla ) {
169- * lb = 1 ;
170- return MP_OKAY ;
171- }
172-
173109 if ((err = mp_init_multi (& La , & Lb , & bn , & t , NULL )) != MP_OKAY ) {
174110 return err ;
175111 }
176112
177- if ((err = s_mp_fp_log (a , & La )) != MP_OKAY ) goto LTM_ERR ;
178- if ((err = s_mp_fp_log (b , & Lb )) != MP_OKAY ) goto LTM_ERR ;
113+ if ((err = s_mp_fp_log (a , & La )) != MP_OKAY ) goto LTM_ERR ;
114+ if ((err = s_mp_fp_log (b , & Lb )) != MP_OKAY ) goto LTM_ERR ;
179115
180- if ((err = mp_add_d (& Lb , 1u , & t )) != MP_OKAY ) goto LTM_ERR ;
181- if ((err = mp_div_2 (& t , & t )) != MP_OKAY ) goto LTM_ERR ;
182- if ((err = mp_sub (& La , & t , & t )) != MP_OKAY ) goto LTM_ERR ;
183- if ((err = mp_div (& t , & Lb , & t , NULL )) != MP_OKAY ) goto LTM_ERR ;
184- if ((err = mp_add_d (& t , 1u , & t )) != MP_OKAY ) goto LTM_ERR ;
116+ if ((err = mp_add_d (& Lb , 1u , & t )) != MP_OKAY ) goto LTM_ERR ;
117+ if ((err = mp_div_2 (& t , & t )) != MP_OKAY ) goto LTM_ERR ;
118+ if ((err = mp_sub (& La , & t , & t )) != MP_OKAY ) goto LTM_ERR ;
119+ if ((err = mp_div (& t , & Lb , & t , NULL )) != MP_OKAY ) goto LTM_ERR ;
120+ if ((err = mp_add_d (& t , 1u , & t )) != MP_OKAY ) goto LTM_ERR ;
185121
186122 n = mp_get_i32 (& t );
187123
188124 if ((err = mp_expt_n (b , n , & bn )) != MP_OKAY ) {
189125 if (err == MP_OVF ) {
190126 n -- ;
191- if ((err = mp_expt_n (b , n , & bn )) != MP_OKAY ) goto LTM_ERR ;
127+ if ((err = mp_expt_n (b , n , & bn )) != MP_OKAY ) goto LTM_ERR ;
192128 } else {
193129 goto LTM_ERR ;
194130 }
@@ -244,8 +180,45 @@ mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
244180 mp_clear_multi (& La , & Lb , & bn , & t , NULL );
245181 return err ;
246182}
247- #endif
248183
249184
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+ }
250223
251224#endif
0 commit comments