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