Botan  2.4.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 = 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 
161  /*
162  This uses a modular inversion algorithm designed by Niels Möller
163  and implemented in Nettle. The same algorithm was later also
164  adapted to GMP in mpn_sec_invert.
165 
166  It can be easily implemented in a way that does not depend on
167  secret branches or memory lookups, providing resistance against
168  some forms of side channel attack.
169 
170  There is also a description of the algorithm in Appendix 5 of "Fast
171  Software Polynomial Multiplication on ARM Processors using the NEON Engine"
172  by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
173  Dahab in LNCS 8182
174  https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
175 
176  Thanks to Niels for creating the algorithm, explaining some things
177  about it, and the reference to the paper.
178  */
179 
180  // todo allow this to be pre-calculated and passed in as arg
181  BigInt mp1o2 = (mod + 1) >> 1;
182 
183  const size_t mod_words = mod.sig_words();
184  BOTAN_ASSERT(mod_words > 0, "Not empty");
185 
186  BigInt a = n;
187  BigInt b = mod;
188  BigInt u = 1, v = 0;
189 
190  a.grow_to(mod_words);
191  u.grow_to(mod_words);
192  v.grow_to(mod_words);
193  mp1o2.grow_to(mod_words);
194 
198  secure_vector<word>& v_w = v.get_word_vector();
199 
200  CT::poison(a_w.data(), a_w.size());
201  CT::poison(b_w.data(), b_w.size());
202  CT::poison(u_w.data(), u_w.size());
203  CT::poison(v_w.data(), v_w.size());
204 
205  // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
206  size_t bits = 2 * mod.bits();
207 
208  while(bits--)
209  {
210  /*
211  const word odd = a.is_odd();
212  a -= odd * b;
213  const word underflow = a.is_negative();
214  b += a * underflow;
215  a.set_sign(BigInt::Positive);
216 
217  a >>= 1;
218 
219  if(underflow)
220  {
221  std::swap(u, v);
222  }
223 
224  u -= odd * v;
225  u += u.is_negative() * mod;
226 
227  const word odd_u = u.is_odd();
228 
229  u >>= 1;
230  u += mp1o2 * odd_u;
231  */
232 
233  const word odd_a = a_w[0] & 1;
234 
235  //if(odd_a) a -= b
236  word underflow = bigint_cnd_sub(odd_a, a_w.data(), b_w.data(), mod_words);
237 
238  //if(underflow) { b -= a; a = abs(a); swap(u, v); }
239  bigint_cnd_add(underflow, b_w.data(), a_w.data(), mod_words);
240  bigint_cnd_abs(underflow, a_w.data(), mod_words);
241  bigint_cnd_swap(underflow, u_w.data(), v_w.data(), mod_words);
242 
243  // a >>= 1
244  bigint_shr1(a_w.data(), mod_words, 0, 1);
245 
246  //if(odd_a) u -= v;
247  word borrow = bigint_cnd_sub(odd_a, u_w.data(), v_w.data(), mod_words);
248 
249  // if(borrow) u += p
250  bigint_cnd_add(borrow, u_w.data(), mod.data(), mod_words);
251 
252  const word odd_u = u_w[0] & 1;
253 
254  // u >>= 1
255  bigint_shr1(u_w.data(), mod_words, 0, 1);
256 
257  //if(odd_u) u += mp1o2;
258  bigint_cnd_add(odd_u, u_w.data(), mp1o2.data(), mod_words);
259  }
260 
261  CT::unpoison(a_w.data(), a_w.size());
262  CT::unpoison(b_w.data(), b_w.size());
263  CT::unpoison(u_w.data(), u_w.size());
264  CT::unpoison(v_w.data(), v_w.size());
265 
266  BOTAN_ASSERT(a.is_zero(), "A is zero");
267 
268  if(b != 1)
269  return 0;
270 
271  return v;
272  }
273 
274 /*
275 * Find the Modular Inverse
276 */
277 BigInt inverse_mod(const BigInt& n, const BigInt& mod)
278  {
279  if(mod.is_zero())
280  throw BigInt::DivideByZero();
281  if(mod.is_negative() || n.is_negative())
282  throw Invalid_Argument("inverse_mod: arguments must be non-negative");
283 
284  if(n.is_zero() || (n.is_even() && mod.is_even()))
285  return 0; // fast fail checks
286 
287  if(mod.is_odd())
288  return ct_inverse_mod_odd_modulus(n, mod);
289 
290  BigInt u = mod, v = n;
291  BigInt A = 1, B = 0, C = 0, D = 1;
292 
293  while(u.is_nonzero())
294  {
295  const size_t u_zero_bits = low_zero_bits(u);
296  u >>= u_zero_bits;
297  for(size_t i = 0; i != u_zero_bits; ++i)
298  {
299  if(A.is_odd() || B.is_odd())
300  { A += n; B -= mod; }
301  A >>= 1; B >>= 1;
302  }
303 
304  const size_t v_zero_bits = low_zero_bits(v);
305  v >>= v_zero_bits;
306  for(size_t i = 0; i != v_zero_bits; ++i)
307  {
308  if(C.is_odd() || D.is_odd())
309  { C += n; D -= mod; }
310  C >>= 1; D >>= 1;
311  }
312 
313  if(u >= v) { u -= v; A -= C; B -= D; }
314  else { v -= u; C -= A; D -= B; }
315  }
316 
317  if(v != 1)
318  return 0; // no modular inverse
319 
320  while(D.is_negative()) D += mod;
321  while(D >= mod) D -= mod;
322 
323  return D;
324  }
325 
326 word monty_inverse(word input)
327  {
328  if(input == 0)
329  throw Exception("monty_inverse: divide by zero");
330 
331  word b = input;
332  word x2 = 1, x1 = 0, y2 = 0, y1 = 1;
333 
334  // First iteration, a = n+1
335  word q = bigint_divop(1, 0, b);
336  word r = (MP_WORD_MAX - q*b) + 1;
337  word x = x2 - q*x1;
338  word y = y2 - q*y1;
339 
340  word a = b;
341  b = r;
342  x2 = x1;
343  x1 = x;
344  y2 = y1;
345  y1 = y;
346 
347  while(b > 0)
348  {
349  q = a / b;
350  r = a - q*b;
351  x = x2 - q*x1;
352  y = y2 - q*y1;
353 
354  a = b;
355  b = r;
356  x2 = x1;
357  x1 = x;
358  y2 = y1;
359  y1 = y;
360  }
361 
362  const word check = y2 * input;
363  BOTAN_ASSERT_EQUAL(check, 1, "monty_inverse result is inverse of input");
364 
365  // Now invert in addition space
366  y2 = (MP_WORD_MAX - y2) + 1;
367 
368  return y2;
369  }
370 
371 /*
372 * Modular Exponentiation
373 */
374 BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
375  {
376  Power_Mod pow_mod(mod);
377 
378  /*
379  * Calling set_base before set_exponent means we end up using a
380  * minimal window. This makes sense given that here we know that any
381  * precomputation is wasted.
382  */
383 
384  if(base.is_negative())
385  {
386  pow_mod.set_base(-base);
387  pow_mod.set_exponent(exp);
388  if(exp.is_even())
389  return pow_mod.execute();
390  else
391  return (mod - pow_mod.execute());
392  }
393  else
394  {
395  pow_mod.set_base(base);
396  pow_mod.set_exponent(exp);
397  return pow_mod.execute();
398  }
399  }
400 
401 namespace {
402 
403 bool mr_witness(BigInt&& y,
404  const Modular_Reducer& reducer_n,
405  const BigInt& n_minus_1, size_t s)
406  {
407  if(y == 1 || y == n_minus_1)
408  return false;
409 
410  for(size_t i = 1; i != s; ++i)
411  {
412  y = reducer_n.square(y);
413 
414  if(y == 1) // found a non-trivial square root
415  return true;
416 
417  if(y == n_minus_1) // -1, trivial square root, so give up
418  return false;
419  }
420 
421  return true; // fails Fermat test
422  }
423 
424 size_t mr_test_iterations(size_t n_bits, size_t prob, bool random)
425  {
426  const size_t base = (prob + 2) / 2; // worst case 4^-t error rate
427 
428  /*
429  * For randomly chosen numbers we can use the estimates from
430  * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf
431  *
432  * These values are derived from the inequality for p(k,t) given on
433  * the second page.
434  */
435  if(random && prob <= 80)
436  {
437  if(n_bits >= 1536)
438  return 2; // < 2^-89
439  if(n_bits >= 1024)
440  return 4; // < 2^-89
441  if(n_bits >= 512)
442  return 5; // < 2^-80
443  if(n_bits >= 256)
444  return 11; // < 2^-80
445  }
446 
447  return base;
448  }
449 
450 }
451 
452 /*
453 * Test for primaility using Miller-Rabin
454 */
456  size_t prob, bool is_random)
457  {
458  if(n == 2)
459  return true;
460  if(n <= 1 || n.is_even())
461  return false;
462 
463  // Fast path testing for small numbers (<= 65521)
464  if(n <= PRIMES[PRIME_TABLE_SIZE-1])
465  {
466  const uint16_t num = static_cast<uint16_t>(n.word_at(0));
467 
468  return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
469  }
470 
471  const size_t test_iterations = mr_test_iterations(n.bits(), prob, is_random);
472 
473  const BigInt n_minus_1 = n - 1;
474  const size_t s = low_zero_bits(n_minus_1);
475 
476  Fixed_Exponent_Power_Mod pow_mod(n_minus_1 >> s, n);
477  Modular_Reducer reducer(n);
478 
479  for(size_t i = 0; i != test_iterations; ++i)
480  {
481  const BigInt a = BigInt::random_integer(rng, 2, n_minus_1);
482  BigInt y = pow_mod(a);
483 
484  if(mr_witness(std::move(y), reducer, n_minus_1, s))
485  return false;
486  }
487 
488  return true;
489  }
490 
491 }
const size_t PRIME_TABLE_SIZE
Definition: numthry.h:240
void bigint_shr1(word x[], size_t x_size, size_t word_shift, size_t bit_shift)
Definition: mp_core.cpp:281
size_t ctz(T n)
Definition: bit_ops.h:97
bool is_negative() const
Definition: bigint.h:353
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:374
size_t bits() const
Definition: bigint.cpp:183
void bigint_cnd_abs(word cnd, word x[], size_t size)
Definition: mp_core.cpp:75
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:432
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:255
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:455
void poison(const T *p, size_t n)
Definition: ct_utils.h:46
bool is_even() const
Definition: bigint.h:237
word bigint_divop(word n1, word n0, word d)
Definition: mp_core.cpp:404
bool is_nonzero() const
Definition: bigint.h:249
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:340
const word * data() const
Definition: bigint.h:430
#define BOTAN_ASSERT(expr, assertion_made)
Definition: assert.h:29
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:55
bool is_odd() const
Definition: bigint.h:243
word bigint_cnd_sub(word cnd, word x[], const word y[], size_t size)
Definition: mp_core.cpp:61
BigInt inverse_mod(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:277
size_t size() const
Definition: bigint.h:392
Definition: alg_id.cpp:13
size_t sig_words() const
Definition: bigint.h:398
const word MP_WORD_MAX
Definition: mp_types.h:29
BigInt execute() const
Definition: pow_mod.cpp:92
void grow_to(size_t n)
Definition: bigint.cpp:260
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.cpp:214
BigInt square(const BigInt &x) const
Definition: reducer.h:39
word monty_inverse(word input)
Definition: numthry.cpp:326
bool is_positive() const
Definition: bigint.h:359