Botan 3.0.0-alpha0
Crypto and TLS for C&
numthry.cpp
Go to the documentation of this file.
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
20namespace Botan {
21
22namespace {
23
24void sub_abs(BigInt& z, const BigInt& x, const BigInt& y)
25 {
26 const size_t x_sw = x.sig_words();
27 const size_t y_sw = y.sig_words();
28 z.resize(std::max(x_sw, y_sw));
29
30 bigint_sub_abs(z.mutable_data(),
31 x.data(), x_sw,
32 y.data(), y_sw);
33 }
34
35}
36
37/*
38* Tonelli-Shanks algorithm
39*/
41 {
42 BOTAN_ARG_CHECK(p > 1, "invalid prime");
43 BOTAN_ARG_CHECK(a < p, "value to solve for must be less than p");
44 BOTAN_ARG_CHECK(a >= 0, "value to solve for must not be negative");
45
46 // some very easy cases
47 if(p == 2 || a <= 1)
48 return a;
49
50 BOTAN_ARG_CHECK(p.is_odd(), "invalid prime");
51
52 if(jacobi(a, p) != 1) // not a quadratic residue
53 return BigInt::from_s32(-1);
54
55 Modular_Reducer mod_p(p);
56 auto monty_p = std::make_shared<Montgomery_Params>(p, mod_p);
57
58 if(p % 4 == 3) // The easy case
59 {
60 return monty_exp_vartime(monty_p, a, ((p+1) >> 2));
61 }
62
63 // Otherwise we have to use Shanks-Tonelli
64 size_t s = low_zero_bits(p - 1);
65 BigInt q = p >> s;
66
67 q -= 1;
68 q >>= 1;
69
70 BigInt r = monty_exp_vartime(monty_p, a, q);
71 BigInt n = mod_p.multiply(a, mod_p.square(r));
72 r = mod_p.multiply(r, a);
73
74 if(n == 1)
75 return r;
76
77 // find random quadratic nonresidue z
78 word z = 2;
79 for(;;)
80 {
81 if(jacobi(BigInt::from_word(z), p) == -1) // found one
82 break;
83
84 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 if(z >= 256)
92 return BigInt::from_s32(-1);
93 }
94
95 BigInt c = monty_exp_vartime(monty_p, BigInt::from_word(z), (q << 1) + 1);
96
97 while(n > 1)
98 {
99 q = n;
100
101 size_t i = 0;
102 while(q != 1)
103 {
104 q = mod_p.square(q);
105 ++i;
106
107 if(i >= s)
108 {
109 return BigInt::from_s32(-1);
110 }
111 }
112
113 BOTAN_ASSERT_NOMSG(s >= (i + 1)); // No underflow!
114 c = monty_exp_vartime(monty_p, c, BigInt::power_of_2(s-i-1));
115 r = mod_p.multiply(r, c);
116 c = mod_p.square(c);
117 n = mod_p.multiply(n, c);
118
119 // s decreases as the algorithm proceeds
120 BOTAN_ASSERT_NOMSG(s >= i);
121 s = i;
122 }
123
124 return r;
125 }
126
127/*
128* Calculate the Jacobi symbol
129*/
130int32_t jacobi(const BigInt& a, const BigInt& n)
131 {
132 if(n.is_even() || n < 2)
133 throw Invalid_Argument("jacobi: second argument must be odd and > 1");
134
135 BigInt x = a % n;
136 BigInt y = n;
137 int32_t J = 1;
138
139 while(y > 1)
140 {
141 x %= y;
142 if(x > y / 2)
143 {
144 x = y - x;
145 if(y % 4 == 3)
146 J = -J;
147 }
148 if(x.is_zero())
149 return 0;
150
151 size_t shifts = low_zero_bits(x);
152 x >>= shifts;
153 if(shifts % 2)
154 {
155 word y_mod_8 = y % 8;
156 if(y_mod_8 == 3 || y_mod_8 == 5)
157 J = -J;
158 }
159
160 if(x % 4 == 3 && y % 4 == 3)
161 J = -J;
162 std::swap(x, y);
163 }
164 return J;
165 }
166
167/*
168* Square a BigInt
169*/
171 {
172 BigInt z = x;
174 z.square(ws);
175 return z;
176 }
177
178/*
179* Return the number of 0 bits at the end of n
180*/
181size_t low_zero_bits(const BigInt& n)
182 {
183 size_t low_zero = 0;
184
185 auto seen_nonempty_word = CT::Mask<word>::cleared();
186
187 for(size_t i = 0; i != n.size(); ++i)
188 {
189 const word x = n.word_at(i);
190
191 // ctz(0) will return sizeof(word)
192 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 low_zero += seen_nonempty_word.if_not_set_return(tz_x);
197
198 seen_nonempty_word |= CT::Mask<word>::expand(x);
199 }
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 return seen_nonempty_word.if_set_return(low_zero);
204 }
205
206
207namespace {
208
209size_t safegcd_loop_bound(size_t f_bits, size_t g_bits)
210 {
211 const size_t d = std::max(f_bits, g_bits);
212 return 4 + 3*d;
213 }
214
215}
216
217/*
218* Calculate the GCD
219*/
220BigInt gcd(const BigInt& a, const BigInt& b)
221 {
222 if(a.is_zero())
223 return abs(b);
224 if(b.is_zero())
225 return abs(a);
226 if(a == 1 || b == 1)
227 return BigInt::one();
228
229 // See https://gcd.cr.yp.to/safegcd-20190413.pdf fig 1.2
230
231 BigInt f = a;
232 BigInt g = b;
235
238
239 const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g));
240 CT::unpoison(common2s);
241
242 f >>= common2s;
243 g >>= common2s;
244
245 f.ct_cond_swap(f.is_even(), g);
246
247 int32_t delta = 1;
248
249 const size_t loop_cnt = safegcd_loop_bound(f.bits(), g.bits());
250
251 BigInt newg, t;
252 for(size_t i = 0; i != loop_cnt; ++i)
253 {
254 sub_abs(newg, f, g);
255
256 const bool need_swap = (g.is_odd() && delta > 0);
257
258 // if(need_swap) { delta *= -1 } else { delta *= 1 }
259 delta *= CT::Mask<uint8_t>::expand(need_swap).if_not_set_return(2) - 1;
260 f.ct_cond_swap(need_swap, g);
261 g.ct_cond_swap(need_swap, newg);
262
263 delta += 1;
264
265 g.ct_cond_add(g.is_odd(), f);
266 g >>= 1;
267 }
268
269 f <<= common2s;
270
273
275
276 return f;
277 }
278
279/*
280* Calculate the LCM
281*/
282BigInt lcm(const BigInt& a, const BigInt& b)
283 {
284 if(a == b)
285 return a;
286
287 auto ab = a * b;
288 ab.set_sign(BigInt::Positive); // ignore the signs of a & b
289 const auto g = gcd(a, b);
290 return ct_divide(ab, g);
291 }
292
293/*
294* Modular Exponentiation
295*/
296BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
297 {
298 if(mod.is_negative() || mod == 1)
299 {
300 return BigInt::zero();
301 }
302
303 if(base.is_zero() || mod.is_zero())
304 {
305 if(exp.is_zero())
306 return BigInt::one();
307 return BigInt::zero();
308 }
309
310 Modular_Reducer reduce_mod(mod);
311
312 const size_t exp_bits = exp.bits();
313
314 if(mod.is_odd())
315 {
316 auto monty_params = std::make_shared<Montgomery_Params>(mod, reduce_mod);
317 return monty_exp(monty_params, reduce_mod.reduce(base), exp, exp_bits);
318 }
319
320 /*
321 Support for even modulus is just a convenience and not considered
322 cryptographically important, so this implementation is slow ...
323 */
324 BigInt accum = BigInt::one();
325 BigInt g = reduce_mod.reduce(base);
326 BigInt t;
327
328 for(size_t i = 0; i != exp_bits; ++i)
329 {
330 t = reduce_mod.multiply(g, accum);
331 g = reduce_mod.square(g);
332 accum.ct_cond_assign(exp.get_bit(i), t);
333 }
334 return accum;
335 }
336
337
339 {
340 if(C < 1)
341 throw Invalid_Argument("is_perfect_square requires C >= 1");
342 if(C == 1)
343 return BigInt::one();
344
345 const size_t n = C.bits();
346 const size_t m = (n + 1) / 2;
347 const BigInt B = C + BigInt::power_of_2(m);
348
349 BigInt X = BigInt::power_of_2(m) - 1;
350 BigInt X2 = (X*X);
351
352 for(;;)
353 {
354 X = (X2 + C) / (2*X);
355 X2 = (X*X);
356
357 if(X2 < B)
358 break;
359 }
360
361 if(X2 == C)
362 return X;
363 else
364 return BigInt::zero();
365 }
366
367/*
368* Test for primality using Miller-Rabin
369*/
370bool is_prime(const BigInt& n,
372 size_t prob,
373 bool is_random)
374 {
375 if(n == 2)
376 return true;
377 if(n <= 1 || n.is_even())
378 return false;
379
380 const size_t n_bits = n.bits();
381
382 // Fast path testing for small numbers (<= 65521)
383 if(n_bits <= 16)
384 {
385 const uint16_t num = static_cast<uint16_t>(n.word_at(0));
386
387 return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
388 }
389
390 Modular_Reducer mod_n(n);
391
392 if(rng.is_seeded())
393 {
394 const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
395
396 if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false)
397 return false;
398
399 if(is_random)
400 return true;
401 else
402 return is_lucas_probable_prime(n, mod_n);
403 }
404 else
405 {
406 return is_bailie_psw_probable_prime(n, mod_n);
407 }
408 }
409
410}
#define BOTAN_ASSERT_NOMSG(expr)
Definition: assert.h:67
#define BOTAN_ARG_CHECK(expr, msg)
Definition: assert.h:36
void ct_cond_add(bool predicate, const BigInt &value)
Definition: bigint.cpp:451
static BigInt zero()
Definition: bigint.h:45
bool is_odd() const
Definition: bigint.h:418
void ct_cond_assign(bool predicate, const BigInt &other)
Definition: bigint.cpp:484
size_t size() const
Definition: bigint.h:594
static BigInt one()
Definition: bigint.h:50
word word_at(size_t n) const
Definition: bigint.h:522
size_t bits() const
Definition: bigint.cpp:309
bool is_even() const
Definition: bigint.h:412
static BigInt power_of_2(size_t n)
Definition: bigint.h:753
void const_time_poison() const
Definition: bigint.h:734
static BigInt from_s32(int32_t n)
Definition: bigint.cpp:51
void ct_cond_swap(bool predicate, BigInt &other)
Definition: bigint.cpp:462
static BigInt from_word(word n)
Definition: bigint.cpp:43
bool is_zero() const
Definition: bigint.h:430
BigInt & square(secure_vector< word > &ws)
Definition: big_ops2.cpp:198
bool is_negative() const
Definition: bigint.h:541
void const_time_unpoison() const
Definition: bigint.h:735
bool get_bit(size_t n) const
Definition: bigint.h:479
void set_sign(Sign sign)
Definition: bigint.h:577
static Mask< T > expand(T v)
Definition: ct_utils.h:121
static Mask< T > cleared()
Definition: ct_utils.h:113
BigInt square(const BigInt &x) const
Definition: reducer.h:46
BigInt multiply(const BigInt &x, const BigInt &y) const
Definition: reducer.h:31
BigInt reduce(const BigInt &x) const
Definition: reducer.cpp:37
virtual bool is_seeded() const =0
fe X
Definition: ge.cpp:26
PolynomialVector b
Definition: kyber.cpp:821
void unpoison(const T *p, size_t n)
Definition: ct_utils.h:58
Definition: alg_id.cpp:13
BigInt power_mod(const BigInt &base, const BigInt &exp, const BigInt &mod)
Definition: numthry.cpp:296
constexpr size_t ctz(T n)
Definition: bit_ops.h:99
BigInt monty_exp(std::shared_ptr< const Montgomery_Params > params_p, const BigInt &g, const BigInt &k, size_t max_k_bits)
Definition: monty_exp.h:44
BigInt lcm(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:282
BigInt square(const BigInt &x)
Definition: numthry.cpp:170
const uint16_t PRIMES[]
Definition: primes.cpp:12
size_t low_zero_bits(const BigInt &n)
Definition: numthry.cpp:181
BigInt abs(const BigInt &n)
Definition: numthry.h:22
CT::Mask< word > bigint_sub_abs(word z[], const word x[], const word y[], size_t N, word ws[])
Definition: mp_core.h:377
const size_t PRIME_TABLE_SIZE
Definition: numthry.h:171
bool is_prime(const BigInt &n, RandomNumberGenerator &rng, size_t prob, bool is_random)
Definition: numthry.cpp:370
bool is_miller_rabin_probable_prime(const BigInt &n, const Modular_Reducer &mod_n, RandomNumberGenerator &rng, size_t test_iterations)
Definition: primality.cpp:147
BigInt monty_exp_vartime(std::shared_ptr< const Montgomery_Params > params_p, const BigInt &g, const BigInt &k)
Definition: monty_exp.h:52
void ct_divide(const BigInt &x, const BigInt &y, BigInt &q_out, BigInt &r_out)
Definition: divide.cpp:51
bool is_bailie_psw_probable_prime(const BigInt &n, const Modular_Reducer &mod_n)
Definition: primality.cpp:89
BigInt gcd(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:220
BigInt sqrt_modulo_prime(const BigInt &a, const BigInt &p)
Definition: numthry.cpp:40
BigInt is_perfect_square(const BigInt &C)
Definition: numthry.cpp:338
size_t miller_rabin_test_iterations(size_t n_bits, size_t prob, bool random)
Definition: primality.cpp:170
int32_t jacobi(const BigInt &a, const BigInt &n)
Definition: numthry.cpp:130
bool is_lucas_probable_prime(const BigInt &C, const Modular_Reducer &mod_C)
Definition: primality.cpp:17
std::vector< T, secure_allocator< T > > secure_vector
Definition: secmem.h:65