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