Botan 3.3.0
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
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
21namespace Botan {
22
23/*
24* Tonelli-Shanks algorithm
25*/
27 BOTAN_ARG_CHECK(p > 1, "invalid prime");
28 BOTAN_ARG_CHECK(a < p, "value to solve for must be less than p");
29 BOTAN_ARG_CHECK(a >= 0, "value to solve for must not be negative");
30
31 // some very easy cases
32 if(p == 2 || a <= 1) {
33 return a;
34 }
35
36 BOTAN_ARG_CHECK(p.is_odd(), "invalid prime");
37
38 if(jacobi(a, p) != 1) { // not a quadratic residue
39 return BigInt::from_s32(-1);
40 }
41
42 Modular_Reducer mod_p(p);
43 auto monty_p = std::make_shared<Montgomery_Params>(p, mod_p);
44
45 // If p == 3 (mod 4) there is a simple solution
46 if(p % 4 == 3) {
47 return monty_exp_vartime(monty_p, a, ((p + 1) >> 2));
48 }
49
50 // Otherwise we have to use Shanks-Tonelli
51 size_t s = low_zero_bits(p - 1);
52 BigInt q = p >> s;
53
54 q -= 1;
55 q >>= 1;
56
57 BigInt r = monty_exp_vartime(monty_p, a, q);
58 BigInt n = mod_p.multiply(a, mod_p.square(r));
59 r = mod_p.multiply(r, a);
60
61 if(n == 1) {
62 return r;
63 }
64
65 // find random quadratic nonresidue z
66 word z = 2;
67 for(;;) {
68 if(jacobi(BigInt::from_word(z), p) == -1) { // found one
69 break;
70 }
71
72 z += 1; // try next z
73
74 /*
75 * The expected number of tests to find a non-residue modulo a
76 * prime is 2. If we have not found one after 256 then almost
77 * certainly we have been given a non-prime p.
78 */
79 if(z >= 256) {
80 return BigInt::from_s32(-1);
81 }
82 }
83
84 BigInt c = monty_exp_vartime(monty_p, BigInt::from_word(z), (q << 1) + 1);
85
86 while(n > 1) {
87 q = n;
88
89 size_t i = 0;
90 while(q != 1) {
91 q = mod_p.square(q);
92 ++i;
93
94 if(i >= s) {
95 return BigInt::from_s32(-1);
96 }
97 }
98
99 BOTAN_ASSERT_NOMSG(s >= (i + 1)); // No underflow!
100 c = monty_exp_vartime(monty_p, c, BigInt::power_of_2(s - i - 1));
101 r = mod_p.multiply(r, c);
102 c = mod_p.square(c);
103 n = mod_p.multiply(n, c);
104
105 // s decreases as the algorithm proceeds
106 BOTAN_ASSERT_NOMSG(s >= i);
107 s = i;
108 }
109
110 return r;
111}
112
113/*
114* Calculate the Jacobi symbol
115*/
116int32_t jacobi(const BigInt& a, const BigInt& n) {
117 if(n.is_even() || n < 2) {
118 throw Invalid_Argument("jacobi: second argument must be odd and > 1");
119 }
120
121 BigInt x = a % n;
122 BigInt y = n;
123 int32_t J = 1;
124
125 while(y > 1) {
126 x %= y;
127 if(x > y / 2) {
128 x = y - x;
129 if(y % 4 == 3) {
130 J = -J;
131 }
132 }
133 if(x.is_zero()) {
134 return 0;
135 }
136
137 size_t shifts = low_zero_bits(x);
138 x >>= shifts;
139 if(shifts % 2) {
140 word y_mod_8 = y % 8;
141 if(y_mod_8 == 3 || y_mod_8 == 5) {
142 J = -J;
143 }
144 }
145
146 if(x % 4 == 3 && y % 4 == 3) {
147 J = -J;
148 }
149 std::swap(x, y);
150 }
151 return J;
152}
153
154/*
155* Square a BigInt
156*/
158 BigInt z = x;
160 z.square(ws);
161 return z;
162}
163
164/*
165* Return the number of 0 bits at the end of n
166*/
167size_t low_zero_bits(const BigInt& n) {
168 size_t low_zero = 0;
169
170 auto seen_nonempty_word = CT::Mask<word>::cleared();
171
172 for(size_t i = 0; i != n.size(); ++i) {
173 const word x = n.word_at(i);
174
175 // ctz(0) will return sizeof(word)
176 const size_t tz_x = ctz(x);
177
178 // if x > 0 we want to count tz_x in total but not any
179 // further words, so set the mask after the addition
180 low_zero += seen_nonempty_word.if_not_set_return(tz_x);
181
182 seen_nonempty_word |= CT::Mask<word>::expand(x);
183 }
184
185 // if we saw no words with x > 0 then n == 0 and the value we have
186 // computed is meaningless. Instead return BigInt::zero() in that case.
187 return static_cast<size_t>(seen_nonempty_word.if_set_return(low_zero));
188}
189
190/*
191* Calculate the GCD in constant time
192*/
193BigInt gcd(const BigInt& a, const BigInt& b) {
194 if(a.is_zero()) {
195 return abs(b);
196 }
197 if(b.is_zero()) {
198 return abs(a);
199 }
200
201 const size_t sz = std::max(a.sig_words(), b.sig_words());
202 auto u = BigInt::with_capacity(sz);
203 auto v = BigInt::with_capacity(sz);
204 u += a;
205 v += b;
206
208 v.const_time_poison();
209
210 u.set_sign(BigInt::Positive);
211 v.set_sign(BigInt::Positive);
212
213 // In the worst case we have two fully populated big ints. After right
214 // shifting so many times, we'll have reached the result for sure.
215 const size_t loop_cnt = u.bits() + v.bits();
216
217 using WordMask = CT::Mask<word>;
218
219 // This temporary is big enough to hold all intermediate results of the
220 // algorithm. No reallocation will happen during the loop.
221 // Note however, that `ct_cond_assign()` will invalidate the 'sig_words'
222 // cache, which _does not_ shrink the capacity of the underlying buffer.
223 auto tmp = BigInt::with_capacity(sz);
224 size_t factors_of_two = 0;
225 for(size_t i = 0; i != loop_cnt; ++i) {
226 auto both_odd = WordMask::expand(u.is_odd()) & WordMask::expand(v.is_odd());
227
228 // Subtract the smaller from the larger if both are odd
229 auto u_gt_v = WordMask::expand(bigint_cmp(u.data(), u.size(), v.data(), v.size()) > 0);
230 bigint_sub_abs(tmp.mutable_data(), u.data(), sz, v.data(), sz);
231 u.ct_cond_assign((u_gt_v & both_odd).as_bool(), tmp);
232 v.ct_cond_assign((~u_gt_v & both_odd).as_bool(), tmp);
233
234 const auto u_is_even = WordMask::expand(u.is_even());
235 const auto v_is_even = WordMask::expand(v.is_even());
236 BOTAN_DEBUG_ASSERT((u_is_even | v_is_even).as_bool());
237
238 // When both are even, we're going to eliminate a factor of 2.
239 // We have to reapply this factor to the final result.
240 factors_of_two += (u_is_even & v_is_even).if_set_return(1);
241
242 // remove one factor of 2, if u is even
243 bigint_shr2(tmp.mutable_data(), u.data(), sz, 0, 1);
244 u.ct_cond_assign(u_is_even.as_bool(), tmp);
245
246 // remove one factor of 2, if v is even
247 bigint_shr2(tmp.mutable_data(), v.data(), sz, 0, 1);
248 v.ct_cond_assign(v_is_even.as_bool(), tmp);
249 }
250
251 // The GCD (without factors of two) is either in u or v, the other one is
252 // zero. The non-zero variable _must_ be odd, because all factors of two were
253 // removed in the loop iterations above.
254 BOTAN_DEBUG_ASSERT(u.is_zero() || v.is_zero());
255 BOTAN_DEBUG_ASSERT(u.is_odd() || v.is_odd());
256
257 // make sure that the GCD (without factors of two) is in u
258 u.ct_cond_assign(u.is_even() /* .is_zero() would not be constant time */, v);
259
260 // re-apply the factors of two
261 u.ct_shift_left(factors_of_two);
262
263 u.const_time_unpoison();
264 v.const_time_unpoison();
265
266 return u;
267}
268
269/*
270* Calculate the LCM
271*/
272BigInt lcm(const BigInt& a, const BigInt& b) {
273 if(a == b) {
274 return a;
275 }
276
277 auto ab = a * b;
278 ab.set_sign(BigInt::Positive); // ignore the signs of a & b
279 const auto g = gcd(a, b);
280 return ct_divide(ab, g);
281}
282
283/*
284* Modular Exponentiation
285*/
286BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) {
287 if(mod.is_negative() || mod == 1) {
288 return BigInt::zero();
289 }
290
291 if(base.is_zero() || mod.is_zero()) {
292 if(exp.is_zero()) {
293 return BigInt::one();
294 }
295 return BigInt::zero();
296 }
297
298 Modular_Reducer reduce_mod(mod);
299
300 const size_t exp_bits = exp.bits();
301
302 if(mod.is_odd()) {
303 auto monty_params = std::make_shared<Montgomery_Params>(mod, reduce_mod);
304 return monty_exp(monty_params, reduce_mod.reduce(base), exp, exp_bits);
305 }
306
307 /*
308 Support for even modulus is just a convenience and not considered
309 cryptographically important, so this implementation is slow ...
310 */
311 BigInt accum = BigInt::one();
312 BigInt g = reduce_mod.reduce(base);
313 BigInt t;
314
315 for(size_t i = 0; i != exp_bits; ++i) {
316 t = reduce_mod.multiply(g, accum);
317 g = reduce_mod.square(g);
318 accum.ct_cond_assign(exp.get_bit(i), t);
319 }
320 return accum;
321}
322
324 if(C < 1) {
325 throw Invalid_Argument("is_perfect_square requires C >= 1");
326 }
327 if(C == 1) {
328 return BigInt::one();
329 }
330
331 const size_t n = C.bits();
332 const size_t m = (n + 1) / 2;
333 const BigInt B = C + BigInt::power_of_2(m);
334
335 BigInt X = BigInt::power_of_2(m) - 1;
336 BigInt X2 = (X * X);
337
338 for(;;) {
339 X = (X2 + C) / (2 * X);
340 X2 = (X * X);
341
342 if(X2 < B) {
343 break;
344 }
345 }
346
347 if(X2 == C) {
348 return X;
349 } else {
350 return BigInt::zero();
351 }
352}
353
354/*
355* Test for primality using Miller-Rabin
356*/
357bool is_prime(const BigInt& n, RandomNumberGenerator& rng, size_t prob, bool is_random) {
358 if(n == 2) {
359 return true;
360 }
361 if(n <= 1 || n.is_even()) {
362 return false;
363 }
364
365 const size_t n_bits = n.bits();
366
367 // Fast path testing for small numbers (<= 65521)
368 if(n_bits <= 16) {
369 const uint16_t num = static_cast<uint16_t>(n.word_at(0));
370
371 return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
372 }
373
374 Modular_Reducer mod_n(n);
375
376 if(rng.is_seeded()) {
377 const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
378
379 if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false) {
380 return false;
381 }
382
383 if(is_random) {
384 return true;
385 } else {
386 return is_lucas_probable_prime(n, mod_n);
387 }
388 } else {
389 return is_bailie_psw_probable_prime(n, mod_n);
390 }
391}
392
393} // namespace Botan
#define BOTAN_ASSERT_NOMSG(expr)
Definition assert.h:59
#define BOTAN_DEBUG_ASSERT(expr)
Definition assert.h:98
#define BOTAN_ARG_CHECK(expr, msg)
Definition assert.h:29
static BigInt zero()
Definition bigint.h:45
size_t sig_words() const
Definition bigint.h:584
bool is_odd() const
Definition bigint.h:416
void ct_cond_assign(bool predicate, const BigInt &other)
Definition bigint.cpp:487
size_t size() const
Definition bigint.h:578
static BigInt one()
Definition bigint.h:50
word word_at(size_t n) const
Definition bigint.h:518
size_t bits() const
Definition bigint.cpp:290
bool is_even() const
Definition bigint.h:410
static BigInt power_of_2(size_t n)
Definition bigint.h:738
void const_time_poison() const
Definition bigint.h:720
static BigInt from_s32(int32_t n)
Definition bigint.cpp:49
static BigInt from_word(word n)
Definition bigint.cpp:42
bool is_zero() const
Definition bigint.h:428
BigInt & square(secure_vector< word > &ws)
Definition big_ops2.cpp:183
bool is_negative() const
Definition bigint.h:528
bool get_bit(size_t n) const
Definition bigint.h:467
static BigInt with_capacity(size_t n)
Definition bigint.cpp:58
void set_sign(Sign sign)
Definition bigint.h:561
static Mask< T > expand(T v)
Definition ct_utils.h:109
static Mask< T > cleared()
Definition ct_utils.h:104
BigInt square(const BigInt &x) const
Definition reducer.h:43
BigInt multiply(const BigInt &x, const BigInt &y) const
Definition reducer.h:30
BigInt reduce(const BigInt &x) const
Definition reducer.cpp:37
virtual bool is_seeded() const =0
FE_25519 X
Definition ge.cpp:25
BigInt monty_exp(const std::shared_ptr< const Montgomery_Params > &params_p, const BigInt &g, const BigInt &k, size_t max_k_bits)
Definition monty_exp.h:41
BigInt power_mod(const BigInt &base, const BigInt &exp, const BigInt &mod)
Definition numthry.cpp:286
BigInt lcm(const BigInt &a, const BigInt &b)
Definition numthry.cpp:272
BigInt square(const BigInt &x)
Definition numthry.cpp:157
const uint16_t PRIMES[]
Definition primes.cpp:12
size_t low_zero_bits(const BigInt &n)
Definition numthry.cpp:167
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:355
const size_t PRIME_TABLE_SIZE
Definition numthry.h:167
bool is_prime(const BigInt &n, RandomNumberGenerator &rng, size_t prob, bool is_random)
Definition numthry.cpp:357
bool is_miller_rabin_probable_prime(const BigInt &n, const Modular_Reducer &mod_n, RandomNumberGenerator &rng, size_t test_iterations)
constexpr size_t ctz(T n)
Definition bit_ops.h:102
void ct_divide(const BigInt &x, const BigInt &y, BigInt &q_out, BigInt &r_out)
Definition divide.cpp:48
bool is_bailie_psw_probable_prime(const BigInt &n, const Modular_Reducer &mod_n)
Definition primality.cpp:89
BigInt monty_exp_vartime(const std::shared_ptr< const Montgomery_Params > &params_p, const BigInt &g, const BigInt &k)
Definition monty_exp.h:49
BigInt gcd(const BigInt &a, const BigInt &b)
Definition numthry.cpp:193
std::vector< T, secure_allocator< T > > secure_vector
Definition secmem.h:61
void bigint_shr2(word y[], const word x[], size_t x_size, size_t word_shift, size_t bit_shift)
Definition mp_core.h:431
BigInt sqrt_modulo_prime(const BigInt &a, const BigInt &p)
Definition numthry.cpp:26
BigInt is_perfect_square(const BigInt &C)
Definition numthry.cpp:323
size_t miller_rabin_test_iterations(size_t n_bits, size_t prob, bool random)
int32_t jacobi(const BigInt &a, const BigInt &n)
Definition numthry.cpp:116
int32_t bigint_cmp(const word x[], size_t x_size, const word y[], size_t y_size)
Definition mp_core.h:490
bool is_lucas_probable_prime(const BigInt &C, const Modular_Reducer &mod_C)
Definition primality.cpp:18