Botan 3.9.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/rng.h>
12#include <botan/internal/barrett.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
43 const Montgomery_Params monty_p(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)).value();
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).value();
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).value();
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)).value();
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 == 1) {
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
207 CT::poison_all(u, v);
208
209 u.set_sign(BigInt::Positive);
210 v.set_sign(BigInt::Positive);
211
212 // In the worst case we have two fully populated big ints. After right
213 // shifting so many times, we'll have reached the result for sure.
214 const size_t loop_cnt = u.bits() + v.bits();
215
216 // This temporary is big enough to hold all intermediate results of the
217 // algorithm. No reallocation will happen during the loop.
218 // Note however, that `ct_cond_assign()` will invalidate the 'sig_words'
219 // cache, which _does not_ shrink the capacity of the underlying buffer.
220 auto tmp = BigInt::with_capacity(sz);
221 secure_vector<word> ws(sz * 2);
222 size_t factors_of_two = 0;
223 for(size_t i = 0; i != loop_cnt; ++i) {
224 auto both_odd = CT::Mask<word>::expand_bool(u.is_odd()) & CT::Mask<word>::expand_bool(v.is_odd());
225
226 // Subtract the smaller from the larger if both are odd
227 auto u_gt_v = CT::Mask<word>::expand_bool(bigint_cmp(u._data(), u.size(), v._data(), v.size()) > 0);
228 bigint_sub_abs(tmp.mutable_data(), u._data(), v._data(), sz, ws.data());
229 u.ct_cond_assign((u_gt_v & both_odd).as_bool(), tmp);
230 v.ct_cond_assign((~u_gt_v & both_odd).as_bool(), tmp);
231
232 const auto u_is_even = CT::Mask<word>::expand_bool(u.is_even());
233 const auto v_is_even = CT::Mask<word>::expand_bool(v.is_even());
234 BOTAN_DEBUG_ASSERT((u_is_even | v_is_even).as_bool());
235
236 // When both are even, we're going to eliminate a factor of 2.
237 // We have to reapply this factor to the final result.
238 factors_of_two += (u_is_even & v_is_even).if_set_return(1);
239
240 // remove one factor of 2, if u is even
241 bigint_shr2(tmp.mutable_data(), u._data(), sz, 1);
242 u.ct_cond_assign(u_is_even.as_bool(), tmp);
243
244 // remove one factor of 2, if v is even
245 bigint_shr2(tmp.mutable_data(), v._data(), sz, 1);
246 v.ct_cond_assign(v_is_even.as_bool(), tmp);
247 }
248
249 // The GCD (without factors of two) is either in u or v, the other one is
250 // zero. The non-zero variable _must_ be odd, because all factors of two were
251 // removed in the loop iterations above.
252 BOTAN_DEBUG_ASSERT(u.is_zero() || v.is_zero());
253 BOTAN_DEBUG_ASSERT(u.is_odd() || v.is_odd());
254
255 // make sure that the GCD (without factors of two) is in u
256 u.ct_cond_assign(u.is_even() /* .is_zero() would not be constant time */, v);
257
258 // re-apply the factors of two
259 u.ct_shift_left(factors_of_two);
260
261 CT::unpoison_all(u, v);
262
263 return u;
264}
265
266/*
267* Calculate the LCM
268*/
269BigInt lcm(const BigInt& a, const BigInt& b) {
270 if(a == b) {
271 return a;
272 }
273
274 auto ab = a * b;
275 ab.set_sign(BigInt::Positive); // ignore the signs of a & b
276 const auto g = gcd(a, b);
277 return ct_divide(ab, g);
278}
279
280/*
281* Modular Exponentiation
282*/
283BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) {
284 if(mod.is_negative() || mod == 1) {
285 return BigInt::zero();
286 }
287
288 if(base.is_zero() || mod.is_zero()) {
289 if(exp.is_zero()) {
290 return BigInt::one();
291 }
292 return BigInt::zero();
293 }
294
296
297 const size_t exp_bits = exp.bits();
298
299 if(mod.is_odd()) {
300 const Montgomery_Params monty_params(mod, reduce_mod);
301 return monty_exp(monty_params, ct_modulo(base, mod), exp, exp_bits).value();
302 }
303
304 /*
305 Support for even modulus is just a convenience and not considered
306 cryptographically important, so this implementation is slow ...
307 */
308 BigInt accum = BigInt::one();
309 BigInt g = ct_modulo(base, mod);
310 BigInt t;
311
312 for(size_t i = 0; i != exp_bits; ++i) {
313 t = reduce_mod.multiply(g, accum);
314 g = reduce_mod.square(g);
315 accum.ct_cond_assign(exp.get_bit(i), t);
316 }
317 return accum;
318}
319
321 if(C < 1) {
322 throw Invalid_Argument("is_perfect_square requires C >= 1");
323 }
324 if(C == 1) {
325 return BigInt::one();
326 }
327
328 const size_t n = C.bits();
329 const size_t m = (n + 1) / 2;
330 const BigInt B = C + BigInt::power_of_2(m);
331
332 BigInt X = BigInt::power_of_2(m) - 1;
333 BigInt X2 = (X * X);
334
335 for(;;) {
336 X = (X2 + C) / (2 * X);
337 X2 = (X * X);
338
339 if(X2 < B) {
340 break;
341 }
342 }
343
344 if(X2 == C) {
345 return X;
346 } else {
347 return BigInt::zero();
348 }
349}
350
351/*
352* Test for primality using Miller-Rabin
353*/
354bool is_prime(const BigInt& n, RandomNumberGenerator& rng, size_t prob, bool is_random) {
355 if(n == 2) {
356 return true;
357 }
358 if(n <= 1 || n.is_even()) {
359 return false;
360 }
361
362 const size_t n_bits = n.bits();
363
364 // Fast path testing for small numbers (<= 65521)
365 if(n_bits <= 16) {
366 const uint16_t num = static_cast<uint16_t>(n.word_at(0));
367
368 return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
369 }
370
372
373 if(rng.is_seeded()) {
374 const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
375
376 if(!is_miller_rabin_probable_prime(n, mod_n, rng, t)) {
377 return false;
378 }
379
380 if(is_random) {
381 return true;
382 } else {
383 return is_lucas_probable_prime(n, mod_n);
384 }
385 } else {
386 return is_bailie_psw_probable_prime(n, mod_n);
387 }
388}
389
390} // namespace Botan
#define BOTAN_ASSERT_NOMSG(expr)
Definition assert.h:75
#define BOTAN_DEBUG_ASSERT(expr)
Definition assert.h:129
#define BOTAN_ARG_CHECK(expr, msg)
Definition assert.h:33
static Barrett_Reduction for_public_modulus(const BigInt &m)
Definition barrett.cpp:33
static Barrett_Reduction for_secret_modulus(const BigInt &m)
Definition barrett.cpp:22
static BigInt zero()
Definition bigint.h:49
size_t sig_words() const
Definition bigint.h:615
bool is_odd() const
Definition bigint.h:445
void ct_cond_assign(bool predicate, const BigInt &other)
Definition bigint.cpp:531
size_t size() const
Definition bigint.h:609
static BigInt one()
Definition bigint.h:54
word word_at(size_t n) const
Definition bigint.h:547
size_t bits() const
Definition bigint.cpp:311
bool is_even() const
Definition bigint.h:439
static BigInt power_of_2(size_t n)
Definition bigint.h:820
static BigInt from_s32(int32_t n)
Definition bigint.cpp:41
static BigInt from_word(word n)
Definition bigint.cpp:34
bool is_zero() const
Definition bigint.h:457
BigInt & square(secure_vector< word > &ws)
Definition big_ops2.cpp:175
bool is_negative() const
Definition bigint.h:559
bool get_bit(size_t n) const
Definition bigint.h:496
static BigInt with_capacity(size_t n)
Definition bigint.cpp:50
void set_sign(Sign sign)
Definition bigint.h:592
static constexpr Mask< T > expand(T v)
Definition ct_utils.h:420
static constexpr Mask< T > expand_bool(bool v)
Definition ct_utils.h:425
static constexpr Mask< T > cleared()
Definition ct_utils.h:415
BigInt value() const
Definition monty.cpp:245
virtual bool is_seeded() const =0
constexpr void poison_all(const Ts &... ts)
Definition ct_utils.h:199
constexpr void unpoison_all(const Ts &... ts)
Definition ct_utils.h:205
BigInt power_mod(const BigInt &base, const BigInt &exp, const BigInt &mod)
Definition numthry.cpp:283
constexpr void bigint_shr2(W y[], const W x[], size_t x_size, size_t shift)
Definition mp_core.h:370
Montgomery_Int monty_exp_vartime(const Montgomery_Params &params_p, const BigInt &g, const BigInt &k)
Definition monty_exp.h:55
bool is_lucas_probable_prime(const BigInt &C, const Barrett_Reduction &mod_C)
Definition primality.cpp:18
BigInt lcm(const BigInt &a, const BigInt &b)
Definition numthry.cpp:269
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:24
const size_t PRIME_TABLE_SIZE
Definition numthry.h:174
bool is_bailie_psw_probable_prime(const BigInt &n, const Barrett_Reduction &mod_n)
Definition primality.cpp:98
bool is_prime(const BigInt &n, RandomNumberGenerator &rng, size_t prob, bool is_random)
Definition numthry.cpp:354
BigInt ct_modulo(const BigInt &x, const BigInt &y)
Definition divide.cpp:192
constexpr int32_t bigint_cmp(const W x[], size_t x_size, const W y[], size_t y_size)
Definition mp_core.h:428
consteval std::array< W, N > reduce_mod(const std::array< W, XN > &x, const std::array< W, N > &p)
void ct_divide(const BigInt &x, const BigInt &y, BigInt &q_out, BigInt &r_out)
Definition divide.cpp:51
constexpr auto bigint_sub_abs(W z[], const W x[], const W y[], size_t N, W ws[]) -> CT::Mask< W >
Definition mp_core.h:281
BigInt gcd(const BigInt &a, const BigInt &b)
Definition numthry.cpp:193
std::vector< T, secure_allocator< T > > secure_vector
Definition secmem.h:69
BigInt sqrt_modulo_prime(const BigInt &a, const BigInt &p)
Definition numthry.cpp:26
BigInt is_perfect_square(const BigInt &C)
Definition numthry.cpp:320
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
bool is_miller_rabin_probable_prime(const BigInt &n, const Barrett_Reduction &mod_n, RandomNumberGenerator &rng, size_t test_iterations)
BOTAN_FORCE_INLINE constexpr size_t ctz(T n)
Definition bit_ops.h:96
std::conditional_t< HasNative64BitRegisters, std::uint64_t, uint32_t > word
Definition types.h:119
Montgomery_Int monty_exp(const Montgomery_Params &params_p, const BigInt &g, const BigInt &k, size_t max_k_bits)
Definition monty_exp.h:47