Coverage Report

Created: 2020-02-14 15:38

/src/botan/src/lib/math/numbertheory/numthry.cpp
Line
Count
Source (jump to first uncovered line)
1
/*
2
* Number Theory Functions
3
* (C) 1999-2011,2016,2018,2019 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/monty.h>
11
#include <botan/divide.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 <botan/internal/primality.h>
18
#include <algorithm>
19
20
namespace Botan {
21
22
namespace {
23
24
void sub_abs(BigInt& z, const BigInt& x, const BigInt& y)
25
0
   {
26
0
   const size_t x_sw = x.sig_words();
27
0
   const size_t y_sw = y.sig_words();
28
0
   z.resize(std::max(x_sw, y_sw));
29
0
30
0
   bigint_sub_abs(z.mutable_data(),
31
0
                  x.data(), x_sw,
32
0
                  y.data(), y_sw);
33
0
   }
34
35
}
36
37
/*
38
* Return the number of 0 bits at the end of n
39
*/
40
size_t low_zero_bits(const BigInt& n)
41
3.27M
   {
42
3.27M
   size_t low_zero = 0;
43
3.27M
44
3.27M
   if(n.is_positive() && n.is_nonzero())
45
3.27M
      {
46
3.27M
      for(size_t i = 0; i != n.size(); ++i)
47
3.27M
         {
48
3.27M
         const word x = n.word_at(i);
49
3.27M
50
3.27M
         if(x)
51
3.27M
            {
52
3.27M
            low_zero += ctz(x);
53
3.27M
            break;
54
3.27M
            }
55
3.21k
         else
56
3.21k
            low_zero += BOTAN_MP_WORD_BITS;
57
3.27M
         }
58
3.27M
      }
59
3.27M
60
3.27M
   return low_zero;
61
3.27M
   }
62
63
/*
64
* Calculate the GCD
65
*/
66
BigInt gcd(const BigInt& a, const BigInt& b)
67
0
   {
68
0
   if(a.is_zero() || b.is_zero())
69
0
      return 0;
70
0
   if(a == 1 || b == 1)
71
0
      return 1;
72
0
73
0
   BigInt f = a;
74
0
   BigInt g = b;
75
0
   f.set_sign(BigInt::Positive);
76
0
   g.set_sign(BigInt::Positive);
77
0
78
0
   const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g));
79
0
80
0
   f >>= common2s;
81
0
   g >>= common2s;
82
0
83
0
   f.ct_cond_swap(f.is_even(), g);
84
0
85
0
   int32_t delta = 1;
86
0
87
0
   const size_t loop_cnt = 4 + 3*std::max(f.bits(), g.bits());
88
0
89
0
   BigInt newg, t;
90
0
   for(size_t i = 0; i != loop_cnt; ++i)
91
0
      {
92
0
      g.shrink_to_fit();
93
0
      f.shrink_to_fit();
94
0
      sub_abs(newg, f, g);
95
0
96
0
      const bool need_swap = (g.is_odd() && delta > 0);
97
0
98
0
      // if(need_swap) delta *= -1
99
0
      delta *= CT::Mask<uint8_t>::expand(need_swap).select(0, 2) - 1;
100
0
      f.ct_cond_swap(need_swap, g);
101
0
      g.ct_cond_swap(need_swap, newg);
102
0
103
0
      delta += 1;
104
0
105
0
      t = g;
106
0
      t += f;
107
0
      g.ct_cond_assign(g.is_odd(), t);
108
0
109
0
      g >>= 1;
110
0
      }
111
0
112
0
   f <<= common2s;
113
0
114
0
   return f;
115
0
   }
116
117
/*
118
* Calculate the LCM
119
*/
120
BigInt lcm(const BigInt& a, const BigInt& b)
121
0
   {
122
0
   return ct_divide(a * b, gcd(a, b));
123
0
   }
124
125
/*
126
Sets result to a^-1 * 2^k mod a
127
with n <= k <= 2n
128
Returns k
129
130
"The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
131
https://citeseerx.ist.psu.edu/viewdoc/citations?doi=10.1.1.75.8377
132
133
A const time implementation of this algorithm is described in
134
"Constant Time Modular Inversion" Joppe W. Bos
135
http://www.joppebos.com/files/CTInversion.pdf
136
*/
137
size_t almost_montgomery_inverse(BigInt& result,
138
                                 const BigInt& a,
139
                                 const BigInt& p)
