Coverage Report

Created: 2021-04-07 06:07

/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
* (C) 2007,2008 Falko Strenzke, FlexSecure GmbH
5
*
6
* Botan is released under the Simplified BSD License (see license.txt)
7
*/
8
9
#include <botan/numthry.h>
10
#include <botan/reducer.h>
11
#include <botan/internal/monty.h>
12
#include <botan/internal/divide.h>
13
#include <botan/rng.h>
14
#include <botan/internal/ct_utils.h>
15
#include <botan/internal/mp_core.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
417k
   {
26
417k
   const size_t x_sw = x.sig_words();
27
417k
   const size_t y_sw = y.sig_words();
28
417k
   z.resize(std::max(x_sw, y_sw));
29
30
417k
   bigint_sub_abs(z.mutable_data(),
31
417k
                  x.data(), x_sw,
32
417k
                  y.data(), y_sw);
33
417k
   }
34
35
}
36
37
/*
38
* Tonelli-Shanks algorithm
39
*/
40
BigInt ressol(const BigInt& a, const BigInt& p)
41
23.6k
   {
42
23.6k
   if(p <= 1 || p.is_even())
43
0
      throw Invalid_Argument("ressol: invalid prime");
44
45
23.6k
   if(a == 0)
46
1
      return 0;
47
23.6k
   else if(a < 0)
48
0
      throw Invalid_Argument("ressol: value to solve for must be positive");
49
23.6k
   else if(a >= p)
50
2
      throw Invalid_Argument("ressol: value to solve for must be less than p");
51
52
23.6k
   if(p == 2)
53
0
      return a;
54
55
23.6k
   if(jacobi(a, p) != 1) // not a quadratic residue
56
9.22k
      return -BigInt(1);
57
58
14.3k
   Modular_Reducer mod_p(p);
59
14.3k
   auto monty_p = std::make_shared<Montgomery_Params>(p, mod_p);
60
61
14.3k
   if(p % 4 == 3) // The easy case
62
13.2k
      {
63
13.2k
      return monty_exp_vartime(monty_p, a, ((p+1) >> 2));
64
13.2k
      }
65
66
1.15k
   size_t s = low_zero_bits(p - 1);
67
1.15k
   BigInt q = p >> s;
68
69
1.15k
   q -= 1;
70
1.15k
   q >>= 1;
71
72
1.15k
   BigInt r = monty_exp_vartime(monty_p, a, q);
73
1.15k
   BigInt n = mod_p.multiply(a, mod_p.square(r));
74
1.15k
   r = mod_p.multiply(r, a);
75
76
1.15k
   if(n == 1)
77
43
      return r;
78
79
   // find random quadratic nonresidue z
80
1.11k
   word z = 2;
81
1.11k
   for(;;)
82
10.8k
      {
83
10.8k
      if(jacobi(z, p) == -1) // found one
84
1.11k
         break;
85
86
9.72k
      z += 1; // try next z
87
88
      /*
89
      * The expected number of tests to find a non-residue modulo a
90
      * prime is 2. If we have not found one after 256 then almost
91
      * certainly we have been given a non-prime p.
92
      */
93
9.72k
      if(z >= 256)
94
0
         return -BigInt(1);
95
9.72k
      }
96
97
1.11k
   BigInt c = monty_exp_vartime(monty_p, z, (q << 1) + 1);
98
99
51.8k
   while(n > 1)
100
50.7k
      {
101
50.7k
      q = n;
102
103
50.7k
      size_t i = 0;
104
2.70M
      while(q != 1)
105
2.65M
         {
106
2.65M
         q = mod_p.square(q);
107
2.65M
         ++i;
108
109
2.65M
         if(i >= s)
110
0
            {
111
0
            return -BigInt(1);
112
0
            }
113
2.65M
         }
114
115
50.7k
      BOTAN_ASSERT_NOMSG(s >= (i + 1)); // No underflow!
116
50.7k
      c = monty_exp_vartime(monty_p, c, BigInt::power_of_2(s-i-1));
117
50.7k
      r = mod_p.multiply(r, c);
118
50.7k
      c = mod_p.square(c);
119
50.7k
      n = mod_p.multiply(n, c);
120
121
      // s decreases as the algorithm proceeds
122
50.7k
      BOTAN_ASSERT_NOMSG(s >= i);
123
50.7k
      s = i;
124
50.7k
      }
125
126
1.11k
   return r;
127
1.11k
   }
