/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 | 646k | { |
26 | 646k | const size_t x_sw = x.sig_words(); |
27 | 646k | const size_t y_sw = y.sig_words(); |
28 | 646k | z.resize(std::max(x_sw, y_sw)); |
29 | | |
30 | 646k | bigint_sub_abs(z.mutable_data(), |
31 | 646k | x.data(), x_sw, |
32 | 646k | y.data(), y_sw); |
33 | 646k | } |
34 | | |
35 | | } |
36 | | |
37 | | /* |
38 | | * Tonelli-Shanks algorithm |
39 | | */ |
40 | | BigInt sqrt_modulo_prime(const BigInt& a, const BigInt& p) |
41 | 33.7k | { |
42 | 33.7k | BOTAN_ARG_CHECK(p > 1, "invalid prime"); |
43 | 33.7k | BOTAN_ARG_CHECK(a < p, "value to solve for must be less than p"); |
44 | 33.7k | BOTAN_ARG_CHECK(a >= 0, "value to solve for must not be negative"); |
45 | | |
46 | | // some very easy cases |
47 | 33.7k | if(p == 2 || a <= 1) |
48 | 14 | return a; |
49 | | |
50 | 33.7k | BOTAN_ARG_CHECK(p.is_odd(), "invalid prime"); |
51 | | |
52 | 33.7k | if(jacobi(a, p) != 1) // not a quadratic residue |
53 | 12.5k | return BigInt::from_s32(-1); |
54 | | |
55 | 21.1k | Modular_Reducer mod_p(p); |
56 | 21.1k | auto monty_p = std::make_shared<Montgomery_Params>(p, mod_p); |
57 | | |
58 | 21.1k | if(p % 4 == 3) // The easy case |
59 | 18.9k | { |
60 | 18.9k | return monty_exp_vartime(monty_p, a, ((p+1) >> 2)); |
61 | 18.9k | } |
62 | | |
63 | | // Otherwise we have to use Shanks-Tonelli |
64 | 2.20k | size_t s = low_zero_bits(p - 1); |
65 | 2.20k | BigInt q = p >> s; |
66 | | |
67 | 2.20k | q -= 1; |
68 | 2.20k | q >>= 1; |
69 | | |
70 | 2.20k | BigInt r = monty_exp_vartime(monty_p, a, q); |
71 | 2.20k | BigInt n = mod_p.multiply(a, mod_p.square(r)); |
72 | 2.20k | r = mod_p.multiply(r, a); |
73 | | |
74 | 2.20k | if(n == 1) |
75 | 139 | return r; |
76 | | |
77 | | // find random quadratic nonresidue z |
78 | 2.06k | word z = 2; |
79 | 2.06k | for(;;) |
80 | 19.3k | { |
81 | 19.3k | if(jacobi(BigInt::from_word(z), p) == -1) // found one |
82 | 2.06k | break; |
83 | | |
84 | 17.2k | 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 | 17.2k | if(z >= 256) |
92 | 0 | return BigInt::from_s32(-1); |
93 | 17.2k | } |
94 | | |
95 | 2.06k | BigInt c = monty_exp_vartime(monty_p, BigInt::from_word(z), (q << 1) + 1); |
96 | | |
97 | 75.8k | while(n > 1) |
98 | 73.7k | { |
99 | 73.7k | q = n; |
100 | | |
101 | 73.7k | size_t i = 0; |
102 | 3.79M | while(q != 1) |
103 | 3.71M | { |
104 | 3.71M | q = mod_p.square(q); |
105 | 3.71M | ++i; |
106 | | |
107 | 3.71M | if(i >= s) |
108 | 0 | { |
109 | 0 | return BigInt::from_s32(-1); |
110 | 0 | } |
111 | 3.71M | } |
112 | | |
113 | 73.7k | BOTAN_ASSERT_NOMSG(s >= (i + 1)); // No underflow! |
114 | 73.7k | c = monty_exp_vartime(monty_p, c, BigInt::power_of_2(s-i-1)); |
115 | 73.7k | r = mod_p.multiply(r, c); |
116 | 73.7k | c = mod_p.square(c); |
117 | 73.7k | n = mod_p.multiply(n, c); |
118 | | |
119 | | // s decreases as the algorithm proceeds |
120 | 73.7k | BOTAN_ASSERT_NOMSG(s >= i); |
121 | 73.7k | s = i; |
122 | 73.7k | } |
123 | | |
124 | 2.06k | return r; |
125 | 2.06k | } |
126 | | |
127 | | /* |
128 | | * Calculate the Jacobi symbol |
129 | | */ |
130 | | int32_t jacobi(const BigInt& a, const BigInt& n) |
131 | 61.2k | { |
132 | 61.2k | if(n.is_even() || n < 2) |
133 | 0 | throw Invalid_Argument("jacobi: second argument must be odd and > 1"); |
134 | | |
135 | 61.2k | BigInt x = a % n; |
136 | 61.2k | BigInt y = n; |
137 | 61.2k | int32_t J = 1; |
138 | | |
139 | 3.63M | while(y > 1) |
140 | 3.57M | { |
141 | 3.57M | x %= y; |
142 | 3.57M | if(x > y / 2) |
143 | 1.65M | { |
144 | 1.65M | x = y - x; |
145 | 1.65M | if(y % 4 == 3) |
146 | 833k | J = -J; |
147 | 1.65M | } |
148 | 3.57M | if(x.is_zero()) |
149 | 0 | return 0; |
150 | | |
151 | 3.57M | size_t shifts = low_zero_bits(x); |
152 | 3.57M | x >>= shifts; |
153 | 3.57M | if(shifts % 2) |
154 | 1.15M | { |
155 | 1.15M | word y_mod_8 = y % 8; |
156 | 1.15M | if(y_mod_8 == 3 || y_mod_8 == 5) |
157 | 566k | J = -J; |
158 | 1.15M | } |
159 | | |
160 | 3.57M | if(x % 4 == 3 && y % 4 == 3) |
161 | 877k | J = -J; |
162 | 3.57M | std::swap(x, y); |
163 | 3.57M | } |
164 | 61.2k | return J; |
165 | 61.2k | } |
166 | | |
167 | | /* |
168 | | * Square a BigInt |
169 | | */ |
170 | | BigInt square(const BigInt& x) |
171 | 6.31M | { |
172 | 6.31M | BigInt z = x; |
173 | 6.31M | secure_vector<word> ws; |
174 | 6.31M | z.square(ws); |
175 | 6.31M | return z; |
176 | 6.31M | } |
177 | | |
178 | | /* |
179 | | * Return the number of 0 bits at the end of n |
180 | | */ |
181 | | size_t low_zero_bits(const BigInt& n) |
182 | 4.65M | { |
183 | 4.65M | size_t low_zero = 0; |
184 | | |
185 | 4.65M | auto seen_nonempty_word = CT::Mask<word>::cleared(); |
186 | | |
187 | 77.6M | for(size_t i = 0; i != n.size(); ++i) |
188 | 73.0M | { |
189 | 73.0M | const word x = n.word_at(i); |
190 | | |
191 | | // ctz(0) will return sizeof(word) |
192 | 73.0M | const size_t tz_x = ctz(x); |
193 | | |
194 | | // if x > 0 we want to count tz_x in total but not any |
195 | | // further words, so set the mask after the addition |
196 | 73.0M | low_zero += seen_nonempty_word.if_not_set_return(tz_x); |
197 | | |
198 | 73.0M | seen_nonempty_word |= CT::Mask<word>::expand(x); |
199 | 73.0M | } |
200 | | |
201 | | // if we saw no words with x > 0 then n == 0 and the value we have |
202 | | // computed is meaningless. Instead return BigInt::zero() in that case. |
203 | 4.65M | return seen_nonempty_word.if_set_return(low_zero); |
204 | 4.65M | } |
205 | | |
206 | | |
207 | | namespace { |
208 | | |
209 | | size_t safegcd_loop_bound(size_t f_bits, size_t g_bits) |
210 | 1.54k | { |
211 | 1.54k | const size_t d = std::max(f_bits, g_bits); |
212 | 1.54k | return 4 + 3*d; |
213 | 1.54k | } |
214 | | |
215 | | } |
216 | | |
217 | | /* |
218 | | * Calculate the GCD |
219 | | */ |
220 | | BigInt gcd(const BigInt& a, const BigInt& b) |
221 | 1.68k | { |
222 | 1.68k | if(a.is_zero()) |
223 | 131 | return abs(b); |
224 | 1.54k | if(b.is_zero()) |
225 | 6 | return abs(a); |
226 | 1.54k | if(a == 1 || b == 1) |
227 | 3 | return BigInt::one(); |
228 | | |
229 | | // See https://gcd.cr.yp.to/safegcd-20190413.pdf fig 1.2 |
230 | | |
231 | 1.54k | BigInt f = a; |
232 | 1.54k | BigInt g = b; |
233 | 1.54k | f.const_time_poison(); |
234 | 1.54k | g.const_time_poison(); |
235 | | |
236 | 1.54k | f.set_sign(BigInt::Positive); |
237 | 1.54k | g.set_sign(BigInt::Positive); |
238 | | |
239 | 1.54k | const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g)); |
240 | 1.54k | CT::unpoison(common2s); |
241 | | |
242 | 1.54k | f >>= common2s; |
243 | 1.54k | g >>= common2s; |
244 | | |
245 | 1.54k | f.ct_cond_swap(f.is_even(), g); |
246 | | |
247 | 1.54k | int32_t delta = 1; |
248 | | |
249 | 1.54k | const size_t loop_cnt = safegcd_loop_bound(f.bits(), g.bits()); |
250 | | |
251 | 1.54k | BigInt newg, t; |
252 | 647k | for(size_t i = 0; i != loop_cnt; ++i) |
253 | 646k | { |
254 | 646k | sub_abs(newg, f, g); |
255 | | |
256 | 646k | const bool need_swap = (g.is_odd() && delta > 0); |
257 | | |
258 | | // if(need_swap) { delta *= -1 } else { delta *= 1 } |
259 | 646k | delta *= CT::Mask<uint8_t>::expand(need_swap).if_not_set_return(2) - 1; |
260 | 646k | f.ct_cond_swap(need_swap, g); |
261 | 646k | g.ct_cond_swap(need_swap, newg); |
262 | | |
263 | 646k | delta += 1; |
264 | | |
265 | 646k | g.ct_cond_add(g.is_odd(), f); |
266 | 646k | g >>= 1; |
267 | 646k | } |
268 | | |
269 | 1.54k | f <<= common2s; |
270 | | |
271 | 1.54k | f.const_time_unpoison(); |
272 | 1.54k | g.const_time_unpoison(); |
273 | | |
274 | 1.54k | BOTAN_ASSERT_NOMSG(g.is_zero()); |
275 | | |
276 | 1.54k | return f; |
277 | 1.54k | } |
278 | | |
279 | | /* |
280 | | * Calculate the LCM |
281 | | */ |
282 | | BigInt lcm(const BigInt& a, const BigInt& b) |
283 | 0 | { |
284 | 0 | if(a == b) |
285 | 0 | return a; |
286 | | |
287 | 0 | auto ab = a * b; |
288 | 0 | ab.set_sign(BigInt::Positive); // ignore the signs of a & b |
289 | 0 | const auto g = gcd(a, b); |
290 | 0 | return ct_divide(ab, g); |
291 | 0 | } |
292 | | |
293 | | /* |
294 | | * Modular Exponentiation |
295 | | */ |
296 | | BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) |
297 | 1.54k | { |
298 | 1.54k | if(mod.is_negative() || mod == 1) |
299 | 12 | { |
300 | 12 | return BigInt::zero(); |
301 | 12 | } |
302 | | |
303 | 1.53k | if(base.is_zero() || mod.is_zero()) |
304 | 26 | { |
305 | 26 | if(exp.is_zero()) |
306 | 1 | return BigInt::one(); |
307 | 25 | return BigInt::zero(); |
308 | 26 | } |
309 | | |
310 | 1.50k | Modular_Reducer reduce_mod(mod); |
311 | | |
312 | 1.50k | const size_t exp_bits = exp.bits(); |
313 | | |
314 | 1.50k | if(mod.is_odd()) |
315 | 782 | { |
316 | 782 | auto monty_params = std::make_shared<Montgomery_Params>(mod, reduce_mod); |
317 | 782 | return monty_exp(monty_params, reduce_mod.reduce(base), exp, exp_bits); |
318 | 782 | } |
319 | | |
320 | | /* |
321 | | Support for even modulus is just a convenience and not considered |
322 | | cryptographically important, so this implementation is slow ... |
323 | | */ |
324 | 726 | BigInt accum = BigInt::one(); |
325 | 726 | BigInt g = reduce_mod.reduce(base); |
326 | 726 | BigInt t; |
327 | | |
328 | 335k | for(size_t i = 0; i != exp_bits; ++i) |
329 | 334k | { |
330 | 334k | t = reduce_mod.multiply(g, accum); |
331 | 334k | g = reduce_mod.square(g); |
332 | 334k | accum.ct_cond_assign(exp.get_bit(i), t); |
333 | 334k | } |
334 | 726 | return accum; |
335 | 1.50k | } |
336 | | |
337 | | |
338 | | BigInt is_perfect_square(const BigInt& C) |
339 | 315 | { |
340 | 315 | if(C < 1) |
341 | 0 | throw Invalid_Argument("is_perfect_square requires C >= 1"); |
342 | 315 | if(C == 1) |
343 | 0 | return BigInt::one(); |
344 | | |
345 | 315 | const size_t n = C.bits(); |
346 | 315 | const size_t m = (n + 1) / 2; |
347 | 315 | const BigInt B = C + BigInt::power_of_2(m); |
348 | | |
349 | 315 | BigInt X = BigInt::power_of_2(m) - 1; |
350 | 315 | BigInt X2 = (X*X); |
351 | | |
352 | 315 | for(;;) |
353 | 1.63k | { |
354 | 1.63k | X = (X2 + C) / (2*X); |
355 | 1.63k | X2 = (X*X); |
356 | | |
357 | 1.63k | if(X2 < B) |
358 | 315 | break; |
359 | 1.63k | } |
360 | | |
361 | 315 | if(X2 == C) |
362 | 0 | return X; |
363 | 315 | else |
364 | 315 | return BigInt::zero(); |
365 | 315 | } |
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 | | } |