140
0
   {
141
0
   size_t k = 0;
142
0
143
0
   BigInt u = p, v = a, r = 0, s = 1;
144
0
145
0
   while(v > 0)
146
0
      {
147
0
      if(u.is_even())
148
0
         {
149
0
         u >>= 1;
150
0
         s <<= 1;
151
0
         }
152
0
      else if(v.is_even())
153
0
         {
154
0
         v >>= 1;
155
0
         r <<= 1;
156
0
         }
157
0
      else if(u > v)
158
0
         {
159
0
         u -= v;
160
0
         u >>= 1;
161
0
         r += s;
162
0
         s <<= 1;
163
0
         }
164
0
      else
165
0
         {
166
0
         v -= u;
167
0
         v >>= 1;
168
0
         s += r;
169
0
         r <<= 1;
170
0
         }
171
0
172
0
      ++k;
173
0
      }
174
0
175
0
   if(r >= p)
176
0
      {
177
0
      r -= p;
178
0
      }
179
0
180
0
   result = p - r;
181
0
182
0
   return k;
183
0
   }
184
185
BigInt normalized_montgomery_inverse(const BigInt& a, const BigInt& p)
186
0
   {
187
0
   BigInt r;
188
0
   size_t k = almost_montgomery_inverse(r, a, p);
189
0
190
0
   for(size_t i = 0; i != k; ++i)
191
0
      {
192
0
      if(r.is_odd())
193
0
         r += p;
194
0
      r >>= 1;
195
0
      }
196
0
197
0
   return r;
198
0
   }