128
129
/*
130
* Calculate the Jacobi symbol
131
*/
132
int32_t jacobi(const BigInt& a, const BigInt& n)
133
37.3k
   {
134
37.3k
   if(n.is_even() || n < 2)
135
0
      throw Invalid_Argument("jacobi: second argument must be odd and > 1");
136
137
37.3k
   BigInt x = a % n;
138
37.3k
   BigInt y = n;
139
37.3k
   int32_t J = 1;
140
141
2.77M
   while(y > 1)
142
2.73M
      {
143
2.73M
      x %= y;
144
2.73M
      if(x > y / 2)
145
1.26M
         {
146
1.26M
         x = y - x;
147
1.26M
         if(y % 4 == 3)
148
632k
            J = -J;
149
1.26M
         }
150
2.73M
      if(x.is_zero())
151
0
         return 0;
152
153
2.73M
      size_t shifts = low_zero_bits(x);
154
2.73M
      x >>= shifts;
155
2.73M
      if(shifts % 2)
156
880k
         {
157
880k
         word y_mod_8 = y % 8;
158
880k
         if(y_mod_8 == 3 || y_mod_8 == 5)
159
428k
            J = -J;
160
880k
         }
161
162
2.73M
      if(x % 4 == 3 && y % 4 == 3)
163
668k
         J = -J;
164
2.73M
      std::swap(x, y);
165
2.73M
      }
166
37.3k
   return J;
167
37.3k
   }
168
169
/*
170
* Square a BigInt
171
*/
172
BigInt square(const BigInt& x)
173
3.87M
   {
174
3.87M
   BigInt z = x;
175
3.87M
   secure_vector<word> ws;
176
3.87M
   z.square(ws);
177
3.87M
   return z;
178
3.87M
   }
179
180
/*
181
* Return the number of 0 bits at the end of n
182
*/
183
size_t low_zero_bits(const BigInt& n)
184
3.60M
   {
185
3.60M
   size_t low_zero = 0;
186
187
3.60M
   auto seen_nonempty_word = CT::Mask<word>::cleared();
188
189
64.3M
   for(size_t i = 0; i != n.size(); ++i)
190
60.7M
      {
191
60.7M
      const word x = n.word_at(i);
192
193
      // ctz(0) will return sizeof(word)
194
60.7M
      const size_t tz_x = ctz(x);
195
196
      // if x > 0 we want to count tz_x in total but not any
197
      // further words, so set the mask after the addition
198
60.7M
      low_zero += seen_nonempty_word.if_not_set_return(tz_x);
199
200
60.7M
      seen_nonempty_word |= CT::Mask<word>::expand(x);
201
60.7M
      }
202
203
   // if we saw no words with x > 0 then n == 0 and the value we have
204
   // computed is meaningless. Instead return 0 in that case.
205
3.60M
   return seen_nonempty_word.if_set_return(low_zero);
206
3.60M
   }
207
208
209
namespace {
210
211
size_t safegcd_loop_bound(size_t f_bits, size_t g_bits)
212
947
   {
213
947
   const size_t d = std::max(f_bits, g_bits);
214
215
947
   if(d < 46)
216
237
      return (49*d + 80) / 17;
217
710
   else
218
710
      return (49*d + 57) / 17;
219
947
   }
220
221
}
222
223
/*
224
* Calculate the GCD
225
*/
226
BigInt gcd(const BigInt& a, const BigInt& b)
227
960
   {
228
960
   if(a.is_zero())
229
9
      return abs(b);
230
951
   if(b.is_zero())
231
1
      return abs(a);
232
950
   if(a == 1 || b == 1)
233
3
      return 1;
234
235
   // See https://gcd.cr.yp.to/safegcd-20190413.pdf fig 1.2
236
237
947
   BigInt f = a;
238
947
   BigInt g = b;
239
947
   f.const_time_poison();
240
947
   g.const_time_poison();
241
242
947
   f.set_sign(BigInt::Positive);
243
947
   g.set_sign(BigInt::Positive);
244
245
947
   const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g));
246
947
   CT::unpoison(common2s);
247
248
947
   f >>= common2s;
249
947
   g >>= common2s;
250
251
947
   f.ct_cond_swap(f.is_even(), g);
252
253
947
   int32_t delta = 1;
254
255
947
   const size_t loop_cnt = safegcd_loop_bound(f.bits(), g.bits());
256
257
947
   BigInt newg, t;
258
418k
   for(size_t i = 0; i != loop_cnt; ++i)
