Coverage Report

Created: 2020-11-21 08:34

/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
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
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
* Tonelli-Shanks algorithm
39
*/
40
BigInt ressol(const BigInt& a, const BigInt& p)
41
25.1k
   {
42
25.1k
   if(p <= 1 || p.is_even())
43
0
      throw Invalid_Argument("ressol: invalid prime");
44
45
25.1k
   if(a == 0)
46
1
      return 0;
47
25.1k
   else if(a < 0)
48
0
      throw Invalid_Argument("ressol: value to solve for must be positive");
49
25.1k
   else if(a >= p)
50
2
      throw Invalid_Argument("ressol: value to solve for must be less than p");
51
52
25.1k
   if(p == 2)
53
0
      return a;
54
55
25.1k
   if(jacobi(a, p) != 1) // not a quadratic residue
56
6.55k
      return -BigInt(1);
57
58
18.5k
   if(p % 4 == 3) // The easy case
59
17.1k
      {
60
17.1k
      return power_mod(a, ((p+1) >> 2), p);
61
17.1k
      }
62
63
1.39k
   size_t s = low_zero_bits(p - 1);
64
1.39k
   BigInt q = p >> s;
65
66
1.39k
   q -= 1;
67
1.39k
   q >>= 1;
68
69
1.39k
   Modular_Reducer mod_p(p);
70
71
1.39k
   BigInt r = power_mod(a, q, p);
72
1.39k
   BigInt n = mod_p.multiply(a, mod_p.square(r));
73
1.39k
   r = mod_p.multiply(r, a);
74
75
1.39k
   if(n == 1)
76
56
      return r;
77
78
   // find random quadratic nonresidue z
79
1.33k
   word z = 2;
80
1.33k
   for(;;)
81
12.5k
      {
82
12.5k
      if(jacobi(z, p) == -1) // found one
83
1.33k
         break;
84
85
11.1k
      z += 1; // try next z
86
87
      /*
88
      * The expected number of tests to find a non-residue modulo a
89
      * prime is 2. If we have not found one after 256 then almost
90
      * certainly we have been given a non-prime p.
91
      */
92
11.1k
      if(z >= 256)
93
0
         return -BigInt(1);
94
11.1k
      }
95
96
1.33k
   BigInt c = power_mod(z, (q << 1) + 1, p);
97
98
62.2k
   while(n > 1)
99
60.8k
      {
100
60.8k
      q = n;
101
102
60.8k
      size_t i = 0;
103
3.20M
      while(q != 1)
104
3.14M
         {
105
3.14M
         q = mod_p.square(q);
106
3.14M
         ++i;
107
108
3.14M
         if(i >= s)
109
0
            {
110
0
            return -BigInt(1);
111
0
            }
112
3.14M
         }
113
114
60.8k
      c = power_mod(c, BigInt::power_of_2(s-i-1), p);
115
60.8k
      r = mod_p.multiply(r, c);
116
60.8k
      c = mod_p.square(c);
117
60.8k
      n = mod_p.multiply(n, c);
118
60.8k
      s = i;
119
60.8k
      }
120
121
1.33k
   return r;
122
1.33k
   }
123
124
/*
125
* Calculate the Jacobi symbol
126
*/
127
int32_t jacobi(const BigInt& a, const BigInt& n)
128
40.1k
   {
129
40.1k
   if(n.is_even() || n < 2)
130
0
      throw Invalid_Argument("jacobi: second argument must be odd and > 1");
131
132
40.1k
   BigInt x = a % n;
133
40.1k
   BigInt y = n;
134
40.1k
   int32_t J = 1;
135
136
2.84M
   while(y > 1)
137
2.80M
      {
138
2.80M
      x %= y;
139
2.80M
      if(x > y / 2)
140
1.29M
         {
141
1.29M
         x = y - x;
142
1.29M
         if(y % 4 == 3)
143
634k
            J = -J;
144
1.29M
         }
145
2.80M
      if(x.is_zero())
146
0
         return 0;
147
148
2.80M
      size_t shifts = low_zero_bits(x);
149
2.80M
      x >>= shifts;
150
2.80M
      if(shifts % 2)
151
898k
         {
152
898k
         word y_mod_8 = y % 8;
153
898k
         if(y_mod_8 == 3 || y_mod_8 == 5)
154
444k
            J = -J;
155
898k
         }
156
157
2.80M
      if(x % 4 == 3 && y % 4 == 3)
158
678k
         J = -J;
159
2.80M
      std::swap(x, y);
160
2.80M
      }
161
40.1k
   return J;
162
40.1k
   }
