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
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); }
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
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;
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
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
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