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