/src/botan/src/lib/math/numbertheory/numthry.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | * Number Theory Functions |
3 | | * (C) 1999-2011,2016,2018,2019 Jack Lloyd |
4 | | * (C) 2007,2008 Falko Strenzke, FlexSecure GmbH |
5 | | * |
6 | | * Botan is released under the Simplified BSD License (see license.txt) |
7 | | */ |
8 | | |
9 | | #include <botan/numthry.h> |
10 | | |
11 | | #include <botan/reducer.h> |
12 | | #include <botan/rng.h> |
13 | | #include <botan/internal/ct_utils.h> |
14 | | #include <botan/internal/divide.h> |
15 | | #include <botan/internal/monty.h> |
16 | | #include <botan/internal/monty_exp.h> |
17 | | #include <botan/internal/mp_core.h> |
18 | | #include <botan/internal/primality.h> |
19 | | #include <algorithm> |
20 | | |
21 | | namespace Botan { |
22 | | |
23 | | namespace { |
24 | | |
25 | 0 | void sub_abs(BigInt& z, const BigInt& x, const BigInt& y) { |
26 | 0 | const size_t x_sw = x.sig_words(); |
27 | 0 | const size_t y_sw = y.sig_words(); |
28 | 0 | z.resize(std::max(x_sw, y_sw)); |
29 | |
|
30 | 0 | bigint_sub_abs(z.mutable_data(), x.data(), x_sw, y.data(), y_sw); |
31 | 0 | } |
32 | | |
33 | | } // namespace |
34 | | |
35 | | /* |
36 | | * Tonelli-Shanks algorithm |
37 | | */ |
38 | 1.37k | BigInt sqrt_modulo_prime(const BigInt& a, const BigInt& p) { |
39 | 1.37k | BOTAN_ARG_CHECK(p > 1, "invalid prime"); |
40 | 1.37k | BOTAN_ARG_CHECK(a < p, "value to solve for must be less than p"); |
41 | 1.37k | BOTAN_ARG_CHECK(a >= 0, "value to solve for must not be negative"); |
42 | | |
43 | | // some very easy cases |
44 | 1.37k | if(p == 2 || a <= 1) { |
45 | 0 | return a; |
46 | 0 | } |
47 | | |
48 | 1.37k | BOTAN_ARG_CHECK(p.is_odd(), "invalid prime"); |
49 | | |
50 | 1.37k | if(jacobi(a, p) != 1) { // not a quadratic residue |
51 | 522 | return BigInt::from_s32(-1); |
52 | 522 | } |
53 | | |
54 | 853 | Modular_Reducer mod_p(p); |
55 | 853 | auto monty_p = std::make_shared<Montgomery_Params>(p, mod_p); |
56 | | |
57 | 853 | if(p % 4 == 3) // The easy case |
58 | 853 | { |
59 | 853 | return monty_exp_vartime(monty_p, a, ((p + 1) >> 2)); |
60 | 853 | } |
61 | | |
62 | | // Otherwise we have to use Shanks-Tonelli |
63 | 0 | size_t s = low_zero_bits(p - 1); |
64 | 0 | BigInt q = p >> s; |
65 | |
|
66 | 0 | q -= 1; |
67 | 0 | q >>= 1; |
68 | |
|
69 | 0 | BigInt r = monty_exp_vartime(monty_p, a, q); |
70 | 0 | BigInt n = mod_p.multiply(a, mod_p.square(r)); |
71 | 0 | r = mod_p.multiply(r, a); |
72 | |
|
73 | 0 | if(n == 1) { |
74 | 0 | return r; |
75 | 0 | } |
76 | | |
77 | | // find random quadratic nonresidue z |
78 | 0 | word z = 2; |
79 | 0 | for(;;) { |
80 | 0 | if(jacobi(BigInt::from_word(z), p) == -1) { // found one |
81 | 0 | break; |
82 | 0 | } |
83 | | |
84 | 0 | z += 1; // try next z |
85 | | |
86 | | /* |
87 | | * The expected number of tests to find a non-residue modulo a |
88 | | * prime is 2. If we have not found one after 256 then almost |
89 | | * certainly we have been given a non-prime p. |
90 | | */ |
91 | 0 | if(z >= 256) { |
92 | 0 | return BigInt::from_s32(-1); |
93 | 0 | } |
94 | 0 | } |
95 | | |
96 | 0 | BigInt c = monty_exp_vartime(monty_p, BigInt::from_word(z), (q << 1) + 1); |
97 | |
|
98 | 0 | while(n > 1) { |
99 | 0 | q = n; |
100 | |
|
101 | 0 | size_t i = 0; |
102 | 0 | while(q != 1) { |
103 | 0 | q = mod_p.square(q); |
104 | 0 | ++i; |
105 | |
|
106 | 0 | if(i >= s) { |
107 | 0 | return BigInt::from_s32(-1); |
108 | 0 | } |
109 | 0 | } |
110 | | |
111 | 0 | BOTAN_ASSERT_NOMSG(s >= (i + 1)); // No underflow! |
112 | 0 | c = monty_exp_vartime(monty_p, c, BigInt::power_of_2(s - i - 1)); |
113 | 0 | r = mod_p.multiply(r, c); |
114 | 0 | c = mod_p.square(c); |
115 | 0 | n = mod_p.multiply(n, c); |
116 | | |
117 | | // s decreases as the algorithm proceeds |
118 | 0 | BOTAN_ASSERT_NOMSG(s >= i); |
119 | 0 | s = i; |
120 | 0 | } |
121 | | |
122 | 0 | return r; |
123 | 0 | } |
124 | | |
125 | | /* |
126 | | * Calculate the Jacobi symbol |
127 | | */ |
128 | 1.37k | int32_t jacobi(const BigInt& a, const BigInt& n) { |
129 | 1.37k | if(n.is_even() || n < 2) { |
130 | 0 | throw Invalid_Argument("jacobi: second argument must be odd and > 1"); |
131 | 0 | } |
132 | | |
133 | 1.37k | BigInt x = a % n; |
134 | 1.37k | BigInt y = n; |
135 | 1.37k | int32_t J = 1; |
136 | | |
137 | 153k | while(y > 1) { |
138 | 152k | x %= y; |
139 | 152k | if(x > y / 2) { |
140 | 72.1k | x = y - x; |
141 | 72.1k | if(y % 4 == 3) { |
142 | 36.4k | J = -J; |
143 | 36.4k | } |
144 | 72.1k | } |
145 | 152k | if(x.is_zero()) { |
146 | 0 | return 0; |
147 | 0 | } |
148 | | |
149 | 152k | size_t shifts = low_zero_bits(x); |
150 | 152k | x >>= shifts; |
151 | 152k | if(shifts % 2) { |
152 | 49.7k | word y_mod_8 = y % 8; |
153 | 49.7k | if(y_mod_8 == 3 || y_mod_8 == 5) { |
154 | 23.7k | J = -J; |
155 | 23.7k | } |
156 | 49.7k | } |
157 | | |
158 | 152k | if(x % 4 == 3 && y % 4 == 3) { |
159 | 36.8k | J = -J; |
160 | 36.8k | } |
161 | 152k | std::swap(x, y); |
162 | 152k | } |
163 | 1.37k | return J; |
164 | 1.37k | } |
165 | | |
166 | | /* |
167 | | * Square a BigInt |
168 | | */ |
169 | 853 | BigInt square(const BigInt& x) { |
170 | 853 | BigInt z = x; |
171 | 853 | secure_vector<word> ws; |
172 | 853 | z.square(ws); |
173 | 853 | return z; |
174 | 853 | } |
175 | | |
176 | | /* |
177 | | * Return the number of 0 bits at the end of n |
178 | | */ |
179 | 152k | size_t low_zero_bits(const BigInt& n) { |
180 | 152k | size_t low_zero = 0; |
181 | | |
182 | 152k | auto seen_nonempty_word = CT::Mask<word>::cleared(); |
183 | | |
184 | 1.37M | for(size_t i = 0; i != n.size(); ++i) { |
185 | 1.22M | const word x = n.word_at(i); |
186 | | |
187 | | // ctz(0) will return sizeof(word) |
188 | 1.22M | const size_t tz_x = ctz(x); |
189 | | |
190 | | // if x > 0 we want to count tz_x in total but not any |
191 | | // further words, so set the mask after the addition |
192 | 1.22M | low_zero += seen_nonempty_word.if_not_set_return(tz_x); |
193 | | |
194 | 1.22M | seen_nonempty_word |= CT::Mask<word>::expand(x); |
195 | 1.22M | } |
196 | | |
197 | | // if we saw no words with x > 0 then n == 0 and the value we have |
198 | | // computed is meaningless. Instead return BigInt::zero() in that case. |
199 | 152k | return seen_nonempty_word.if_set_return(low_zero); |
200 | 152k | } |
201 | | |
202 | | namespace { |
203 | | |
204 | 0 | size_t safegcd_loop_bound(size_t f_bits, size_t g_bits) { |
205 | 0 | const size_t d = std::max(f_bits, g_bits); |
206 | 0 | return 4 + 3 * d; |
207 | 0 | } |
208 | | |
209 | | } // namespace |
210 | | |
211 | | /* |
212 | | * Calculate the GCD |
213 | | */ |
214 | 0 | BigInt gcd(const BigInt& a, const BigInt& b) { |
215 | 0 | if(a.is_zero()) { |
216 | 0 | return abs(b); |
217 | 0 | } |
218 | 0 | if(b.is_zero()) { |
219 | 0 | return abs(a); |
220 | 0 | } |
221 | 0 | if(a == 1 || b == 1) { |
222 | 0 | return BigInt::one(); |
223 | 0 | } |
224 | | |
225 | | // See https://gcd.cr.yp.to/safegcd-20190413.pdf fig 1.2 |
226 | | |
227 | 0 | BigInt f = a; |
228 | 0 | BigInt g = b; |
229 | 0 | f.const_time_poison(); |
230 | 0 | g.const_time_poison(); |
231 | |
|
232 | 0 | f.set_sign(BigInt::Positive); |
233 | 0 | g.set_sign(BigInt::Positive); |
234 | |
|
235 | 0 | const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g)); |
236 | 0 | CT::unpoison(common2s); |
237 | |
|
238 | 0 | f >>= common2s; |
239 | 0 | g >>= common2s; |
240 | |
|
241 | 0 | f.ct_cond_swap(f.is_even(), g); |
242 | |
|
243 | 0 | int32_t delta = 1; |
244 | |
|
245 | 0 | const size_t loop_cnt = safegcd_loop_bound(f.bits(), g.bits()); |
246 | |
|
247 | 0 | BigInt newg, t; |
248 | 0 | for(size_t i = 0; i != loop_cnt; ++i) { |
249 | 0 | sub_abs(newg, f, g); |
250 | |
|
251 | 0 | const bool need_swap = (g.is_odd() && delta > 0); |
252 | | |
253 | | // if(need_swap) { delta *= -1 } else { delta *= 1 } |
254 | 0 | delta *= CT::Mask<uint8_t>::expand(need_swap).if_not_set_return(2) - 1; |
255 | 0 | f.ct_cond_swap(need_swap, g); |
256 | 0 | g.ct_cond_swap(need_swap, newg); |
257 | |
|
258 | 0 | delta += 1; |
259 | |
|
260 | 0 | g.ct_cond_add(g.is_odd(), f); |
261 | 0 | g >>= 1; |
262 | 0 | } |
263 | |
|
264 | 0 | f <<= common2s; |
265 | |
|
266 | 0 | f.const_time_unpoison(); |
267 | 0 | g.const_time_unpoison(); |
268 | |
|
269 | 0 | BOTAN_ASSERT_NOMSG(g.is_zero()); |
270 | |
|
271 | 0 | return f; |
272 | 0 | } |
273 | | |
274 | | /* |
275 | | * Calculate the LCM |
276 | | */ |
277 | 0 | BigInt lcm(const BigInt& a, const BigInt& b) { |
278 | 0 | if(a == b) { |
279 | 0 | return a; |
280 | 0 | } |
281 | | |
282 | 0 | auto ab = a * b; |
283 | 0 | ab.set_sign(BigInt::Positive); // ignore the signs of a & b |
284 | 0 | const auto g = gcd(a, b); |
285 | 0 | return ct_divide(ab, g); |
286 | 0 | } |
287 | | |
288 | | /* |
289 | | * Modular Exponentiation |
290 | | */ |
291 | 0 | BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) { |
292 | 0 | if(mod.is_negative() || mod == 1) { |
293 | 0 | return BigInt::zero(); |
294 | 0 | } |
295 | | |
296 | 0 | if(base.is_zero() || mod.is_zero()) { |
297 | 0 | if(exp.is_zero()) { |
298 | 0 | return BigInt::one(); |
299 | 0 | } |
300 | 0 | return BigInt::zero(); |
301 | 0 | } |
302 | | |
303 | 0 | Modular_Reducer reduce_mod(mod); |
304 | |
|
305 | 0 | const size_t exp_bits = exp.bits(); |
306 | |
|
307 | 0 | if(mod.is_odd()) { |
308 | 0 | auto monty_params = std::make_shared<Montgomery_Params>(mod, reduce_mod); |
309 | 0 | return monty_exp(monty_params, reduce_mod.reduce(base), exp, exp_bits); |
310 | 0 | } |
311 | | |
312 | | /* |
313 | | Support for even modulus is just a convenience and not considered |
314 | | cryptographically important, so this implementation is slow ... |
315 | | */ |
316 | 0 | BigInt accum = BigInt::one(); |
317 | 0 | BigInt g = reduce_mod.reduce(base); |
318 | 0 | BigInt t; |
319 | |
|
320 | 0 | for(size_t i = 0; i != exp_bits; ++i) { |
321 | 0 | t = reduce_mod.multiply(g, accum); |
322 | 0 | g = reduce_mod.square(g); |
323 | 0 | accum.ct_cond_assign(exp.get_bit(i), t); |
324 | 0 | } |
325 | 0 | return accum; |
326 | 0 | } |
327 | | |
328 | 0 | BigInt is_perfect_square(const BigInt& C) { |
329 | 0 | if(C < 1) { |
330 | 0 | throw Invalid_Argument("is_perfect_square requires C >= 1"); |
331 | 0 | } |
332 | 0 | if(C == 1) { |
333 | 0 | return BigInt::one(); |
334 | 0 | } |
335 | | |
336 | 0 | const size_t n = C.bits(); |
337 | 0 | const size_t m = (n + 1) / 2; |
338 | 0 | const BigInt B = C + BigInt::power_of_2(m); |
339 | |
|
340 | 0 | BigInt X = BigInt::power_of_2(m) - 1; |
341 | 0 | BigInt X2 = (X * X); |
342 | |
|
343 | 0 | for(;;) { |
344 | 0 | X = (X2 + C) / (2 * X); |
345 | 0 | X2 = (X * X); |
346 | |
|
347 | 0 | if(X2 < B) { |
348 | 0 | break; |
349 | 0 | } |
350 | 0 | } |
351 | |
|
352 | 0 | if(X2 == C) { |
353 | 0 | return X; |
354 | 0 | } else { |
355 | 0 | return BigInt::zero(); |
356 | 0 | } |
357 | 0 | } |
358 | | |
359 | | /* |
360 | | * Test for primality using Miller-Rabin |
361 | | */ |
362 | 0 | bool is_prime(const BigInt& n, RandomNumberGenerator& rng, size_t prob, bool is_random) { |
363 | 0 | if(n == 2) { |
364 | 0 | return true; |
365 | 0 | } |
366 | 0 | if(n <= 1 || n.is_even()) { |
367 | 0 | return false; |
368 | 0 | } |
369 | | |
370 | 0 | const size_t n_bits = n.bits(); |
371 | | |
372 | | // Fast path testing for small numbers (<= 65521) |
373 | 0 | if(n_bits <= 16) { |
374 | 0 | const uint16_t num = static_cast<uint16_t>(n.word_at(0)); |
375 | |
|
376 | 0 | return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num); |
377 | 0 | } |
378 | | |
379 | 0 | Modular_Reducer mod_n(n); |
380 | |
|
381 | 0 | if(rng.is_seeded()) { |
382 | 0 | const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random); |
383 | |
|
384 | 0 | if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false) { |
385 | 0 | return false; |
386 | 0 | } |
387 | | |
388 | 0 | if(is_random) { |
389 | 0 | return true; |
390 | 0 | } else { |
391 | 0 | return is_lucas_probable_prime(n, mod_n); |
392 | 0 | } |
393 | 0 | } else { |
394 | 0 | return is_bailie_psw_probable_prime(n, mod_n); |
395 | 0 | } |
396 | 0 | } |
397 | | |
398 | | } // namespace Botan |