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