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