199
200
BigInt ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
201
50.5k
   {
202
50.5k
   if(n.is_negative() || mod.is_negative())
203
0
      throw Invalid_Argument("ct_inverse_mod_odd_modulus: arguments must be non-negative");
204
50.5k
   if(mod < 3 || mod.is_even())
205
0
      throw Invalid_Argument("Bad modulus to ct_inverse_mod_odd_modulus");
206
50.5k
   if(n >= mod)
207
0
      throw Invalid_Argument("ct_inverse_mod_odd_modulus n >= mod not supported");
208
50.5k
209
50.5k
   /*
210
50.5k
   This uses a modular inversion algorithm designed by Niels Möller
211
50.5k
   and implemented in Nettle. The same algorithm was later also
212
50.5k
   adapted to GMP in mpn_sec_invert.
213
50.5k
214
50.5k
   It can be easily implemented in a way that does not depend on
215
50.5k
   secret branches or memory lookups, providing resistance against
216
50.5k
   some forms of side channel attack.
217
50.5k
218
50.5k
   There is also a description of the algorithm in Appendix 5 of "Fast
219
50.5k
   Software Polynomial Multiplication on ARM Processors using the NEON Engine"
220
50.5k
   by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
221
50.5k
   Dahab in LNCS 8182
222
50.5k
      https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
223
50.5k
224
50.5k
   Thanks to Niels for creating the algorithm, explaining some things
225
50.5k
   about it, and the reference to the paper.
226
50.5k
   */
227
50.5k
228
50.5k
   // todo allow this to be pre-calculated and passed in as arg
229
50.5k
   BigInt mp1o2 = (mod + 1) >> 1;
230
50.5k
231
50.5k
   const size_t mod_words = mod.sig_words();
232
50.5k
   BOTAN_ASSERT(mod_words > 0, "Not empty");
233
50.5k
234
50.5k
   BigInt a = n;
235
50.5k
   BigInt b = mod;
236
50.5k
   BigInt u = 1, v = 0;
237
50.5k
238
50.5k
   a.grow_to(mod_words);
239
50.5k
   u.grow_to(mod_words);
240
50.5k
   v.grow_to(mod_words);
241
50.5k
   mp1o2.grow_to(mod_words);
242
50.5k
243
50.5k
   secure_vector<word>& a_w = a.get_word_vector();
244
50.5k
   secure_vector<word>& b_w = b.get_word_vector();
245
50.5k
   secure_vector<word>& u_w = u.get_word_vector();
246
50.5k
   secure_vector<word>& v_w = v.get_word_vector();
247
50.5k
248
50.5k
   CT::poison(a_w.data(), a_w.size());
249
50.5k
   CT::poison(b_w.data(), b_w.size());
250
50.5k
   CT::poison(u_w.data(), u_w.size());
251
50.5k
   CT::poison(v_w.data(), v_w.size());
252
50.5k
253
50.5k
   // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
254
50.5k
   size_t bits = 2 * mod.bits();
255
50.5k
256
42.5M
   while(bits--)
257
42.5M
      {
258
42.5M
      /*
259
42.5M
      const word odd = a.is_odd();
260
42.5M
      a -= odd * b;
261
42.5M
      const word underflow = a.is_negative();
262
42.5M
      b += a * underflow;
263
42.5M
      a.set_sign(BigInt::Positive);
264
42.5M
265
42.5M
      a >>= 1;
266
42.5M
267
42.5M
      if(underflow)
268
42.5M
         {
269
42.5M
         std::swap(u, v);
270
42.5M
         }
271
42.5M
272
42.5M
      u -= odd * v;
273
42.5M
      u += u.is_negative() * mod;
274
42.5M
275
42.5M
      const word odd_u = u.is_odd();
276
42.5M
277
42.5M
      u >>= 1;
278
42.5M
      u += mp1o2 * odd_u;
279
42.5M
      */
280
42.5M
281
42.5M
      const word odd_a = a_w[0] & 1;
282
42.5M
283
42.5M
      //if(odd_a) a -= b
284
42.5M
      word underflow = bigint_cnd_sub(odd_a, a_w.data(), b_w.data(), mod_words);
285
42.5M
286
42.5M
      //if(underflow) { b -= a; a = abs(a); swap(u, v); }
287
42.5M
      bigint_cnd_add(underflow, b_w.data(), a_w.data(), mod_words);
288
42.5M
      bigint_cnd_abs(underflow, a_w.data(), mod_words);
289
42.5M
      bigint_cnd_swap(underflow, u_w.data(), v_w.data(), mod_words);
290
42.5M
291
42.5M
      // a >>= 1
292
42.5M
      bigint_shr1(a_w.data(), mod_words, 0, 1);
293
42.5M
294
42.5M
      //if(odd_a) u -= v;
295
42.5M
      word borrow = bigint_cnd_sub(odd_a, u_w.data(), v_w.data(), mod_words);
296
42.5M
297
42.5M
      // if(borrow) u += p
298
42.5M
      bigint_cnd_add(borrow, u_w.data(), mod.data(), mod_words);
299
42.5M
300
42.5M
      const word odd_u = u_w[0] & 1;
301
42.5M
302
42.5M
      // u >>= 1
303
42.5M
      bigint_shr1(u_w.data(), mod_words, 0, 1);
304
42.5M
305
42.5M
      //if(odd_u) u += mp1o2;
306
42.5M
      bigint_cnd_add(odd_u, u_w.data(), mp1o2.data(), mod_words);
307
42.5M
      }
308
50.5k
309
50.5k
   CT::unpoison(a_w.data(), a_w.size());
310
50.5k
   CT::unpoison(b_w.data(), b_w.size());
311
50.5k
   CT::unpoison(u_w.data(), u_w.size());
312
50.5k
   CT::unpoison(v_w.data(), v_w.size());
313
50.5k
314
50.5k
   BOTAN_ASSERT(a.is_zero(), "A is zero");
315
50.5k
316
50.5k
   if(b != 1)
317
96
      return 0;
318
50.4k
319
50.4k
   return v;
320
50.4k
   }
321
322
/*
323
* Find the Modular Inverse
324
*/
325
BigInt inverse_mod(const BigInt& n, const BigInt& mod)
326
50.1k
   {
327
50.1k
   if(mod.is_zero())
328
0
      throw BigInt::DivideByZero();
329
50.1k
   if(mod.is_negative() || n.is_negative())
330
0
      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
331
50.1k
   if(n.is_zero())
332
144
      return 0;
333
49.9k
334
49.9k
   if(mod.is_odd() && n < mod)
335
49.8k
      return ct_inverse_mod_odd_modulus(n, mod);
336
148
337
148
   return inverse_euclid(n, mod);
338
148
   }
