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