163
164
/*
165
* Square a BigInt
166
*/
167
BigInt square(const BigInt& x)
168
4.96M
   {
169
4.96M
   BigInt z = x;
170
4.96M
   secure_vector<word> ws;
171
4.96M
   z.square(ws);
172
4.96M
   return z;
173
4.96M
   }
174
175
/*
176
* Return the number of 0 bits at the end of n
177
*/
178
size_t low_zero_bits(const BigInt& n)
179
3.65M
   {
180
3.65M
   size_t low_zero = 0;
181
182
3.65M
   auto seen_nonempty_word = CT::Mask<word>::cleared();
183
184
64.0M
   for(size_t i = 0; i != n.size(); ++i)
185
60.4M
      {
186
60.4M
      const word x = n.word_at(i);
187
188
      // ctz(0) will return sizeof(word)
189
60.4M
      const size_t tz_x = ctz(x);
190
191
      // if x > 0 we want to count tz_x in total but not any
192
      // further words, so set the mask after the addition
193
60.4M
      low_zero += seen_nonempty_word.if_not_set_return(tz_x);
194
195
60.4M
      seen_nonempty_word |= CT::Mask<word>::expand(x);
196
60.4M
      }
197
198
   // if we saw no words with x > 0 then n == 0 and the value we have
199
   // computed is meaningless. Instead return 0 in that case.
200
3.65M
   return seen_nonempty_word.if_set_return(low_zero);
201
3.65M
   }
202
203
/*
204
* Calculate the GCD
205
*/
206
BigInt gcd(const BigInt& a, const BigInt& b)
207
0
   {
208
0
   if(a.is_zero() || b.is_zero())
209
0
      return 0;
210
0
   if(a == 1 || b == 1)
211
0
      return 1;
212
213
   // See https://gcd.cr.yp.to/safegcd-20190413.pdf fig 1.2
214
215
0
   BigInt f = a;
216
0
   BigInt g = b;
217
0
   f.const_time_poison();
218
0
   g.const_time_poison();
219
220
0
   f.set_sign(BigInt::Positive);
221
0
   g.set_sign(BigInt::Positive);
222
223
0
   const size_t common2s = std::min(low_zero_bits(f), low_zero_bits(g));
224
0
   CT::unpoison(common2s);
225
226
0
   f >>= common2s;
227
0
   g >>= common2s;
228
229
0
   f.ct_cond_swap(f.is_even(), g);
230
231
0
   int32_t delta = 1;
232
233
0
   const size_t loop_cnt = 4 + 3*std::max(f.bits(), g.bits());
234
235
0
   BigInt newg, t;
236
0
   for(size_t i = 0; i != loop_cnt; ++i)
237
0
      {
238
0
      sub_abs(newg, f, g);
239
240
0
      const bool need_swap = (g.is_odd() && delta > 0);
241
242
      // if(need_swap) delta *= -1
243
0
      delta *= CT::Mask<uint8_t>::expand(need_swap).select(0, 2) - 1;
244
0
      f.ct_cond_swap(need_swap, g);
245
0
      g.ct_cond_swap(need_swap, newg);
246
247
0
      delta += 1;
248
249
0
      g.ct_cond_add(g.is_odd(), f);
250
0
      g >>= 1;
251
0
      }
252
253
0
   f <<= common2s;
254
255
0
   f.const_time_unpoison();
256
0
   g.const_time_unpoison();
257
258
0
   return f;
259
0
   }
