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