259
417k
      {
260
417k
      sub_abs(newg, f, g);
261
262
417k
      const bool need_swap = (g.is_odd() && delta > 0);
263
264
      // if(need_swap) { delta *= -1 } else { delta *= 1 }
265
417k
      delta *= CT::Mask<uint8_t>::expand(need_swap).if_not_set_return(2) - 1;
266
417k
      f.ct_cond_swap(need_swap, g);
267
417k
      g.ct_cond_swap(need_swap, newg);
268
269
417k
      delta += 1;
270
271
417k
      g.ct_cond_add(g.is_odd(), f);
272
417k
      g >>= 1;
273
417k
      }
274
275
947
   f <<= common2s;
276
277
947
   f.const_time_unpoison();
278
947
   g.const_time_unpoison();
279
280
947
   BOTAN_ASSERT_NOMSG(g.is_zero());
281
282
947
   return f;
283
947
   }
284
285
/*
286
* Calculate the LCM
287
*/
288
BigInt lcm(const BigInt& a, const BigInt& b)
289
0
   {
290
0
   return ct_divide(a * b, gcd(a, b));
291
0
   }
292
293
/*
294
* Modular Exponentiation
295
*/
296
BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
297
1.03k
   {
298
1.03k
   if(mod.is_negative() || mod == 1)
299
22
      {
300
22
      return 0;
301
22
      }
302
303
1.01k
   if(base.is_zero() || mod.is_zero())
304
18
      {
305
18
      if(exp.is_zero())
306
2
         return 1;
307
16
      return 0;
308
16
      }
309
310
994
   Modular_Reducer reduce_mod(mod);
311
312
994
   const size_t exp_bits = exp.bits();
313
314
994
   if(mod.is_odd())
315
393
      {
316
393
      auto monty_params = std::make_shared<Montgomery_Params>(mod, reduce_mod);
317
393
      return monty_exp(monty_params, reduce_mod.reduce(base), exp, exp_bits);
318
393
      }
319
320
   /*
321
   Support for even modulus is just a convenience and not considered
322
   cryptographically important, so this implementation is slow ...
323
   */
324
601
   BigInt accum = 1;
325
601
   BigInt g = reduce_mod.reduce(base);
326
601
   BigInt t;
327
328
276k
   for(size_t i = 0; i != exp_bits; ++i)
329
275k
      {
330
275k
      t = reduce_mod.multiply(g, accum);
331
275k
      g = reduce_mod.square(g);
332
275k
      accum.ct_cond_assign(exp.get_bit(i), t);
333
275k
      }
334
601
   return accum;
335
601
   }
336
337
338
BigInt is_perfect_square(const BigInt& C)
339
90
   {
340
90
   if(C < 1)
341
0
      throw Invalid_Argument("is_perfect_square requires C >= 1");
342
90
   if(C == 1)
343
0
      return 1;
344
345
90
   const size_t n = C.bits();
346
90
   const size_t m = (n + 1) / 2;
347
90
   const BigInt B = C + BigInt::power_of_2(m);
348
349
90
   BigInt X = BigInt::power_of_2(m) - 1;
350
90
   BigInt X2 = (X*X);
351
352
90
   for(;;)
353
385
      {
354
385
      X = (X2 + C) / (2*X);
355
385
      X2 = (X*X);
356
357
385
      if(X2 < B)
358
90
         break;
359
385
      }
360
361
90
   if(X2 == C)
362
0
      return X;
363
90
   else
364
90
      return 0;
365
90
   }
366
367
/*
368
* Test for primality using Miller-Rabin
369
*/
370
bool is_prime(const BigInt& n,
371
              RandomNumberGenerator& rng,
372
              size_t prob,
373
              bool is_random)
374
0
   {
375
0
   if(n == 2)
376
0
      return true;
377
0
   if(n <= 1 || n.is_even())
378
0
      return false;
379
380
0
   const size_t n_bits = n.bits();
381
382
   // Fast path testing for small numbers (<= 65521)
383
0
   if(n_bits <= 16)
384
0
      {
385
0
      const uint16_t num = static_cast<uint16_t>(n.word_at(0));
386
387
0
      return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
388
0
      }
389
390
0
   Modular_Reducer mod_n(n);
391
392
0
   if(rng.is_seeded())
393
0
      {
394
0
      const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
395
396
0
      if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false)
397
0
         return false;
398
399
0
      if(is_random)
400
0
         return true;
401
0
      else
402
0
         return is_lucas_probable_prime(n, mod_n);
403
0
      }
404
0
   else
405
0
      {
406
0
      return is_bailie_psw_probable_prime(n, mod_n);
407
0
      }
408
0
   }
409
410
}