339
340
BigInt inverse_euclid(const BigInt& n, const BigInt& mod)
341
831
   {
342
831
   if(mod.is_zero())
343
0
      throw BigInt::DivideByZero();
344
831
   if(mod.is_negative() || n.is_negative())
345
0
      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
346
831
347
831
   if(n.is_zero() || (n.is_even() && mod.is_even()))
348
6
      return 0; // fast fail checks
349
825
350
825
   BigInt u = mod, v = n;
351
825
   BigInt A = 1, B = 0, C = 0, D = 1;
352
825
   BigInt T0, T1, T2;
353
825
354
342k
   while(u.is_nonzero())
355
341k
      {
356
341k
      const size_t u_zero_bits = low_zero_bits(u);
357
341k
      u >>= u_zero_bits;
358
341k
359
341k
      const size_t v_zero_bits = low_zero_bits(v);
360
341k
      v >>= v_zero_bits;
361
341k
362
341k
      const bool u_gte_v = (u >= v);
363
341k
364
754k
      for(size_t i = 0; i != u_zero_bits; ++i)
365
412k
         {
366
412k
         const bool needs_adjust = A.is_odd() || B.is_odd();
367
412k
368
412k
         T0 = A + n;
369
412k
         T1 = B - mod;
370
412k
371
412k
         A.ct_cond_assign(needs_adjust, T0);
372
412k
         B.ct_cond_assign(needs_adjust, T1);
373
412k
374
412k
         A >>= 1;
375
412k
         B >>= 1;
376
412k
         }
377
341k
378
604k
      for(size_t i = 0; i != v_zero_bits; ++i)
379
262k
         {
380
262k
         const bool needs_adjust = C.is_odd() || D.is_odd();
381
262k
         T0 = C + n;
382
262k
         T1 = D - mod;
383
262k
384
262k
         C.ct_cond_assign(needs_adjust, T0);
385
262k
         D.ct_cond_assign(needs_adjust, T1);
386
262k
387
262k
         C >>= 1;
388
262k
         D >>= 1;
389
262k
         }
390
341k
391
341k
      T0 = u - v;
392
341k
      T1 = A - C;
393
341k
      T2 = B - D;
394
341k
395
341k
      T0.cond_flip_sign(!u_gte_v);
396
341k
      T1.cond_flip_sign(!u_gte_v);
397
341k
      T2.cond_flip_sign(!u_gte_v);
398
341k
399
341k
      u.ct_cond_assign(u_gte_v, T0);
400
341k
      A.ct_cond_assign(u_gte_v, T1);
401
341k
      B.ct_cond_assign(u_gte_v, T2);
402
341k
403
341k
      v.ct_cond_assign(!u_gte_v, T0);
404
341k
      C.ct_cond_assign(!u_gte_v, T1);
405
341k
      D.ct_cond_assign(!u_gte_v, T2);
406
341k
      }
407
825
408
825
   if(v != 1)
409
113
      return 0; // no modular inverse
410
712
411
2.00k
   while(D.is_negative())
412
1.29k
      D += mod;
413
1.36k
   while(D >= mod)
414
652
      D -= mod;
415
712
416
712
   return D;
417
712
   }
418
419
word monty_inverse(word a)
420
88.5k
   {
421
88.5k
   if(a % 2 == 0)
422
0
      throw Invalid_Argument("monty_inverse only valid for odd integers");
423
88.5k
424
88.5k
   /*
425
88.5k
   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
426
88.5k
   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
427
88.5k
   */
428
88.5k
429
88.5k
   word b = 1;
430
88.5k
   word r = 0;
431
88.5k
432
5.75M
   for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
433
5.66M
      {
434
5.66M
      const word bi = b % 2;
435
5.66M
      r >>= 1;
436
5.66M
      r += bi << (BOTAN_MP_WORD_BITS - 1);
437
5.66M
438
5.66M
      b -= a * bi;
439
5.66M
      b >>= 1;
440
5.66M
      }
441
88.5k
442
88.5k
   // Now invert in addition space
443
88.5k
   r = (MP_WORD_MAX - r) + 1;
444
88.5k
445
88.5k
   return r;
446
88.5k
   }
