comparison libtommath/bn_mp_prime_is_prime.c @ 1655:f52919ffd3b1

update ltm to 1.1.0 and enable FIPS 186.4 compliant key-generation (#79) * make key-generation compliant to FIPS 186.4 * fix includes in tommath_class.h * update fuzzcorpus instead of error-out * fixup fuzzing make-targets * update Makefile.in * apply necessary patches to ltm sources * clean-up not required ltm files * update to vanilla ltm 1.1.0 this already only contains the required files * remove set/get double
author Steffen Jaeckel <s_jaeckel@gmx.de>
date Mon, 16 Sep 2019 15:50:38 +0200
parents 8bba51a55704
children a36e545fb43d
comparison
equal deleted inserted replaced
1654:cc0fc5131c5c 1655:f52919ffd3b1
1 #include <tommath_private.h> 1 #include "tommath_private.h"
2 #ifdef BN_MP_PRIME_IS_PRIME_C 2 #ifdef BN_MP_PRIME_IS_PRIME_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis 3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4 * 4 *
5 * LibTomMath is a library that provides multiple-precision 5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality. 6 * integer arithmetic as well as number theoretic functionality.
7 * 7 *
8 * The library was designed directly after the MPI library by 8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with 9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place. 10 * additional optimizations in place.
11 * 11 *
12 * The library is free for all purposes without any express 12 * SPDX-License-Identifier: Unlicense
13 * guarantee it works.
14 *
15 * Tom St Denis, [email protected], http://libtom.org
16 */ 13 */
17 14
18 /* performs a variable number of rounds of Miller-Rabin 15 /* portable integer log of two with small footprint */
19 * 16 static unsigned int s_floor_ilog2(int value)
20 * Probability of error after t rounds is no more than
21
22 *
23 * Sets result to 1 if probably prime, 0 otherwise
24 */
25 int mp_prime_is_prime (mp_int * a, int t, int *result)
26 { 17 {
27 mp_int b; 18 unsigned int r = 0;
28 int ix, err, res; 19 while ((value >>= 1) != 0) {
29 20 r++;
30 /* default to no */ 21 }
31 *result = MP_NO; 22 return r;
32 23 }
33 /* valid value of t? */ 24
34 if ((t <= 0) || (t > PRIME_SIZE)) { 25
35 return MP_VAL; 26 int mp_prime_is_prime(const mp_int *a, int t, int *result)
36 } 27 {
37 28 mp_int b;
38 /* is the input equal to one of the primes in the table? */ 29 int ix, err, res, p_max = 0, size_a, len;
39 for (ix = 0; ix < PRIME_SIZE; ix++) { 30 unsigned int fips_rand, mask;
40 if (mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) { 31
32 /* default to no */
33 *result = MP_NO;
34
35 /* valid value of t? */
36 if (t > PRIME_SIZE) {
37 return MP_VAL;
38 }
39
40 /* Some shortcuts */
41 /* N > 3 */
42 if (a->used == 1) {
43 if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
44 *result = 0;
45 return MP_OKAY;
46 }
47 if (a->dp[0] == 2u) {
41 *result = 1; 48 *result = 1;
42 return MP_OKAY; 49 return MP_OKAY;
43 } 50 }
44 } 51 }
45 52
46 /* first perform trial division */ 53 /* N must be odd */
47 if ((err = mp_prime_is_divisible (a, &res)) != MP_OKAY) { 54 if (mp_iseven(a) == MP_YES) {
48 return err; 55 return MP_OKAY;
49 } 56 }
50 57 /* N is not a perfect square: floor(sqrt(N))^2 != N */
51 /* return if it was trivially divisible */ 58 if ((err = mp_is_square(a, &res)) != MP_OKAY) {
52 if (res == MP_YES) { 59 return err;
53 return MP_OKAY; 60 }
54 } 61 if (res != 0) {
55 62 return MP_OKAY;
56 /* now perform the miller-rabin rounds */ 63 }
57 if ((err = mp_init (&b)) != MP_OKAY) { 64
58 return err; 65 /* is the input equal to one of the primes in the table? */
59 } 66 for (ix = 0; ix < PRIME_SIZE; ix++) {
60 67 if (mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) {
61 for (ix = 0; ix < t; ix++) { 68 *result = MP_YES;
62 /* set the prime */ 69 return MP_OKAY;
63 mp_set (&b, ltm_prime_tab[ix]); 70 }
64 71 }
65 if ((err = mp_prime_miller_rabin (a, &b, &res)) != MP_OKAY) { 72 #ifdef MP_8BIT
73 /* The search in the loop above was exhaustive in this case */
74 if ((a->used == 1) && (PRIME_SIZE >= 31)) {
75 return MP_OKAY;
76 }
77 #endif
78
79 /* first perform trial division */
80 if ((err = mp_prime_is_divisible(a, &res)) != MP_OKAY) {
81 return err;
82 }
83
84 /* return if it was trivially divisible */
85 if (res == MP_YES) {
86 return MP_OKAY;
87 }
88
89 /*
90 Run the Miller-Rabin test with base 2 for the BPSW test.
91 */
92 if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
93 return err;
94 }
95
96 if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
66 goto LBL_B; 97 goto LBL_B;
67 } 98 }
68 99 if (res == MP_NO) {
69 if (res == MP_NO) {
70 goto LBL_B; 100 goto LBL_B;
71 } 101 }
72 } 102 /*
73 103 Rumours have it that Mathematica does a second M-R test with base 3.
74 /* passed the test */ 104 Other rumours have it that their strong L-S test is slightly different.
75 *result = MP_YES; 105 It does not hurt, though, beside a bit of extra runtime.
76 LBL_B:mp_clear (&b); 106 */
77 return err; 107 b.dp[0]++;
108 if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
109 goto LBL_B;
110 }
111 if (res == MP_NO) {
112 goto LBL_B;
113 }
114
115 /*
116 * Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
117 * slow so if speed is an issue, define LTM_USE_FIPS_ONLY to use M-R tests with
118 * bases 2, 3 and t random bases.
119 */
120 #ifndef LTM_USE_FIPS_ONLY
121 if (t >= 0) {
122 /*
123 * Use a Frobenius-Underwood test instead of the Lucas-Selfridge test for
124 * MP_8BIT (It is unknown if the Lucas-Selfridge test works with 16-bit
125 * integers but the necesssary analysis is on the todo-list).
126 */
127 #if defined (MP_8BIT) || defined (LTM_USE_FROBENIUS_TEST)
128 err = mp_prime_frobenius_underwood(a, &res);
129 if ((err != MP_OKAY) && (err != MP_ITER)) {
130 goto LBL_B;
131 }
132 if (res == MP_NO) {
133 goto LBL_B;
134 }
135 #else
136 if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
137 goto LBL_B;
138 }
139 if (res == MP_NO) {
140 goto LBL_B;
141 }
142 #endif
143 }
144 #endif
145
146 /* run at least one Miller-Rabin test with a random base */
147 if (t == 0) {
148 t = 1;
149 }
150
151 /*
152 abs(t) extra rounds of M-R to extend the range of primes it can find if t < 0.
153 Only recommended if the input range is known to be < 3317044064679887385961981
154
155 It uses the bases for a deterministic M-R test if input < 3317044064679887385961981
156 The caller has to check the size.
157
158 Not for cryptographic use because with known bases strong M-R pseudoprimes can
159 be constructed. Use at least one M-R test with a random base (t >= 1).
160
161 The 1119 bit large number
162
163 80383745745363949125707961434194210813883768828755814583748891752229742737653\
164 33652186502336163960045457915042023603208766569966760987284043965408232928738\
165 79185086916685732826776177102938969773947016708230428687109997439976544144845\
166 34115587245063340927902227529622941498423068816854043264575340183297861112989\
167 60644845216191652872597534901
168
169 has been constructed by F. Arnault (F. Arnault, "Rabin-Miller primality test:
170 composite numbers which pass it.", Mathematics of Computation, 1995, 64. Jg.,
171 Nr. 209, S. 355-361), is a semiprime with the two factors
172
173 40095821663949960541830645208454685300518816604113250877450620473800321707011\
174 96242716223191597219733582163165085358166969145233813917169287527980445796800\
175 452592031836601
176
177 20047910831974980270915322604227342650259408302056625438725310236900160853505\
178 98121358111595798609866791081582542679083484572616906958584643763990222898400\
179 226296015918301
180
181 and it is a strong pseudoprime to all forty-six prime M-R bases up to 200
182
183 It does not fail the strong Bailley-PSP test as implemented here, it is just
184 given as an example, if not the reason to use the BPSW-test instead of M-R-tests
185 with a sequence of primes 2...n.
186
187 */
188 if (t < 0) {
189 t = -t;
190 /*
191 Sorenson, Jonathan; Webster, Jonathan (2015).
192 "Strong Pseudoprimes to Twelve Prime Bases".
193 */
194 /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
195 if ((err = mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
196 goto LBL_B;
197 }
198
199 if (mp_cmp(a, &b) == MP_LT) {
200 p_max = 12;
201 } else {
202 /* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
203 if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
204 goto LBL_B;
205 }
206
207 if (mp_cmp(a, &b) == MP_LT) {
208 p_max = 13;
209 } else {
210 err = MP_VAL;
211 goto LBL_B;
212 }
213 }
214
215 /* for compatibility with the current API (well, compatible within a sign's width) */
216 if (p_max < t) {
217 p_max = t;
218 }
219
220 if (p_max > PRIME_SIZE) {
221 err = MP_VAL;
222 goto LBL_B;
223 }
224 /* we did bases 2 and 3 already, skip them */
225 for (ix = 2; ix < p_max; ix++) {
226 mp_set(&b, ltm_prime_tab[ix]);
227 if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
228 goto LBL_B;
229 }
230 if (res == MP_NO) {
231 goto LBL_B;
232 }
233 }
234 }
235 /*
236 Do "t" M-R tests with random bases between 3 and "a".
237 See Fips 186.4 p. 126ff
238 */
239 else if (t > 0) {
240 /*
241 * The mp_digit's have a defined bit-size but the size of the
242 * array a.dp is a simple 'int' and this library can not assume full
243 * compliance to the current C-standard (ISO/IEC 9899:2011) because
244 * it gets used for small embeded processors, too. Some of those MCUs
245 * have compilers that one cannot call standard compliant by any means.
246 * Hence the ugly type-fiddling in the following code.
247 */
248 size_a = mp_count_bits(a);
249 mask = (1u << s_floor_ilog2(size_a)) - 1u;
250 /*
251 Assuming the General Rieman hypothesis (never thought to write that in a
252 comment) the upper bound can be lowered to 2*(log a)^2.
253 E. Bach, "Explicit bounds for primality testing and related problems,"
254 Math. Comp. 55 (1990), 355-380.
255
256 size_a = (size_a/10) * 7;
257 len = 2 * (size_a * size_a);
258
259 E.g.: a number of size 2^2048 would be reduced to the upper limit
260
261 floor(2048/10)*7 = 1428
262 2 * 1428^2 = 4078368
263
264 (would have been ~4030331.9962 with floats and natural log instead)
265 That number is smaller than 2^28, the default bit-size of mp_digit.
266 */
267
268 /*
269 How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
270 does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
271 Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.
272
273 The function mp_rand() goes to some length to use a cryptographically
274 good PRNG. That also means that the chance to always get the same base
275 in the loop is non-zero, although very low.
276 If the BPSW test and/or the addtional Frobenious test have been
277 performed instead of just the Miller-Rabin test with the bases 2 and 3,
278 a single extra test should suffice, so such a very unlikely event
279 will not do much harm.
280
281 To preemptivly answer the dangling question: no, a witness does not
282 need to be prime.
283 */
284 for (ix = 0; ix < t; ix++) {
285 /* mp_rand() guarantees the first digit to be non-zero */
286 if ((err = mp_rand(&b, 1)) != MP_OKAY) {
287 goto LBL_B;
288 }
289 /*
290 * Reduce digit before casting because mp_digit might be bigger than
291 * an unsigned int and "mask" on the other side is most probably not.
292 */
293 fips_rand = (unsigned int)(b.dp[0] & (mp_digit) mask);
294 #ifdef MP_8BIT
295 /*
296 * One 8-bit digit is too small, so concatenate two if the size of
297 * unsigned int allows for it.
298 */
299 if (((sizeof(unsigned int) * CHAR_BIT)/2) >= (sizeof(mp_digit) * CHAR_BIT)) {
300 if ((err = mp_rand(&b, 1)) != MP_OKAY) {
301 goto LBL_B;
302 }
303 fips_rand <<= sizeof(mp_digit) * CHAR_BIT;
304 fips_rand |= (unsigned int) b.dp[0];
305 fips_rand &= mask;
306 }
307 #endif
308 if (fips_rand > (unsigned int)(INT_MAX - DIGIT_BIT)) {
309 len = INT_MAX / DIGIT_BIT;
310 } else {
311 len = (((int)fips_rand + DIGIT_BIT) / DIGIT_BIT);
312 }
313 /* Unlikely. */
314 if (len < 0) {
315 ix--;
316 continue;
317 }
318 /*
319 * As mentioned above, one 8-bit digit is too small and
320 * although it can only happen in the unlikely case that
321 * an "unsigned int" is smaller than 16 bit a simple test
322 * is cheap and the correction even cheaper.
323 */
324 #ifdef MP_8BIT
325 /* All "a" < 2^8 have been caught before */
326 if (len == 1) {
327 len++;
328 }
329 #endif
330 if ((err = mp_rand(&b, len)) != MP_OKAY) {
331 goto LBL_B;
332 }
333 /*
334 * That number might got too big and the witness has to be
335 * smaller than or equal to "a"
336 */
337 len = mp_count_bits(&b);
338 if (len > size_a) {
339 len = len - size_a;
340 if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
341 goto LBL_B;
342 }
343 }
344
345 /* Although the chance for b <= 3 is miniscule, try again. */
346 if (mp_cmp_d(&b, 3uL) != MP_GT) {
347 ix--;
348 continue;
349 }
350 if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
351 goto LBL_B;
352 }
353 if (res == MP_NO) {
354 goto LBL_B;
355 }
356 }
357 }
358
359 /* passed the test */
360 *result = MP_YES;
361 LBL_B:
362 mp_clear(&b);
363 return err;
78 } 364 }
79 #endif 365
80 366 #endif
81 /* ref: $Format:%D$ */ 367
82 /* git commit: $Format:%H$ */ 368 /* ref: HEAD -> master, tag: v1.1.0 */
83 /* commit time: $Format:%ai$ */ 369 /* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */
370 /* commit time: 2019-01-28 20:32:32 +0100 */