Coverage Report

Created: 2022-06-23 06:44

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