447
448
/*
449
* Modular Exponentiation
450
*/
451
BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
452
74.0k
   {
453
74.0k
   if(mod.is_negative() || mod == 1)
454
22
      {
455
22
      return 0;
456
22
      }
457
74.0k
458
74.0k
   if(base.is_zero() || mod.is_zero())
459
19
      {
460
19
      if(exp.is_zero())
461
1
         return 1;
462
18
      return 0;
463
18
      }
464
74.0k
465
74.0k
   Modular_Reducer reduce_mod(mod);
466
74.0k
467
74.0k
   const size_t exp_bits = exp.bits();
468
74.0k
469
74.0k
   if(mod.is_odd())
470
73.4k
      {
471
73.4k
      const size_t powm_window = 4;
472
73.4k
473
73.4k
      auto monty_mod = std::make_shared<Montgomery_Params>(mod, reduce_mod);
474
73.4k
      auto powm_base_mod = monty_precompute(monty_mod, reduce_mod.reduce(base), powm_window);
475
73.4k
      return monty_execute(*powm_base_mod, exp, exp_bits);
476
73.4k
      }
477
629
478
629
   /*
479
629
   Support for even modulus is just a convenience and not considered
480
629
   cryptographically important, so this implementation is slow ...
481
629
   */
482
629
   BigInt accum = 1;
483
629
   BigInt g = reduce_mod.reduce(base);
484
629
   BigInt t;
485
629
486
282k
   for(size_t i = 0; i != exp_bits; ++i)
487
281k
      {
488
281k
      t = reduce_mod.multiply(g, accum);
489
281k
      g = reduce_mod.square(g);
490
281k
      accum.ct_cond_assign(exp.get_bit(i), t);
491
281k
      }
492
629
   return accum;
493
629
   }
494
495
496
BigInt is_perfect_square(const BigInt& C)
497
88
   {
498
88
   if(C < 1)
499
0
      throw Invalid_Argument("is_perfect_square requires C >= 1");
500
88
   if(C == 1)
501
0
      return 1;
502
88
503
88
   const size_t n = C.bits();
504
88
   const size_t m = (n + 1) / 2;
505
88
   const BigInt B = C + BigInt::power_of_2(m);
506
88
507
88
   BigInt X = BigInt::power_of_2(m) - 1;
508
88
   BigInt X2 = (X*X);
509
88
510
88
   for(;;)
511
453
      {
512
453
      X = (X2 + C) / (2*X);
513
453
      X2 = (X*X);
514
453
515
453
      if(X2 < B)
516
88
         break;
517
453
      }
518
88
519
88
   if(X2 == C)
520
0
      return X;
521
88
   else
522
88
      return 0;
523
88
   }
524
525
/*
526
* Test for primality using Miller-Rabin
527
*/
528
bool is_prime(const BigInt& n,
529
              RandomNumberGenerator& rng,
530
              size_t prob,
531
              bool is_random)
532
195
   {
533
195
   if(n == 2)
534
0
      return true;
535
195
   if(n <= 1 || n.is_even())
536
0
      return false;
537
195
538
195
   const size_t n_bits = n.bits();
539
195
540
195
   // Fast path testing for small numbers (<= 65521)
541
195
   if(n_bits <= 16)
542
0
      {
543
0
      const uint16_t num = static_cast<uint16_t>(n.word_at(0));
544
0
545
0
      return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
546
0
      }
547
195
548
195
   Modular_Reducer mod_n(n);
549
195
550
195
   if(rng.is_seeded())
551
195
      {
552
195
      const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
553
195
554
195
      if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false)
555
102
         return false;
556
93
557
93
      return is_lucas_probable_prime(n, mod_n);
558
93
      }
559
0
   else
560
0
      {
561
0
      return is_bailie_psw_probable_prime(n, mod_n);
562
0
      }
563
195
   }
564
565
}