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