Line data Source code
1 : #include "tommath_private.h"
2 : #ifdef BN_MP_PRIME_IS_PRIME_C
3 : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
4 : /* SPDX-License-Identifier: Unlicense */
5 :
6 : /* portable integer log of two with small footprint */
7 0 : static unsigned int s_floor_ilog2(int value)
8 : {
9 0 : unsigned int r = 0;
10 0 : while ((value >>= 1) != 0) {
11 0 : r++;
12 : }
13 0 : return r;
14 : }
15 :
16 :
17 0 : mp_err mp_prime_is_prime(const mp_int *a, int t, mp_bool *result)
18 : {
19 0 : mp_int b;
20 0 : int ix, p_max = 0, size_a, len;
21 0 : mp_bool res;
22 0 : mp_err err;
23 0 : unsigned int fips_rand, mask;
24 :
25 : /* default to no */
26 0 : *result = MP_NO;
27 :
28 : /* Some shortcuts */
29 : /* N > 3 */
30 0 : if (a->used == 1) {
31 0 : if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
32 0 : *result = MP_NO;
33 0 : return MP_OKAY;
34 : }
35 0 : if (a->dp[0] == 2u) {
36 0 : *result = MP_YES;
37 0 : return MP_OKAY;
38 : }
39 : }
40 :
41 : /* N must be odd */
42 0 : if (MP_IS_EVEN(a)) {
43 0 : return MP_OKAY;
44 : }
45 : /* N is not a perfect square: floor(sqrt(N))^2 != N */
46 0 : if ((err = mp_is_square(a, &res)) != MP_OKAY) {
47 0 : return err;
48 : }
49 0 : if (res != MP_NO) {
50 0 : return MP_OKAY;
51 : }
52 :
53 : /* is the input equal to one of the primes in the table? */
54 0 : for (ix = 0; ix < PRIVATE_MP_PRIME_TAB_SIZE; ix++) {
55 0 : if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
56 0 : *result = MP_YES;
57 0 : return MP_OKAY;
58 : }
59 : }
60 : #ifdef MP_8BIT
61 : /* The search in the loop above was exhaustive in this case */
62 : if ((a->used == 1) && (PRIVATE_MP_PRIME_TAB_SIZE >= 31)) {
63 : return MP_OKAY;
64 : }
65 : #endif
66 :
67 : /* first perform trial division */
68 0 : if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
69 0 : return err;
70 : }
71 :
72 : /* return if it was trivially divisible */
73 0 : if (res == MP_YES) {
74 0 : return MP_OKAY;
75 : }
76 :
77 : /*
78 : Run the Miller-Rabin test with base 2 for the BPSW test.
79 : */
80 0 : if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
81 0 : return err;
82 : }
83 :
84 0 : if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
85 0 : goto LBL_B;
86 : }
87 0 : if (res == MP_NO) {
88 0 : goto LBL_B;
89 : }
90 : /*
91 : Rumours have it that Mathematica does a second M-R test with base 3.
92 : Other rumours have it that their strong L-S test is slightly different.
93 : It does not hurt, though, beside a bit of extra runtime.
94 : */
95 0 : b.dp[0]++;
96 0 : if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
97 0 : goto LBL_B;
98 : }
99 0 : if (res == MP_NO) {
100 0 : goto LBL_B;
101 : }
102 :
103 : /*
104 : * Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
105 : * slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
106 : * bases 2, 3 and t random bases.
107 : */
108 : #ifndef LTM_USE_ONLY_MR
109 0 : if (t >= 0) {
110 : /*
111 : * Use a Frobenius-Underwood test instead of the Lucas-Selfridge test for
112 : * MP_8BIT (It is unknown if the Lucas-Selfridge test works with 16-bit
113 : * integers but the necesssary analysis is on the todo-list).
114 : */
115 : #if defined (MP_8BIT) || defined (LTM_USE_FROBENIUS_TEST)
116 : err = mp_prime_frobenius_underwood(a, &res);
117 : if ((err != MP_OKAY) && (err != MP_ITER)) {
118 : goto LBL_B;
119 : }
120 : if (res == MP_NO) {
121 : goto LBL_B;
122 : }
123 : #else
124 0 : if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
125 0 : goto LBL_B;
126 : }
127 0 : if (res == MP_NO) {
128 0 : goto LBL_B;
129 : }
130 : #endif
131 : }
132 : #endif
133 :
134 : /* run at least one Miller-Rabin test with a random base */
135 0 : if (t == 0) {
136 0 : t = 1;
137 : }
138 :
139 : /*
140 : Only recommended if the input range is known to be < 3317044064679887385961981
141 :
142 : It uses the bases necessary for a deterministic M-R test if the input is
143 : smaller than 3317044064679887385961981
144 : The caller has to check the size.
145 : TODO: can be made a bit finer grained but comparing is not free.
146 : */
147 0 : if (t < 0) {
148 : /*
149 : Sorenson, Jonathan; Webster, Jonathan (2015).
150 : "Strong Pseudoprimes to Twelve Prime Bases".
151 : */
152 : /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
153 0 : if ((err = mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
154 0 : goto LBL_B;
155 : }
156 :
157 0 : if (mp_cmp(a, &b) == MP_LT) {
158 0 : p_max = 12;
159 : } else {
160 : /* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
161 0 : if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
162 0 : goto LBL_B;
163 : }
164 :
165 0 : if (mp_cmp(a, &b) == MP_LT) {
166 0 : p_max = 13;
167 : } else {
168 0 : err = MP_VAL;
169 0 : goto LBL_B;
170 : }
171 : }
172 :
173 : /* we did bases 2 and 3 already, skip them */
174 0 : for (ix = 2; ix < p_max; ix++) {
175 0 : mp_set(&b, s_mp_prime_tab[ix]);
176 0 : if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
177 0 : goto LBL_B;
178 : }
179 0 : if (res == MP_NO) {
180 0 : goto LBL_B;
181 : }
182 : }
183 : }
184 : /*
185 : Do "t" M-R tests with random bases between 3 and "a".
186 : See Fips 186.4 p. 126ff
187 : */
188 0 : else if (t > 0) {
189 : /*
190 : * The mp_digit's have a defined bit-size but the size of the
191 : * array a.dp is a simple 'int' and this library can not assume full
192 : * compliance to the current C-standard (ISO/IEC 9899:2011) because
193 : * it gets used for small embeded processors, too. Some of those MCUs
194 : * have compilers that one cannot call standard compliant by any means.
195 : * Hence the ugly type-fiddling in the following code.
196 : */
197 0 : size_a = mp_count_bits(a);
198 0 : mask = (1u << s_floor_ilog2(size_a)) - 1u;
199 : /*
200 : Assuming the General Rieman hypothesis (never thought to write that in a
201 : comment) the upper bound can be lowered to 2*(log a)^2.
202 : E. Bach, "Explicit bounds for primality testing and related problems,"
203 : Math. Comp. 55 (1990), 355-380.
204 :
205 : size_a = (size_a/10) * 7;
206 : len = 2 * (size_a * size_a);
207 :
208 : E.g.: a number of size 2^2048 would be reduced to the upper limit
209 :
210 : floor(2048/10)*7 = 1428
211 : 2 * 1428^2 = 4078368
212 :
213 : (would have been ~4030331.9962 with floats and natural log instead)
214 : That number is smaller than 2^28, the default bit-size of mp_digit.
215 : */
216 :
217 : /*
218 : How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
219 : does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
220 : Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.
221 :
222 : The function mp_rand() goes to some length to use a cryptographically
223 : good PRNG. That also means that the chance to always get the same base
224 : in the loop is non-zero, although very low.
225 : If the BPSW test and/or the addtional Frobenious test have been
226 : performed instead of just the Miller-Rabin test with the bases 2 and 3,
227 : a single extra test should suffice, so such a very unlikely event
228 : will not do much harm.
229 :
230 : To preemptivly answer the dangling question: no, a witness does not
231 : need to be prime.
232 : */
233 0 : for (ix = 0; ix < t; ix++) {
234 : /* mp_rand() guarantees the first digit to be non-zero */
235 0 : if ((err = mp_rand(&b, 1)) != MP_OKAY) {
236 0 : goto LBL_B;
237 : }
238 : /*
239 : * Reduce digit before casting because mp_digit might be bigger than
240 : * an unsigned int and "mask" on the other side is most probably not.
241 : */
242 0 : fips_rand = (unsigned int)(b.dp[0] & (mp_digit) mask);
243 : #ifdef MP_8BIT
244 : /*
245 : * One 8-bit digit is too small, so concatenate two if the size of
246 : * unsigned int allows for it.
247 : */
248 : if ((MP_SIZEOF_BITS(unsigned int)/2) >= MP_SIZEOF_BITS(mp_digit)) {
249 : if ((err = mp_rand(&b, 1)) != MP_OKAY) {
250 : goto LBL_B;
251 : }
252 : fips_rand <<= MP_SIZEOF_BITS(mp_digit);
253 : fips_rand |= (unsigned int) b.dp[0];
254 : fips_rand &= mask;
255 : }
256 : #endif
257 0 : if (fips_rand > (unsigned int)(INT_MAX - MP_DIGIT_BIT)) {
258 0 : len = INT_MAX / MP_DIGIT_BIT;
259 : } else {
260 0 : len = (((int)fips_rand + MP_DIGIT_BIT) / MP_DIGIT_BIT);
261 : }
262 : /* Unlikely. */
263 0 : if (len < 0) {
264 0 : ix--;
265 0 : continue;
266 : }
267 : /*
268 : * As mentioned above, one 8-bit digit is too small and
269 : * although it can only happen in the unlikely case that
270 : * an "unsigned int" is smaller than 16 bit a simple test
271 : * is cheap and the correction even cheaper.
272 : */
273 : #ifdef MP_8BIT
274 : /* All "a" < 2^8 have been caught before */
275 : if (len == 1) {
276 : len++;
277 : }
278 : #endif
279 0 : if ((err = mp_rand(&b, len)) != MP_OKAY) {
280 0 : goto LBL_B;
281 : }
282 : /*
283 : * That number might got too big and the witness has to be
284 : * smaller than "a"
285 : */
286 0 : len = mp_count_bits(&b);
287 0 : if (len >= size_a) {
288 0 : len = (len - size_a) + 1;
289 0 : if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
290 0 : goto LBL_B;
291 : }
292 : }
293 : /* Although the chance for b <= 3 is miniscule, try again. */
294 0 : if (mp_cmp_d(&b, 3uL) != MP_GT) {
295 0 : ix--;
296 0 : continue;
297 : }
298 0 : if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
299 0 : goto LBL_B;
300 : }
301 0 : if (res == MP_NO) {
302 0 : goto LBL_B;
303 : }
304 : }
305 : }
306 :
307 : /* passed the test */
308 0 : *result = MP_YES;
309 0 : LBL_B:
310 0 : mp_clear(&b);
311 0 : return err;
312 : }
313 :
314 : #endif
|
