Coverage Report

Created: 2021-02-21 07:20

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