Botan 2.13.0 Crypto and TLS for C++11

- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

numthry.cpp
Go to the documentation of this file.
1 /*
2 * Number Theory Functions
3 * (C) 1999-2011,2016,2018,2019 Jack Lloyd
4 *
5 * Botan is released under the Simplified BSD License (see license.txt)
6 */
7
8 #include <botan/numthry.h>
9 #include <botan/reducer.h>
10 #include <botan/monty.h>
11 #include <botan/divide.h>
12 #include <botan/rng.h>
13 #include <botan/internal/bit_ops.h>
14 #include <botan/internal/mp_core.h>
15 #include <botan/internal/ct_utils.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 {
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 * Return the number of 0 bits at the end of n
39 */
40 size_t low_zero_bits(const BigInt& n)
41 {
42 size_t low_zero = 0;
43
44 if(n.is_positive() && n.is_nonzero())
45 {
46 for(size_t i = 0; i != n.size(); ++i)
47 {
48 const word x = n.word_at(i);
49
50 if(x)
51 {
52 low_zero += ctz(x);
53 break;
54 }
55 else
56 low_zero += BOTAN_MP_WORD_BITS;
57 }
58 }
59
60 return low_zero;
61 }
62
63 /*
64 * Calculate the GCD
65 */
66 BigInt gcd(const BigInt& a, const BigInt& b)
67 {
68 if(a.is_zero() || b.is_zero())
69 return 0;
70 if(a == 1 || b == 1)
71 return 1;
72
73 BigInt f = a;
74 BigInt g = b;
75 f.set_sign(BigInt::Positive);
76 g.set_sign(BigInt::Positive);
77
78 const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g));
79
80 f >>= common2s;
81 g >>= common2s;
82
83 f.ct_cond_swap(f.is_even(), g);
84
85 int32_t delta = 1;
86
87 const size_t loop_cnt = 4 + 3*std::max(f.bits(), g.bits());
88
89 BigInt newg, t;
90 for(size_t i = 0; i != loop_cnt; ++i)
91 {
92 g.shrink_to_fit();
93 f.shrink_to_fit();
94 sub_abs(newg, f, g);
95
96 const bool need_swap = (g.is_odd() && delta > 0);
97
98 // if(need_swap) delta *= -1
99 delta *= CT::Mask<uint8_t>::expand(need_swap).select(0, 2) - 1;
100 f.ct_cond_swap(need_swap, g);
101 g.ct_cond_swap(need_swap, newg);
102
103 delta += 1;
104
105 t = g;
106 t += f;
107 g.ct_cond_assign(g.is_odd(), t);
108
109 g >>= 1;
110 }
111
112 f <<= common2s;
113
114 return f;
115 }
116
117 /*
118 * Calculate the LCM
119 */
120 BigInt lcm(const BigInt& a, const BigInt& b)
121 {
122 return ct_divide(a * b, gcd(a, b));
123 }
124
125 /*
126 Sets result to a^-1 * 2^k mod a
127 with n <= k <= 2n
128 Returns k
129
130 "The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
131 https://citeseerx.ist.psu.edu/viewdoc/citations?doi=10.1.1.75.8377
132
133 A const time implementation of this algorithm is described in
134 "Constant Time Modular Inversion" Joppe W. Bos
135 http://www.joppebos.com/files/CTInversion.pdf
136 */
137 size_t almost_montgomery_inverse(BigInt& result,
138 const BigInt& a,
139 const BigInt& p)
140 {
141 size_t k = 0;
142
143 BigInt u = p, v = a, r = 0, s = 1;
144
145 while(v > 0)
146 {
147 if(u.is_even())
148 {
149 u >>= 1;
150 s <<= 1;
151 }
152 else if(v.is_even())
153 {
154 v >>= 1;
155 r <<= 1;
156 }
157 else if(u > v)
158 {
159 u -= v;
160 u >>= 1;
161 r += s;
162 s <<= 1;
163 }
164 else
165 {
166 v -= u;
167 v >>= 1;
168 s += r;
169 r <<= 1;
170 }
171
172 ++k;
173 }
174
175 if(r >= p)
176 {
177 r -= p;
178 }
179
180 result = p - r;
181
182 return k;
183 }
184
185 BigInt normalized_montgomery_inverse(const BigInt& a, const BigInt& p)
186 {
187 BigInt r;
188 size_t k = almost_montgomery_inverse(r, a, p);
189
190 for(size_t i = 0; i != k; ++i)
191 {
192 if(r.is_odd())
193 r += p;
194 r >>= 1;
195 }
196
197 return r;
198 }
199
200 BigInt ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
201 {
202 if(n.is_negative() || mod.is_negative())
203 throw Invalid_Argument("ct_inverse_mod_odd_modulus: arguments must be non-negative");
204 if(mod < 3 || mod.is_even())
205 throw Invalid_Argument("Bad modulus to ct_inverse_mod_odd_modulus");
206 if(n >= mod)
207 throw Invalid_Argument("ct_inverse_mod_odd_modulus n >= mod not supported");
208
209 /*
210 This uses a modular inversion algorithm designed by Niels Möller
211 and implemented in Nettle. The same algorithm was later also
212 adapted to GMP in mpn_sec_invert.
213
214 It can be easily implemented in a way that does not depend on
215 secret branches or memory lookups, providing resistance against
216 some forms of side channel attack.
217
218 There is also a description of the algorithm in Appendix 5 of "Fast
219 Software Polynomial Multiplication on ARM Processors using the NEON Engine"
220 by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
221 Dahab in LNCS 8182
222 https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
223
224 Thanks to Niels for creating the algorithm, explaining some things
225 about it, and the reference to the paper.
226 */
227
228 // todo allow this to be pre-calculated and passed in as arg
229 BigInt mp1o2 = (mod + 1) >> 1;
230
231 const size_t mod_words = mod.sig_words();
232 BOTAN_ASSERT(mod_words > 0, "Not empty");
233
234 BigInt a = n;
235 BigInt b = mod;
236 BigInt u = 1, v = 0;
237
238 a.grow_to(mod_words);
239 u.grow_to(mod_words);
240 v.grow_to(mod_words);
241 mp1o2.grow_to(mod_words);
242
243 secure_vector<word>& a_w = a.get_word_vector();
244 secure_vector<word>& b_w = b.get_word_vector();
245 secure_vector<word>& u_w = u.get_word_vector();
246 secure_vector<word>& v_w = v.get_word_vector();
247
248 CT::poison(a_w.data(), a_w.size());
249 CT::poison(b_w.data(), b_w.size());
250 CT::poison(u_w.data(), u_w.size());
251 CT::poison(v_w.data(), v_w.size());
252
253 // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
254 size_t bits = 2 * mod.bits();
255
256 while(bits--)
257 {
258 /*
259 const word odd = a.is_odd();
260 a -= odd * b;
261 const word underflow = a.is_negative();
262 b += a * underflow;
263 a.set_sign(BigInt::Positive);
264
265 a >>= 1;
266
267 if(underflow)
268 {
269 std::swap(u, v);
270 }
271
272 u -= odd * v;
273 u += u.is_negative() * mod;
274
275 const word odd_u = u.is_odd();
276
277 u >>= 1;
278 u += mp1o2 * odd_u;
279 */
280
281 const word odd_a = a_w[0] & 1;
282
283 //if(odd_a) a -= b
284 word underflow = bigint_cnd_sub(odd_a, a_w.data(), b_w.data(), mod_words);
285
286 //if(underflow) { b -= a; a = abs(a); swap(u, v); }
287 bigint_cnd_add(underflow, b_w.data(), a_w.data(), mod_words);
288 bigint_cnd_abs(underflow, a_w.data(), mod_words);
289 bigint_cnd_swap(underflow, u_w.data(), v_w.data(), mod_words);
290
291 // a >>= 1
292 bigint_shr1(a_w.data(), mod_words, 0, 1);
293
294 //if(odd_a) u -= v;
295 word borrow = bigint_cnd_sub(odd_a, u_w.data(), v_w.data(), mod_words);
296
297 // if(borrow) u += p
298 bigint_cnd_add(borrow, u_w.data(), mod.data(), mod_words);
299
300 const word odd_u = u_w[0] & 1;
301
302 // u >>= 1
303 bigint_shr1(u_w.data(), mod_words, 0, 1);
304
305 //if(odd_u) u += mp1o2;
306 bigint_cnd_add(odd_u, u_w.data(), mp1o2.data(), mod_words);
307 }
308
309 CT::unpoison(a_w.data(), a_w.size());
310 CT::unpoison(b_w.data(), b_w.size());
311 CT::unpoison(u_w.data(), u_w.size());
312 CT::unpoison(v_w.data(), v_w.size());
313
314 BOTAN_ASSERT(a.is_zero(), "A is zero");
315
316 if(b != 1)
317 return 0;
318
319 return v;
320 }
321
322 /*
323 * Find the Modular Inverse
324 */
325 BigInt inverse_mod(const BigInt& n, const BigInt& mod)
326 {
327 if(mod.is_zero())
328 throw BigInt::DivideByZero();
329 if(mod.is_negative() || n.is_negative())
330 throw Invalid_Argument("inverse_mod: arguments must be non-negative");
331 if(n.is_zero())
332 return 0;
333
334 if(mod.is_odd() && n < mod)
335 return ct_inverse_mod_odd_modulus(n, mod);
336
337 return inverse_euclid(n, mod);
338 }
339
340 BigInt inverse_euclid(const BigInt& n, const BigInt& mod)
341 {
342 if(mod.is_zero())
343 throw BigInt::DivideByZero();
344 if(mod.is_negative() || n.is_negative())
345 throw Invalid_Argument("inverse_mod: arguments must be non-negative");
346
347 if(n.is_zero() || (n.is_even() && mod.is_even()))
348 return 0; // fast fail checks
349
350 BigInt u = mod, v = n;
351 BigInt A = 1, B = 0, C = 0, D = 1;
352 BigInt T0, T1, T2;
353
354 while(u.is_nonzero())
355 {
356 const size_t u_zero_bits = low_zero_bits(u);
357 u >>= u_zero_bits;
358
359 const size_t v_zero_bits = low_zero_bits(v);
360 v >>= v_zero_bits;
361
362 const bool u_gte_v = (u >= v);
363
364 for(size_t i = 0; i != u_zero_bits; ++i)
365 {
366 const bool needs_adjust = A.is_odd() || B.is_odd();
367
368 T0 = A + n;
369 T1 = B - mod;
370
371 A.ct_cond_assign(needs_adjust, T0);
372 B.ct_cond_assign(needs_adjust, T1);
373
374 A >>= 1;
375 B >>= 1;
376 }
377
378 for(size_t i = 0; i != v_zero_bits; ++i)
379 {
380 const bool needs_adjust = C.is_odd() || D.is_odd();
381 T0 = C + n;
382 T1 = D - mod;
383
384 C.ct_cond_assign(needs_adjust, T0);
385 D.ct_cond_assign(needs_adjust, T1);
386
387 C >>= 1;
388 D >>= 1;
389 }
390
391 T0 = u - v;
392 T1 = A - C;
393 T2 = B - D;
394
395 T0.cond_flip_sign(!u_gte_v);
396 T1.cond_flip_sign(!u_gte_v);
397 T2.cond_flip_sign(!u_gte_v);
398
399 u.ct_cond_assign(u_gte_v, T0);
400 A.ct_cond_assign(u_gte_v, T1);
401 B.ct_cond_assign(u_gte_v, T2);
402
403 v.ct_cond_assign(!u_gte_v, T0);
404 C.ct_cond_assign(!u_gte_v, T1);
405 D.ct_cond_assign(!u_gte_v, T2);
406 }
407
408 if(v != 1)
409 return 0; // no modular inverse
410
411 while(D.is_negative())
412 D += mod;
413 while(D >= mod)
414 D -= mod;
415
416 return D;
417 }
418
419 word monty_inverse(word a)
420 {
421 if(a % 2 == 0)
422 throw Invalid_Argument("monty_inverse only valid for odd integers");
423
424 /*
425 * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
426 * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
427 */
428
429 word b = 1;
430 word r = 0;
431
432 for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
433 {
434 const word bi = b % 2;
435 r >>= 1;
436 r += bi << (BOTAN_MP_WORD_BITS - 1);
437
438 b -= a * bi;
439 b >>= 1;
440 }
441
442 // Now invert in addition space
443 r = (MP_WORD_MAX - r) + 1;
444
445 return r;
446 }
447
448 /*
449 * Modular Exponentiation
450 */
451 BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
452 {
453 if(mod.is_negative() || mod == 1)
454 {
455 return 0;
456 }
457
458 if(base.is_zero() || mod.is_zero())
459 {
460 if(exp.is_zero())
461 return 1;
462 return 0;
463 }
464
465 Modular_Reducer reduce_mod(mod);
466
467 const size_t exp_bits = exp.bits();
468
469 if(mod.is_odd())
470 {
471 const size_t powm_window = 4;
472
473 auto monty_mod = std::make_shared<Montgomery_Params>(mod, reduce_mod);
474 auto powm_base_mod = monty_precompute(monty_mod, reduce_mod.reduce(base), powm_window);
475 return monty_execute(*powm_base_mod, exp, exp_bits);
476 }
477
478 /*
479 Support for even modulus is just a convenience and not considered
480 cryptographically important, so this implementation is slow ...
481 */
482 BigInt accum = 1;
483 BigInt g = reduce_mod.reduce(base);
484 BigInt t;
485
486 for(size_t i = 0; i != exp_bits; ++i)
487 {
488 t = reduce_mod.multiply(g, accum);
489 g = reduce_mod.square(g);
490 accum.ct_cond_assign(exp.get_bit(i), t);
491 }
492 return accum;
493 }
494
495
496 BigInt is_perfect_square(const BigInt& C)
497 {
498 if(C < 1)
499 throw Invalid_Argument("is_perfect_square requires C >= 1");
500 if(C == 1)
501 return 1;
502
503 const size_t n = C.bits();
504 const size_t m = (n + 1) / 2;
505 const BigInt B = C + BigInt::power_of_2(m);
506
507 BigInt X = BigInt::power_of_2(m) - 1;
508 BigInt X2 = (X*X);
509
510 for(;;)
511 {
512 X = (X2 + C) / (2*X);
513 X2 = (X*X);
514
515 if(X2 < B)
516 break;
517 }
518
519 if(X2 == C)
520 return X;
521 else
522 return 0;
523 }
524
525 /*
526 * Test for primality using Miller-Rabin
527 */
528 bool is_prime(const BigInt& n,
529 RandomNumberGenerator& rng,
530 size_t prob,
531 bool is_random)
532 {
533 if(n == 2)
534 return true;
535 if(n <= 1 || n.is_even())
536 return false;
537
538 const size_t n_bits = n.bits();
539
540 // Fast path testing for small numbers (<= 65521)
541 if(n_bits <= 16)
542 {
543 const uint16_t num = static_cast<uint16_t>(n.word_at(0));
544
545 return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
546 }
547
548 Modular_Reducer mod_n(n);
549
550 if(rng.is_seeded())
551 {
552 const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
553
554 if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false)
555 return false;
556
557 return is_lucas_probable_prime(n, mod_n);
558 }
559 else
560 {
561 return is_bailie_psw_probable_prime(n, mod_n);
562 }
563 }
564
565 }
Botan::PRIME_TABLE_SIZE
const size_t PRIME_TABLE_SIZE
Definition: numthry.h:276
Botan::bigint_shr1
void bigint_shr1(word x[], size_t x_size, size_t word_shift, size_t bit_shift)
Definition: mp_core.h:427
X
fe X
Definition: ge.cpp:27
Botan::BigInt::word_at
word word_at(size_t n) const
Definition: bigint.h:506
Botan::ctz
size_t ctz(T n)
Definition: bit_ops.h:99
Botan::BigInt::sig_words
size_t sig_words() const
Definition: bigint.h:584
Botan::Invalid_Argument
Definition: exceptn.h:136
Botan::is_lucas_probable_prime
bool is_lucas_probable_prime(const BigInt &C, const Modular_Reducer &mod_C)
Definition: primality.cpp:17
Botan::BigInt::is_odd
bool is_odd() const
Definition: bigint.h:407
Botan::BigInt::is_even
bool is_even() const
Definition: bigint.h:401
Botan::PRIMES
const uint16_t PRIMES[]
Definition: primes.cpp:12
Botan::low_zero_bits
size_t low_zero_bits(const BigInt &n)
Definition: numthry.cpp:40
Botan::BigInt
Definition: bigint.h:24
Botan::gcd
BigInt gcd(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:66
Botan::BigInt::Positive
Definition: bigint.h:35
Botan::miller_rabin_test_iterations
size_t miller_rabin_test_iterations(size_t n_bits, size_t prob, bool random)
Definition: primality.cpp:168
Botan::power_mod
BigInt power_mod(const BigInt &base, const BigInt &exp, const BigInt &mod)
Definition: numthry.cpp:451
Botan::RandomNumberGenerator
Definition: rng.h:24
Botan::inverse_euclid
BigInt inverse_euclid(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:340
Botan::BigInt::get_word_vector
secure_vector< word > & get_word_vector()
Definition: bigint.h:623
Botan::ct_inverse_mod_odd_modulus
BigInt ct_inverse_mod_odd_modulus(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:200
Botan::normalized_montgomery_inverse
BigInt normalized_montgomery_inverse(const BigInt &a, const BigInt &p)
Definition: numthry.cpp:185
Botan::is_prime
bool is_prime(const BigInt &n, RandomNumberGenerator &rng, size_t prob, bool is_random)
Definition: numthry.cpp:528
Botan::BigInt::is_negative
bool is_negative() const
Definition: bigint.h:525
Botan::CT::poison
void poison(const T *p, size_t n)
Definition: ct_utils.h:48
Botan::BigInt::size
size_t size() const
Definition: bigint.h:578
Botan::BigInt::is_nonzero
bool is_nonzero() const
Definition: bigint.h:413
BOTAN_ASSERT
#define BOTAN_ASSERT(expr, assertion_made)
Definition: assert.h:55
Botan::is_bailie_psw_probable_prime
bool is_bailie_psw_probable_prime(const BigInt &n, const Modular_Reducer &mod_n)
Definition: primality.cpp:95
Botan::BigInt::bits
size_t bits() const
Definition: bigint.cpp:288
Botan::CT::Mask::expand
static Mask< T > expand(T v)
Definition: ct_utils.h:123
Botan::bigint_cnd_sub
word bigint_cnd_sub(word cnd, word x[], size_t x_size, const word y[], size_t y_size)
Definition: mp_core.h:88
Botan::secure_vector
std::vector< T, secure_allocator< T >> secure_vector
Definition: secmem.h:65
Botan::BigInt::DivideByZero
Definition: bigint.h:43
Botan::bigint_sub_abs
CT::Mask< word > bigint_sub_abs(word z[], const word x[], const word y[], size_t N, word ws[])
Definition: mp_core.h:377
Botan::lcm
BigInt lcm(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:120
Botan::Modular_Reducer::multiply
BigInt multiply(const BigInt &x, const BigInt &y) const
Definition: reducer.h:31
Botan::ct_divide
void ct_divide(const BigInt &x, const BigInt &y, BigInt &q_out, BigInt &r_out)
Definition: divide.cpp:52
Botan::is_perfect_square
BigInt is_perfect_square(const BigInt &C)
Definition: numthry.cpp:496
Botan::BigInt::ct_cond_swap
void ct_cond_swap(bool predicate, BigInt &other)
Definition: bigint.cpp:447
Botan::bigint_cnd_swap
void bigint_cnd_swap(word cnd, word x[], word y[], size_t size)
Definition: mp_core.h:29
Botan::bigint_cnd_add
word bigint_cnd_add(word cnd, word x[], word x_size, const word y[], size_t y_size)
Definition: mp_core.h:42
Botan::inverse_mod
BigInt inverse_mod(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:325
Botan::BigInt::data
const word * data() const
Definition: bigint.h:618
Botan::monty_precompute
std::shared_ptr< const Montgomery_Exponentation_State > monty_precompute(std::shared_ptr< const Montgomery_Params > params, const BigInt &g, size_t window_bits, bool const_time)
Definition: monty_exp.cpp:157
Botan::BigInt::power_of_2
static BigInt power_of_2(size_t n)
Definition: bigint.h:751
Botan
Definition: alg_id.cpp:13
Botan::Modular_Reducer::reduce
BigInt reduce(const BigInt &x) const
Definition: reducer.cpp:37
Botan::BigInt::cond_flip_sign
void cond_flip_sign(bool predicate)
Definition: bigint.cpp:456
Botan::MP_WORD_MAX
const word MP_WORD_MAX
Definition: mp_core.h:22
Botan::monty_inverse
word monty_inverse(word a)
Definition: numthry.cpp:419
Botan::RandomNumberGenerator::is_seeded
virtual bool is_seeded() const =0
Botan::almost_montgomery_inverse
size_t almost_montgomery_inverse(BigInt &result, const BigInt &a, const BigInt &p)
Definition: numthry.cpp:137
Botan::BigInt::grow_to
void grow_to(size_t n) const
Definition: bigint.h:634
Botan::BigInt::is_positive
bool is_positive() const
Definition: bigint.h:531
Botan::is_miller_rabin_probable_prime
bool is_miller_rabin_probable_prime(const BigInt &n, const Modular_Reducer &mod_n, RandomNumberGenerator &rng, size_t test_iterations)
Definition: primality.cpp:146
Botan::CT::unpoison
void unpoison(const T *p, size_t n)
Definition: ct_utils.h:59
Botan::BigInt::get_bit
bool get_bit(size_t n) const
Definition: bigint.h:463
Botan::bigint_cnd_abs
void bigint_cnd_abs(word cnd, word x[], size_t size)
Definition: mp_core.h:212
Botan::BigInt::is_zero
bool is_zero() const
Definition: bigint.h:419
Botan::BigInt::shrink_to_fit
void shrink_to_fit(size_t min_size=0)
Definition: bigint.h:640
Botan::monty_execute
BigInt monty_execute(const Montgomery_Exponentation_State &precomputed_state, const BigInt &k, size_t max_k_bits)
Definition: monty_exp.cpp:165
Botan::Modular_Reducer::square
BigInt square(const BigInt &x) const
Definition: reducer.h:39
Botan::BigInt::set_sign
void set_sign(Sign sign)
Definition: bigint.h:561
Botan::Modular_Reducer
Definition: reducer.h:18
Botan::BigInt::ct_cond_assign
void ct_cond_assign(bool predicate, const BigInt &other)
Definition: bigint.cpp:469