Mercurial > dropbear

comparison libtommath/bn_mp_log_u32.c @ 1692:1051e4eea25a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update LibTomMath to 1.2.0 (#84) * update C files * update other files * update headers * update makefiles * remove mp_set/get_double() * use ltm 1.2.0 API * update ltm_desc * use bundled tommath if system-tommath is too old * XMALLOC etc. were changed to MP_MALLOC etc.
author Steffen Jaeckel <s@jaeckel.eu>
date Tue, 26 May 2020 17:36:47 +0200
parents
children
comparison
equal deleted inserted replaced
1691:2d3745d58843 1692:1051e4eea25a
1 #include "tommath_private.h"
2 #ifdef BN_MP_LOG_U32_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
4 /* SPDX-License-Identifier: Unlicense */
5
6 /* Compute log_{base}(a) */
7 static mp_word s_pow(mp_word base, mp_word exponent)
8 {
9 mp_word result = 1uLL;
10 while (exponent != 0u) {
11 if ((exponent & 1u) == 1u) {
12 result *= base;
13 }
14 exponent >>= 1;
15 base *= base;
16 }
17
18 return result;
19 }
20
21 static mp_digit s_digit_ilogb(mp_digit base, mp_digit n)
22 {
23 mp_word bracket_low = 1uLL, bracket_mid, bracket_high, N;
24 mp_digit ret, high = 1uL, low = 0uL, mid;
25
26 if (n < base) {
27 return 0uL;
28 }
29 if (n == base) {
30 return 1uL;
31 }
32
33 bracket_high = (mp_word) base ;
34 N = (mp_word) n;
35
36 while (bracket_high < N) {
37 low = high;
38 bracket_low = bracket_high;
39 high <<= 1;
40 bracket_high *= bracket_high;
41 }
42
43 while (((mp_digit)(high - low)) > 1uL) {
44 mid = (low + high) >> 1;
45 bracket_mid = bracket_low * s_pow(base, (mp_word)(mid - low));
46
47 if (N < bracket_mid) {
48 high = mid ;
49 bracket_high = bracket_mid ;
50 }
51 if (N > bracket_mid) {
52 low = mid ;
53 bracket_low = bracket_mid ;
54 }
55 if (N == bracket_mid) {
56 return (mp_digit) mid;
57 }
58 }
59
60 if (bracket_high == N) {
61 ret = high;
62 } else {
63 ret = low;
64 }
65
66 return ret;
67 }
68
69 /* TODO: output could be "int" because the output of mp_radix_size is int, too,
70 as is the output of mp_bitcount.
71 With the same problem: max size is INT_MAX * MP_DIGIT not INT_MAX only!
72 */
73 mp_err mp_log_u32(const mp_int *a, uint32_t base, uint32_t *c)
74 {
75 mp_err err;
76 mp_ord cmp;
77 uint32_t high, low, mid;
78 mp_int bracket_low, bracket_high, bracket_mid, t, bi_base;
79
80 err = MP_OKAY;
81
82 if (a->sign == MP_NEG) {
83 return MP_VAL;
84 }
85
86 if (MP_IS_ZERO(a)) {
87 return MP_VAL;
88 }
89
90 if (base < 2u) {
91 return MP_VAL;
92 }
93
94 /* A small shortcut for bases that are powers of two. */
95 if ((base & (base - 1u)) == 0u) {
96 int y, bit_count;
97 for (y=0; (y < 7) && ((base & 1u) == 0u); y++) {
98 base >>= 1;
99 }
100 bit_count = mp_count_bits(a) - 1;
101 *c = (uint32_t)(bit_count/y);
102 return MP_OKAY;
103 }
104
105 if (a->used == 1) {
106 *c = (uint32_t)s_digit_ilogb(base, a->dp[0]);
107 return err;
108 }
109
110 cmp = mp_cmp_d(a, base);
111 if ((cmp == MP_LT) || (cmp == MP_EQ)) {
112 *c = cmp == MP_EQ;
113 return err;
114 }
115
116 if ((err =
117 mp_init_multi(&bracket_low, &bracket_high,
118 &bracket_mid, &t, &bi_base, NULL)) != MP_OKAY) {
119 return err;
120 }
121
122 low = 0u;
123 mp_set(&bracket_low, 1uL);
124 high = 1u;
125
126 mp_set(&bracket_high, base);
127
128 /*
129 A kind of Giant-step/baby-step algorithm.
130 Idea shamelessly stolen from https://programmingpraxis.com/2010/05/07/integer-logarithms/2/
131 The effect is asymptotic, hence needs benchmarks to test if the Giant-step should be skipped
132 for small n.
133 */
134 while (mp_cmp(&bracket_high, a) == MP_LT) {
135 low = high;
136 if ((err = mp_copy(&bracket_high, &bracket_low)) != MP_OKAY) {
137 goto LBL_ERR;
138 }
139 high <<= 1;
140 if ((err = mp_sqr(&bracket_high, &bracket_high)) != MP_OKAY) {
141 goto LBL_ERR;
142 }
143 }
144 mp_set(&bi_base, base);
145
146 while ((high - low) > 1u) {
147 mid = (high + low) >> 1;
148
149 if ((err = mp_expt_u32(&bi_base, (uint32_t)(mid - low), &t)) != MP_OKAY) {
150 goto LBL_ERR;
151 }
152 if ((err = mp_mul(&bracket_low, &t, &bracket_mid)) != MP_OKAY) {
153 goto LBL_ERR;
154 }
155 cmp = mp_cmp(a, &bracket_mid);
156 if (cmp == MP_LT) {
157 high = mid;
158 mp_exch(&bracket_mid, &bracket_high);
159 }
160 if (cmp == MP_GT) {
161 low = mid;
162 mp_exch(&bracket_mid, &bracket_low);
163 }
164 if (cmp == MP_EQ) {
165 *c = mid;
166 goto LBL_END;
167 }
168 }
169
170 *c = (mp_cmp(&bracket_high, a) == MP_EQ) ? high : low;
171
172 LBL_END:
173 LBL_ERR:
174 mp_clear_multi(&bracket_low, &bracket_high, &bracket_mid,
175 &t, &bi_base, NULL);
176 return err;
177 }
178
179
180 #endif