260
261
/*
262
* Calculate the LCM
263
*/
264
BigInt lcm(const BigInt& a, const BigInt& b)
265
0
   {
266
0
   return ct_divide(a * b, gcd(a, b));
267
0
   }
268
269
/*
270
* Modular Exponentiation
271
*/
272
BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
273
81.9k
   {
274
81.9k
   if(mod.is_negative() || mod == 1)
275
29
      {
276
29
      return 0;
277
29
      }
278
279
81.9k
   if(base.is_zero() || mod.is_zero())
280
24
      {
281
24
      if(exp.is_zero())
282
4
         return 1;
283
20
      return 0;
284
20
      }
285
286
81.9k
   Modular_Reducer reduce_mod(mod);
287
288
81.9k
   const size_t exp_bits = exp.bits();
289
290
81.9k
   if(mod.is_odd())
291
81.2k
      {
292
81.2k
      const size_t powm_window = 4;
293
294
81.2k
      auto monty_mod = std::make_shared<Montgomery_Params>(mod, reduce_mod);
295
81.2k
      auto powm_base_mod = monty_precompute(monty_mod, reduce_mod.reduce(base), powm_window);
296
81.2k
      return monty_execute(*powm_base_mod, exp, exp_bits);
297
81.2k
      }
298
299
   /*
300
   Support for even modulus is just a convenience and not considered
301
   cryptographically important, so this implementation is slow ...
302
   */
303
673
   BigInt accum = 1;
304
673
   BigInt g = reduce_mod.reduce(base);
305
673
   BigInt t;
306
307
306k
   for(size_t i = 0; i != exp_bits; ++i)
308
305k
      {
309
305k
      t = reduce_mod.multiply(g, accum);
310
305k
      g = reduce_mod.square(g);
311
305k
      accum.ct_cond_assign(exp.get_bit(i), t);
312
305k
      }
313
673
   return accum;
314
673
   }
315
316
317
BigInt is_perfect_square(const BigInt& C)
318
89
   {
319
89
   if(C < 1)
320
0
      throw Invalid_Argument("is_perfect_square requires C >= 1");
321
89
   if(C == 1)
322
0
      return 1;
323
324
89
   const size_t n = C.bits();
325
89
   const size_t m = (n + 1) / 2;
326
89
   const BigInt B = C + BigInt::power_of_2(m);
327
328
89
   BigInt X = BigInt::power_of_2(m) - 1;
329
89
   BigInt X2 = (X*X);
330
331
89
   for(;;)
332
453
      {
333
453
      X = (X2 + C) / (2*X);
334
453
      X2 = (X*X);
335
336
453
      if(X2 < B)
337
89
         break;
338
453
      }
339
340
89
   if(X2 == C)
341
0
      return X;
342
89
   else
343
89
      return 0;
344
89
   }
345
346
/*
347
* Test for primality using Miller-Rabin
348
*/
349
bool is_prime(const BigInt& n,
350
              RandomNumberGenerator& rng,
351
              size_t prob,
352
              bool is_random)
353
210
   {
354
210
   if(n == 2)
355
0
      return true;
356
210
   if(n <= 1 || n.is_even())
357
0
      return false;
358
359
210
   const size_t n_bits = n.bits();
360
361
   // Fast path testing for small numbers (<= 65521)
362
210
   if(n_bits <= 16)
363
0
      {
364
0
      const uint16_t num = static_cast<uint16_t>(n.word_at(0));
365
366
0
      return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
367
0
      }
368
369
210
   Modular_Reducer mod_n(n);
370
371
210
   if(rng.is_seeded())
372
210
      {
373
210
      const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
374
375
210
      if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false)
376
103
         return false;
377
378
107
      if(is_random)
379
0
         return true;
380
107
      else
381
107
         return is_lucas_probable_prime(n, mod_n);
382
0
      }
383
0
   else
384
0
      {
385
0
      return is_bailie_psw_probable_prime(n, mod_n);
386
0
      }
387
210
   }
388
389
}