xref: /dragonfly/crypto/libressl/crypto/bn/bn_bpsw.c (revision 6f5ec8b5) 
1 /*	$OpenBSD: bn_bpsw.c,v 1.7 2022/08/31 21:34:14 tb Exp $ */
2 /*
3  * Copyright (c) 2022 Martin Grenouilloux <martin.grenouilloux@lse.epita.fr>
4  * Copyright (c) 2022 Theo Buehler <tb@openbsd.org>
5  *
6  * Permission to use, copy, modify, and distribute this software for any
7  * purpose with or without fee is hereby granted, provided that the above
8  * copyright notice and this permission notice appear in all copies.
9  *
10  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17  */
18 
19 #include <openssl/bn.h>
20 
21 #include "bn_lcl.h"
22 #include "bn_prime.h"
23 
24 /*
25  * For an odd n compute a / 2 (mod n). If a is even, we can do a plain
26  * division, otherwise calculate (a + n) / 2. Then reduce (mod n).
27  */
28 
29 static int
30 bn_div_by_two_mod_odd_n(BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
31 {
32 	if (!BN_is_odd(n))
33 		return 0;
34 
35 	if (BN_is_odd(a)) {
36 		if (!BN_add(a, a, n))
37 			return 0;
38 	}
39 	if (!BN_rshift1(a, a))
40 		return 0;
41 	if (!BN_mod_ct(a, a, n, ctx))
42 		return 0;
43 
44 	return 1;
45 }
46 
47 /*
48  * Given the next binary digit of k and the current Lucas terms U and V, this
49  * helper computes the next terms in the Lucas sequence defined as follows:
50  *
51  *   U' = U * V                  (mod n)
52  *   V' = (V^2 + D * U^2) / 2    (mod n)
53  *
54  * If digit == 0, bn_lucas_step() returns U' and V'. If digit == 1, it returns
55  *
56  *   U'' = (U' + V') / 2         (mod n)
57  *   V'' = (V' + D * U') / 2     (mod n)
58  *
59  * Compare with FIPS 186-4, Appendix C.3.3, step 6.
60  */
61 
62 static int
63 bn_lucas_step(BIGNUM *U, BIGNUM *V, int digit, const BIGNUM *D,
64     const BIGNUM *n, BN_CTX *ctx)
65 {
66 	BIGNUM *tmp;
67 	int ret = 0;
68 
69 	BN_CTX_start(ctx);
70 
71 	if ((tmp = BN_CTX_get(ctx)) == NULL)
72 		goto err;
73 
74 	/* Calculate D * U^2 before computing U'. */
75 	if (!BN_sqr(tmp, U, ctx))
76 		goto err;
77 	if (!BN_mul(tmp, D, tmp, ctx))
78 		goto err;
79 
80 	/* U' = U * V (mod n). */
81 	if (!BN_mod_mul(U, U, V, n, ctx))
82 		goto err;
83 
84 	/* V' = (V^2 + D * U^2) / 2 (mod n). */
85 	if (!BN_sqr(V, V, ctx))
86 		goto err;
87 	if (!BN_add(V, V, tmp))
88 		goto err;
89 	if (!bn_div_by_two_mod_odd_n(V, n, ctx))
90 		goto err;
91 
92 	if (digit == 1) {
93 		/* Calculate D * U' before computing U''. */
94 		if (!BN_mul(tmp, D, U, ctx))
95 			goto err;
96 
97 		/* U'' = (U' + V') / 2 (mod n). */
98 		if (!BN_add(U, U, V))
99 			goto err;
100 		if (!bn_div_by_two_mod_odd_n(U, n, ctx))
101 			goto err;
102 
103 		/* V'' = (V' + D * U') / 2 (mod n). */
104 		if (!BN_add(V, V, tmp))
105 			goto err;
106 		if (!bn_div_by_two_mod_odd_n(V, n, ctx))
107 			goto err;
108 	}
109 
110 	ret = 1;
111 
112  err:
113 	BN_CTX_end(ctx);
114 
115 	return ret;
116 }
117 
118 /*
119  * Compute the Lucas terms U_k, V_k, see FIPS 186-4, Appendix C.3.3, steps 4-6.
120  */
121 
122 static int
123 bn_lucas(BIGNUM *U, BIGNUM *V, const BIGNUM *k, const BIGNUM *D,
124     const BIGNUM *n, BN_CTX *ctx)
125 {
126 	int digit, i;
127 	int ret = 0;
128 
129 	if (!BN_one(U))
130 		goto err;
131 	if (!BN_one(V))
132 		goto err;
133 
134 	/*
135 	 * Iterate over the digits of k from MSB to LSB. Start at digit 2
136 	 * since the first digit is dealt with by setting U = 1 and V = 1.
137 	 */
138 
139 	for (i = BN_num_bits(k) - 2; i >= 0; i--) {
140 		digit = BN_is_bit_set(k, i);
141 
142 		if (!bn_lucas_step(U, V, digit, D, n, ctx))
143 			goto err;
144 	}
145 
146 	ret = 1;
147 
148  err:
149 	return ret;
150 }
151 
152 /*
153  * This is a stronger variant of the Lucas test in FIPS 186-4, Appendix C.3.3.
154  * Every strong Lucas pseudoprime n is also a Lucas pseudoprime since
155  * U_{n+1} == 0 follows from U_k == 0 or V_{k * 2^r} == 0 for 0 <= r < s.
156  */
157 
158 static int
159 bn_strong_lucas_test(int *is_prime, const BIGNUM *n, const BIGNUM *D,
160     BN_CTX *ctx)
161 {
162 	BIGNUM *k, *U, *V;
163 	int r, s;
164 	int ret = 0;
165 
166 	BN_CTX_start(ctx);
167 
168 	if ((k = BN_CTX_get(ctx)) == NULL)
169 		goto err;
170 	if ((U = BN_CTX_get(ctx)) == NULL)
171 		goto err;
172 	if ((V = BN_CTX_get(ctx)) == NULL)
173 		goto err;
174 
175 	/*
176 	 * Factorize n + 1 = k * 2^s with odd k: shift away the s trailing ones
177 	 * of n and set the lowest bit of the resulting number k.
178 	 */
179 
180 	s = 0;
181 	while (BN_is_bit_set(n, s))
182 		s++;
183 	if (!BN_rshift(k, n, s))
184 		goto err;
185 	if (!BN_set_bit(k, 0))
186 		goto err;
187 
188 	/*
189 	 * Calculate the Lucas terms U_k and V_k. If either of them is zero,
190 	 * then n is a strong Lucas pseudoprime.
191 	 */
192 
193 	if (!bn_lucas(U, V, k, D, n, ctx))
194 		goto err;
195 
196 	if (BN_is_zero(U) || BN_is_zero(V)) {
197 		*is_prime = 1;
198 		goto done;
199 	}
200 
201 	/*
202 	 * Calculate the Lucas terms U_{k * 2^r}, V_{k * 2^r} for 1 <= r < s.
203 	 * If any V_{k * 2^r} is zero then n is a strong Lucas pseudoprime.
204 	 */
205 
206 	for (r = 1; r < s; r++) {
207 		if (!bn_lucas_step(U, V, 0, D, n, ctx))
208 			goto err;
209 
210 		if (BN_is_zero(V)) {
211 			*is_prime = 1;
212 			goto done;
213 		}
214 	}
215 
216 	/*
217 	 * If we got here, n is definitely composite.
218 	 */
219 
220 	*is_prime = 0;
221 
222  done:
223 	ret = 1;
224 
225  err:
226 	BN_CTX_end(ctx);
227 
228 	return ret;
229 }
230 
231 /*
232  * Test n for primality using the strong Lucas test with Selfridge's Method A.
233  * Returns 1 if n is prime or a strong Lucas-Selfridge pseudoprime.
234  * If it returns 0 then n is definitely composite.
235  */
236 
237 static int
238 bn_strong_lucas_selfridge(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
239 {
240 	BIGNUM *D, *two;
241 	int is_perfect_square, jacobi_symbol, sign;
242 	int ret = 0;
243 
244 	BN_CTX_start(ctx);
245 
246 	/* If n is a perfect square, it is composite. */
247 	if (!bn_is_perfect_square(&is_perfect_square, n, ctx))
248 		goto err;
249 	if (is_perfect_square) {
250 		*is_prime = 0;
251 		goto done;
252 	}
253 
254 	/*
255 	 * Find the first D in the Selfridge sequence 5, -7, 9, -11, 13, ...
256 	 * such that the Jacobi symbol (D/n) is -1.
257 	 */
258 
259 	if ((D = BN_CTX_get(ctx)) == NULL)
260 		goto err;
261 	if ((two = BN_CTX_get(ctx)) == NULL)
262 		goto err;
263 
264 	sign = 1;
265 	if (!BN_set_word(D, 5))
266 		goto err;
267 	if (!BN_set_word(two, 2))
268 		goto err;
269 
270 	while (1) {
271 		/* For odd n the Kronecker symbol computes the Jacobi symbol. */
272 		if ((jacobi_symbol = BN_kronecker(D, n, ctx)) == -2)
273 			goto err;
274 
275 		/* We found the value for D. */
276 		if (jacobi_symbol == -1)
277 			break;
278 
279 		/* n and D have prime factors in common. */
280 		if (jacobi_symbol == 0) {
281 			*is_prime = 0;
282 			goto done;
283 		}
284 
285 		sign = -sign;
286 		if (!BN_uadd(D, D, two))
287 			goto err;
288 		BN_set_negative(D, sign == -1);
289 	}
290 
291 	if (!bn_strong_lucas_test(is_prime, n, D, ctx))
292 		goto err;
293 
294  done:
295 	ret = 1;
296 
297  err:
298 	BN_CTX_end(ctx);
299 
300 	return ret;
301 }
302 
303 /*
304  * Miller-Rabin primality test for base 2.
305  */
306 
307 static int
308 bn_miller_rabin_base_2(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
309 {
310 	BIGNUM *n_minus_one, *k, *x;
311 	int i, s;
312 	int ret = 0;
313 
314 	BN_CTX_start(ctx);
315 
316 	if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
317 		goto err;
318 	if ((k = BN_CTX_get(ctx)) == NULL)
319 		goto err;
320 	if ((x = BN_CTX_get(ctx)) == NULL)
321 		goto err;
322 
323 	if (BN_is_word(n, 2) || BN_is_word(n, 3)) {
324 		*is_prime = 1;
325 		goto done;
326 	}
327 
328 	if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
329 		*is_prime = 0;
330 		goto done;
331 	}
332 
333 	if (!BN_sub(n_minus_one, n, BN_value_one()))
334 		goto err;
335 
336 	/*
337 	 * Factorize n - 1 = k * 2^s.
338 	 */
339 
340 	s = 0;
341 	while (!BN_is_bit_set(n_minus_one, s))
342 		s++;
343 	if (!BN_rshift(k, n_minus_one, s))
344 		goto err;
345 
346 	/*
347 	 * If 2^k is 1 or -1 (mod n) then n is a 2-pseudoprime.
348 	 */
349 
350 	if (!BN_set_word(x, 2))
351 		goto err;
352 	if (!BN_mod_exp_ct(x, x, k, n, ctx))
353 		goto err;
354 
355 	if (BN_is_one(x) || BN_cmp(x, n_minus_one) == 0) {
356 		*is_prime = 1;
357 		goto done;
358 	}
359 
360 	/*
361 	 * If 2^{2^i k} == -1 (mod n) for some 1 <= i < s, then n is a
362 	 * 2-pseudoprime.
363 	 */
364 
365 	for (i = 1; i < s; i++) {
366 		if (!BN_mod_sqr(x, x, n, ctx))
367 			goto err;
368 		if (BN_cmp(x, n_minus_one) == 0) {
369 			*is_prime = 1;
370 			goto done;
371 		}
372 	}
373 
374 	/*
375 	 * If we got here, n is definitely composite.
376 	 */
377 
378 	*is_prime = 0;
379 
380  done:
381 	ret = 1;
382 
383  err:
384 	BN_CTX_end(ctx);
385 
386 	return ret;
387 }
388 
389 /*
390  * The Baillie-Pomerance-Selfridge-Wagstaff algorithm combines a Miller-Rabin
391  * test for base 2 with a Strong Lucas pseudoprime test.
392  */
393 
394 int
395 bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx)
396 {
397 	BN_CTX *ctx = NULL;
398 	BN_ULONG mod;
399 	int i;
400 	int ret = 0;
401 
402 	if (BN_is_word(n, 2)) {
403 		*is_prime = 1;
404 		goto done;
405 	}
406 
407 	if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
408 		*is_prime = 0;
409 		goto done;
410 	}
411 
412 	/* Trial divisions with the first 2048 primes. */
413 	for (i = 0; i < NUMPRIMES; i++) {
414 		if ((mod = BN_mod_word(n, primes[i])) == (BN_ULONG)-1)
415 			goto err;
416 		if (mod == 0) {
417 			*is_prime = BN_is_word(n, primes[i]);
418 			goto done;
419 		}
420 	}
421 
422 	if ((ctx = in_ctx) == NULL)
423 		ctx = BN_CTX_new();
424 	if (ctx == NULL)
425 		goto err;
426 
427 	if (!bn_miller_rabin_base_2(is_prime, n, ctx))
428 		goto err;
429 	if (!*is_prime)
430 		goto done;
431 
432 	/* XXX - Miller-Rabin for random bases? See FIPS 186-4, Table C.1. */
433 
434 	if (!bn_strong_lucas_selfridge(is_prime, n, ctx))
435 		goto err;
436 
437  done:
438 	ret = 1;
439 
440  err:
441 	if (ctx != in_ctx)
442 		BN_CTX_free(ctx);
443 
444 	return ret;
445 }
446