Botan  2.7.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/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 <algorithm>
18 
19 namespace Botan {
20 
21 /*
22 * Return the number of 0 bits at the end of n
23 */
24 size_t low_zero_bits(const BigInt& n)
25  {
26  size_t low_zero = 0;
27 
28  if(n.is_positive() && n.is_nonzero())
29  {
30  for(size_t i = 0; i != n.size(); ++i)
31  {
32  const word x = n.word_at(i);
33 
34  if(x)
35  {
36  low_zero += ctz(x);
37  break;
38  }
39  else
40  low_zero += BOTAN_MP_WORD_BITS;
41  }
42  }
43 
44  return low_zero;
45  }
46 
47 /*
48 * Calculate the GCD
49 */
50 BigInt gcd(const BigInt& a, const BigInt& b)
51  {
52  if(a.is_zero() || b.is_zero())
53  return 0;
54  if(a == 1 || b == 1)
55  return 1;
56 
57  BigInt X[2] = { a, b };
58  X[0].set_sign(BigInt::Positive);
59  X[1].set_sign(BigInt::Positive);
60 
61  const size_t shift = std::min(low_zero_bits(X[0]), low_zero_bits(X[1]));
62 
63  X[0] >>= shift;
64  X[1] >>= shift;
65 
66  while(X[0].is_nonzero())
67  {
68  X[0] >>= low_zero_bits(X[0]);
69  X[1] >>= low_zero_bits(X[1]);
70 
71  const uint8_t sel = static_cast<uint8_t>(X[0] >= X[1]);
72 
73  X[sel^1] -= X[sel];
74  X[sel^1] >>= 1;
75  }
76 
77  return (X[1] << shift);
78  }
79 
80 /*
81 * Calculate the LCM
82 */
83 BigInt lcm(const BigInt& a, const BigInt& b)
84  {
85  return ((a * b) / gcd(a, b));
86  }
87 
88 /*
89 Sets result to a^-1 * 2^k mod a
90 with n <= k <= 2n
91 Returns k
92 
93 "The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
94 https://citeseerx.ist.psu.edu/viewdoc/citations?doi=10.1.1.75.8377
95 
96 A const time implementation of this algorithm is described in
97 "Constant Time Modular Inversion" Joppe W. Bos
98 http://www.joppebos.com/files/CTInversion.pdf
99 */
101  const BigInt& a,
102  const BigInt& p)
103  {
104  size_t k = 0;
105 
106  BigInt u = p, v = a, r = 0, s = 1;
107 
108  while(v > 0)
109  {
110  if(u.is_even())
111  {
112  u >>= 1;
113  s <<= 1;
114  }
115  else if(v.is_even())
116  {
117  v >>= 1;
118  r <<= 1;
119  }
120  else if(u > v)
121  {
122  u -= v;
123  u >>= 1;
124  r += s;
125  s <<= 1;
126  }
127  else
128  {
129  v -= u;
130  v >>= 1;
131  s += r;
132  r <<= 1;
133  }
134 
135  ++k;
136  }
137 
138  if(r >= p)
139  {
140  r -= p;
141  }
142 
143  result = p - r;
144 
145  return k;
146  }
147 
149  {
150  BigInt r;
151  size_t k = almost_montgomery_inverse(r, a, p);
152 
153  for(size_t i = 0; i != k; ++i)
154  {
155  if(r.is_odd())
156  r += p;
157  r >>= 1;
158  }
159 
160  return r;
161  }
162 
164  {
165  if(n.is_negative() || mod.is_negative())
166  throw Invalid_Argument("ct_inverse_mod_odd_modulus: arguments must be non-negative");
167  if(mod < 3 || mod.is_even())
168  throw Invalid_Argument("Bad modulus to ct_inverse_mod_odd_modulus");
169  if(n >= mod)
170  throw Invalid_Argument("ct_inverse_mod_odd_modulus n >= mod not supported");
171 
172  /*
173  This uses a modular inversion algorithm designed by Niels Möller
174  and implemented in Nettle. The same algorithm was later also
175  adapted to GMP in mpn_sec_invert.
176 
177  It can be easily implemented in a way that does not depend on
178  secret branches or memory lookups, providing resistance against
179  some forms of side channel attack.
180 
181  There is also a description of the algorithm in Appendix 5 of "Fast
182  Software Polynomial Multiplication on ARM Processors using the NEON Engine"
183  by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
184  Dahab in LNCS 8182
185  https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
186 
187  Thanks to Niels for creating the algorithm, explaining some things
188  about it, and the reference to the paper.
189  */
190 
191  // todo allow this to be pre-calculated and passed in as arg
192  BigInt mp1o2 = (mod + 1) >> 1;
193 
194  const size_t mod_words = mod.sig_words();
195  BOTAN_ASSERT(mod_words > 0, "Not empty");
196 
197  BigInt a = n;
198  BigInt b = mod;
199  BigInt u = 1, v = 0;
200 
201  a.grow_to(mod_words);
202  u.grow_to(mod_words);
203  v.grow_to(mod_words);
204  mp1o2.grow_to(mod_words);
205 
209  secure_vector<word>& v_w = v.get_word_vector();
210 
211  CT::poison(a_w.data(), a_w.size());
212  CT::poison(b_w.data(), b_w.size());
213  CT::poison(u_w.data(), u_w.size());
214  CT::poison(v_w.data(), v_w.size());
215 
216  // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
217  size_t bits = 2 * mod.bits();
218 
219  while(bits--)
220  {
221  /*
222  const word odd = a.is_odd();
223  a -= odd * b;
224  const word underflow = a.is_negative();
225  b += a * underflow;
226  a.set_sign(BigInt::Positive);
227 
228  a >>= 1;
229 
230  if(underflow)
231  {
232  std::swap(u, v);
233  }
234 
235  u -= odd * v;
236  u += u.is_negative() * mod;
237 
238  const word odd_u = u.is_odd();
239 
240  u >>= 1;
241  u += mp1o2 * odd_u;
242  */
243 
244  const word odd_a = a_w[0] & 1;
245 
246  //if(odd_a) a -= b
247  word underflow = bigint_cnd_sub(odd_a, a_w.data(), b_w.data(), mod_words);
248 
249  //if(underflow) { b -= a; a = abs(a); swap(u, v); }
250  bigint_cnd_add(underflow, b_w.data(), a_w.data(), mod_words);
251  bigint_cnd_abs(underflow, a_w.data(), mod_words);
252  bigint_cnd_swap(underflow, u_w.data(), v_w.data(), mod_words);
253 
254  // a >>= 1
255  bigint_shr1(a_w.data(), mod_words, 0, 1);
256 
257  //if(odd_a) u -= v;
258  word borrow = bigint_cnd_sub(odd_a, u_w.data(), v_w.data(), mod_words);
259 
260  // if(borrow) u += p
261  bigint_cnd_add(borrow, u_w.data(), mod.data(), mod_words);
262 
263  const word odd_u = u_w[0] & 1;
264 
265  // u >>= 1
266  bigint_shr1(u_w.data(), mod_words, 0, 1);
267 
268  //if(odd_u) u += mp1o2;
269  bigint_cnd_add(odd_u, u_w.data(), mp1o2.data(), mod_words);
270  }
271 
272  CT::unpoison(a_w.data(), a_w.size());
273  CT::unpoison(b_w.data(), b_w.size());
274  CT::unpoison(u_w.data(), u_w.size());
275  CT::unpoison(v_w.data(), v_w.size());
276 
277  BOTAN_ASSERT(a.is_zero(), "A is zero");
278 
279  if(b != 1)
280  return 0;
281 
282  return v;
283  }
284 
285 /*
286 * Find the Modular Inverse
287 */
288 BigInt inverse_mod(const BigInt& n, const BigInt& mod)
289  {
290  if(mod.is_zero())
291  throw BigInt::DivideByZero();
292  if(mod.is_negative() || n.is_negative())
293  throw Invalid_Argument("inverse_mod: arguments must be non-negative");
294 
295  if(n.is_zero() || (n.is_even() && mod.is_even()))
296  return 0; // fast fail checks
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 
317  while(u.is_nonzero())
318  {
319  const size_t u_zero_bits = low_zero_bits(u);
320  u >>= u_zero_bits;
321  for(size_t i = 0; i != u_zero_bits; ++i)
322  {
323  if(A.is_odd() || B.is_odd())
324  { A += n; B -= mod; }
325  A >>= 1; B >>= 1;
326  }
327 
328  const size_t v_zero_bits = low_zero_bits(v);
329  v >>= v_zero_bits;
330  for(size_t i = 0; i != v_zero_bits; ++i)
331  {
332  if(C.is_odd() || D.is_odd())
333  { C += n; D -= mod; }
334  C >>= 1; D >>= 1;
335  }
336 
337  if(u >= v) { u -= v; A -= C; B -= D; }
338  else { v -= u; C -= A; D -= B; }
339  }
340 
341  if(v != 1)
342  return 0; // no modular inverse
343 
344  while(D.is_negative()) D += mod;
345  while(D >= mod) D -= mod;
346 
347  return D;
348  }
349 
350 word monty_inverse(word input)
351  {
352  if(input == 0)
353  throw Exception("monty_inverse: divide by zero");
354 
355  word b = input;
356  word x2 = 1, x1 = 0, y2 = 0, y1 = 1;
357 
358  // First iteration, a = n+1
359  word q = bigint_divop(1, 0, b);
360  word r = (MP_WORD_MAX - q*b) + 1;
361  word x = x2 - q*x1;
362  word y = y2 - q*y1;
363 
364  word a = b;
365  b = r;
366  x2 = x1;
367  x1 = x;
368  y2 = y1;
369  y1 = y;
370 
371  while(b > 0)
372  {
373  q = a / b;
374  r = a - q*b;
375  x = x2 - q*x1;
376  y = y2 - q*y1;
377 
378  a = b;
379  b = r;
380  x2 = x1;
381  x1 = x;
382  y2 = y1;
383  y1 = y;
384  }
385 
386  const word check = y2 * input;
387  BOTAN_ASSERT_EQUAL(check, 1, "monty_inverse result is inverse of input");
388 
389  // Now invert in addition space
390  y2 = (MP_WORD_MAX - y2) + 1;
391 
392  return y2;
393  }
394 
395 /*
396 * Modular Exponentiation
397 */
398 BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
399  {
400  if(mod.is_negative() || mod == 1)
401  {
402  return 0;
403  }
404 
405  if(base.is_zero() || mod.is_zero())
406  {
407  if(exp.is_zero())
408  return 1;
409  return 0;
410  }
411 
412  Power_Mod pow_mod(mod);
413 
414  /*
415  * Calling set_base before set_exponent means we end up using a
416  * minimal window. This makes sense given that here we know that any
417  * precomputation is wasted.
418  */
419 
420  if(base.is_negative())
421  {
422  pow_mod.set_base(-base);
423  pow_mod.set_exponent(exp);
424  if(exp.is_even())
425  return pow_mod.execute();
426  else
427  return (mod - pow_mod.execute());
428  }
429  else
430  {
431  pow_mod.set_base(base);
432  pow_mod.set_exponent(exp);
433  return pow_mod.execute();
434  }
435  }
436 
437 namespace {
438 
439 bool mr_witness(BigInt&& y,
440  const Modular_Reducer& reducer_n,
441  const BigInt& n_minus_1, size_t s)
442  {
443  if(y == 1 || y == n_minus_1)
444  return false;
445 
446  for(size_t i = 1; i != s; ++i)
447  {
448  y = reducer_n.square(y);
449 
450  if(y == 1) // found a non-trivial square root
451  return true;
452 
453  /*
454  -1 is the trivial square root of unity, so ``a`` is not a
455  witness for this number - give up
456  */
457  if(y == n_minus_1)
458  return false;
459  }
460 
461  return true; // is a witness
462  }
463 
464 size_t mr_test_iterations(size_t n_bits, size_t prob, bool random)
465  {
466  const size_t base = (prob + 2) / 2; // worst case 4^-t error rate
467 
468  /*
469  * If the candidate prime was maliciously constructed, we can't rely
470  * on arguments based on p being random.
471  */
472  if(random == false)
473  return base;
474 
475  /*
476  * For randomly chosen numbers we can use the estimates from
477  * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf
478  *
479  * These values are derived from the inequality for p(k,t) given on
480  * the second page.
481  */
482  if(prob <= 128)
483  {
484  if(n_bits >= 1536)
485  return 4; // < 2^-133
486  if(n_bits >= 1024)
487  return 6; // < 2^-133
488  if(n_bits >= 512)
489  return 12; // < 2^-129
490  if(n_bits >= 256)
491  return 29; // < 2^-128
492  }
493 
494  /*
495  If the user desires a smaller error probability than we have
496  precomputed error estimates for, just fall back to using the worst
497  case error rate.
498  */
499  return base;
500  }
501 
502 }
503 
504 /*
505 * Test for primality using Miller-Rabin
506 */
508  size_t prob, bool is_random)
509  {
510  if(n == 2)
511  return true;
512  if(n <= 1 || n.is_even())
513  return false;
514 
515  // Fast path testing for small numbers (<= 65521)
516  if(n <= PRIMES[PRIME_TABLE_SIZE-1])
517  {
518  const uint16_t num = static_cast<uint16_t>(n.word_at(0));
519 
520  return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
521  }
522 
523  const size_t test_iterations =
524  mr_test_iterations(n.bits(), prob, is_random && rng.is_seeded());
525 
526  const BigInt n_minus_1 = n - 1;
527  const size_t s = low_zero_bits(n_minus_1);
528  const BigInt nm1_s = n_minus_1 >> s;
529  const size_t n_bits = n.bits();
530 
531  const Modular_Reducer mod_n(n);
532  auto monty_n = std::make_shared<Montgomery_Params>(n, mod_n);
533 
534  const size_t powm_window = 4;
535 
536  for(size_t i = 0; i != test_iterations; ++i)
537  {
538  BigInt a;
539 
540  if(rng.is_seeded())
541  {
542  a = BigInt::random_integer(rng, 2, n_minus_1);
543  }
544  else
545  {
546  /*
547  * If passed a null RNG just use 2,3,5, ... as bases
548  *
549  * This is not ideal but in certain circumstances we need to
550  * test for primality but have no RNG available.
551  */
552  a = PRIMES[i];
553  }
554 
555  auto powm_a_n = monty_precompute(monty_n, a, powm_window);
556 
557  BigInt y = monty_execute(*powm_a_n, nm1_s, n_bits);
558 
559  if(mr_witness(std::move(y), mod_n, n_minus_1, s))
560  return false;
561  }
562 
563  return true;
564  }
565 
566 }
const size_t PRIME_TABLE_SIZE
Definition: numthry.h:269
void bigint_shr1(word x[], size_t x_size, size_t word_shift, size_t bit_shift)
Definition: mp_core.cpp:359
fe X
Definition: ge.cpp:27
size_t ctz(T n)
Definition: bit_ops.h:97
bool is_negative() const
Definition: bigint.h:460
const uint16_t PRIMES[]
Definition: primes.cpp:12
size_t low_zero_bits(const BigInt &n)
Definition: numthry.cpp:24
BigInt gcd(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:50
BigInt power_mod(const BigInt &base, const BigInt &exp, const BigInt &mod)
Definition: numthry.cpp:398
size_t bits() const
Definition: bigint.cpp:228
BigInt inverse_euclid(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:304
void bigint_cnd_abs(word cnd, word x[], size_t size)
Definition: mp_core.cpp:123
void bigint_cnd_swap(word cnd, word x[], word y[], size_t size)
Definition: mp_core.cpp:22
secure_vector< word > & get_word_vector()
Definition: bigint.h:553
word bigint_cnd_add(word cnd, word x[], const word y[], size_t size)
Definition: mp_core.cpp:39
bool is_zero() const
Definition: bigint.h:355
BigInt ct_inverse_mod_odd_modulus(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:163
BigInt normalized_montgomery_inverse(const BigInt &a, const BigInt &p)
Definition: numthry.cpp:148
bool is_prime(const BigInt &n, RandomNumberGenerator &rng, size_t prob, bool is_random)
Definition: numthry.cpp:507
void poison(const T *p, size_t n)
Definition: ct_utils.h:46
bool is_even() const
Definition: bigint.h:337
word bigint_divop(word n1, word n0, word d)
Definition: mp_core.cpp:482
bool is_nonzero() const
Definition: bigint.h:349
static BigInt random_integer(RandomNumberGenerator &rng, const BigInt &min, const BigInt &max)
Definition: big_rand.cpp:45
word word_at(size_t n) const
Definition: bigint.h:440
const word * data() const
Definition: bigint.h:551
#define BOTAN_ASSERT(expr, assertion_made)
Definition: assert.h:43
void set_exponent(const BigInt &exponent) const
Definition: pow_mod.cpp:79
BigInt lcm(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:83
#define BOTAN_ASSERT_EQUAL(expr1, expr2, assertion_made)
Definition: assert.h:69
bool is_odd() const
Definition: bigint.h:343
word bigint_cnd_sub(word cnd, word x[], const word y[], size_t size)
Definition: mp_core.cpp:69
BigInt inverse_mod(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:288
size_t size() const
Definition: bigint.h:513
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:151
Definition: alg_id.cpp:13
BigInt & square(secure_vector< word > &ws)
Definition: big_ops2.cpp:248
size_t sig_words() const
Definition: bigint.h:519
const word MP_WORD_MAX
Definition: mp_core.h:19
BigInt execute() const
Definition: pow_mod.cpp:92
void grow_to(size_t n)
Definition: bigint.cpp:302
virtual bool is_seeded() const =0
size_t almost_montgomery_inverse(BigInt &result, const BigInt &a, const BigInt &p)
Definition: numthry.cpp:100
void unpoison(const T *p, size_t n)
Definition: ct_utils.h:57
std::vector< T, secure_allocator< T > > secure_vector
Definition: secmem.h:88
void set_base(const BigInt &base) const
Definition: pow_mod.cpp:66
BigInt monty_execute(const Montgomery_Exponentation_State &precomputed_state, const BigInt &k, size_t max_k_bits)
Definition: monty_exp.cpp:159
word monty_inverse(word input)
Definition: numthry.cpp:350
bool is_positive() const
Definition: bigint.h:466