Coverage Report

Created: 2021-05-04 09:02

/src/botan/src/lib/math/numbertheory/primality.cpp
Line
Count
Source (jump to first uncovered line)
1
/*
2
* (C) 2016,2018 Jack Lloyd
3
*
4
* Botan is released under the Simplified BSD License (see license.txt)
5
*/
6
7
#include <botan/internal/primality.h>
8
#include <botan/internal/monty_exp.h>
9
#include <botan/bigint.h>
10
#include <botan/internal/monty.h>
11
#include <botan/reducer.h>
12
#include <botan/rng.h>
13
#include <algorithm>
14
15
namespace Botan {
16
17
bool is_lucas_probable_prime(const BigInt& C, const Modular_Reducer& mod_C)
18
905
   {
19
905
   if(C <= 1)
20
0
      return false;
21
905
   else if(C == 2)
22
0
      return true;
23
905
   else if(C.is_even())
24
0
      return false;
25
905
   else if(C == 3 || C == 5 || C == 7 || C == 11 || C == 13)
26
7
      return true;
27
28
898
   BigInt D = BigInt::from_word(5);
29
30
898
   for(;;)
31
3.22k
      {
32
3.22k
      int32_t j = jacobi(D, C);
33
3.22k
      if(j == 0)
34
0
         return false;
35
36
3.22k
      if(j == -1)
37
898
         break;
38
39
      // Check 5, -7, 9, -11, 13, -15, 17, ...
40
2.32k
      if(D.is_negative())
41
933
         {
42
933
         D.flip_sign();
43
933
         D += 2;
44
933
         }
45
1.39k
      else
46
1.39k
         {
47
1.39k
         D += 2;
48
1.39k
         D.flip_sign();
49
1.39k
         }
50
51
2.32k
      if(D == 17 && is_perfect_square(C).is_nonzero())
52
0
         return false;
53
2.32k
      }
54
55
898
   const BigInt K = C + 1;
56
898
   const size_t K_bits = K.bits() - 1;
57
58
898
   BigInt U = BigInt::one();
59
898
   BigInt V = BigInt::one();
60
61
898
   BigInt Ut, Vt, U2, V2;
62
63
175k
   for(size_t i = 0; i != K_bits; ++i)
64
174k
      {
65
174k
      const bool k_bit = K.get_bit(K_bits - 1 - i);
66
67
174k
      Ut = mod_C.multiply(U, V);
68
69
174k
      Vt = mod_C.reduce(mod_C.square(V) + mod_C.multiply(D, mod_C.square(U)));
70
174k
      Vt.ct_cond_add(Vt.is_odd(), C);
71
174k
      Vt >>= 1;
72
174k
      Vt = mod_C.reduce(Vt);
73
74
174k
      U = Ut;
75
174k
      V = Vt;
76
77
174k
      U2 = mod_C.reduce(Ut + Vt);
78
174k
      U2.ct_cond_add(U2.is_odd(), C);
79
174k
      U2 >>= 1;
80
81
174k
      V2 = mod_C.reduce(Vt + Ut*D);
82
174k
      V2.ct_cond_add(V2.is_odd(), C);
83
174k
      V2 >>= 1;
84
85
174k
      U.ct_cond_assign(k_bit, U2);
86
174k
      V.ct_cond_assign(k_bit, V2);
87
174k
      }
88
89
898
   return (U == 0);
90
898
   }
91
92
bool is_bailie_psw_probable_prime(const BigInt& n, const Modular_Reducer& mod_n)
93
986
   {
94
986
   if(n < 3 || n.is_even())
95
6
      return false;
96
97
980
   auto monty_n = std::make_shared<Montgomery_Params>(n, mod_n);
98
980
   const auto base = BigInt::from_word(2);
99
980
   return passes_miller_rabin_test(n, mod_n, monty_n, base) && is_lucas_probable_prime(n, mod_n);
100
980
   }
101
102
bool is_bailie_psw_probable_prime(const BigInt& n)
103
986
   {
104
986
   Modular_Reducer mod_n(n);
105
986
   return is_bailie_psw_probable_prime(n, mod_n);
106
986
   }
107
108
bool passes_miller_rabin_test(const BigInt& n,
109
                              const Modular_Reducer& mod_n,
110
                              const std::shared_ptr<Montgomery_Params>& monty_n,
111
                              const BigInt& a)
112
1.02k
   {
113
1.02k
   if(n < 3 || n.is_even())
114
0
      return false;
115
116
1.02k
   BOTAN_ASSERT_NOMSG(n > 1);
117
118
1.02k
   const BigInt n_minus_1 = n - 1;
119
1.02k
   const size_t s = low_zero_bits(n_minus_1);
120
1.02k
   const BigInt nm1_s = n_minus_1 >> s;
121
1.02k
   const size_t n_bits = n.bits();
122
123
1.02k
   const size_t powm_window = 4;
124
125
1.02k
   auto powm_a_n = monty_precompute(monty_n, a, powm_window);
126
127
1.02k
   BigInt y = monty_execute(*powm_a_n, nm1_s, n_bits);
128
129
1.02k
   if(y == 1 || y == n_minus_1)
130
392
      return true;
131
132
21.0k
   for(size_t i = 1; i != s; ++i)
133
20.9k
      {
134
20.9k
      y = mod_n.square(y);
135
136
20.9k
      if(y == 1) // found a non-trivial square root
137
1
         return false;
138
139
      /*
140
      -1 is the trivial square root of unity, so ``a`` is not a
141
      witness for this number - give up
142
      */
143
20.9k
      if(y == n_minus_1)
144
541
         return true;
145
20.9k
      }
146
147
89
   return false;
148
631
   }
149
150
bool is_miller_rabin_probable_prime(const BigInt& n,
151
                                    const Modular_Reducer& mod_n,
152
                                    RandomNumberGenerator& rng,
153
                                    size_t test_iterations)
154
15
   {
155
15
   if(n < 3 || n.is_even())
156
0
      return false;
157
158
15
   auto monty_n = std::make_shared<Montgomery_Params>(n, mod_n);
159
160
44
   for(size_t i = 0; i != test_iterations; ++i)
161
43
      {
162
43
      const BigInt a = BigInt::random_integer(rng, BigInt::from_word(2), n);
163
164
43
      if(!passes_miller_rabin_test(n, mod_n, monty_n, a))
165
14
         return false;
166
43
      }
167
168
   // Failed to find a counterexample
169
1
   return true;
170
15
   }
171
172
173
size_t miller_rabin_test_iterations(size_t n_bits, size_t prob, bool random)
174
1
   {
175
1
   const size_t base = (prob + 2) / 2; // worst case 4^-t error rate
176
177
   /*
178
   * If the candidate prime was maliciously constructed, we can't rely
179
   * on arguments based on p being random.
180
   */
181
1
   if(random == false)
182
0
      return base;
183
184
   /*
185
   * For randomly chosen numbers we can use the estimates from
186
   * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf
187
   *
188
   * These values are derived from the inequality for p(k,t) given on
189
   * the second page.
190
   */
191
1
   if(prob <= 128)
192
1
      {
193
1
      if(n_bits >= 1536)
194
0
         return 4; // < 2^-133
195
1
      if(n_bits >= 1024)
196
0
         return 6; // < 2^-133
197
1
      if(n_bits >= 512)
198
0
         return 12; // < 2^-129
199
1
      if(n_bits >= 256)
200
1
         return 29; // < 2^-128
201
0
      }
202
203
   /*
204
   If the user desires a smaller error probability than we have
205
   precomputed error estimates for, just fall back to using the worst
206
   case error rate.
207
   */
208
0
   return base;
209
0
   }
210
211
}