1 /* $OpenBSD: bn_bpsw.c,v 1.11 2023/08/03 18:53:55 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_local.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
bn_div_by_two_mod_odd_n(BIGNUM * a,const BIGNUM * n,BN_CTX * ctx)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
bn_lucas_step(BIGNUM * U,BIGNUM * V,int digit,const BIGNUM * D,const BIGNUM * n,BN_CTX * ctx)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
bn_lucas(BIGNUM * U,BIGNUM * V,const BIGNUM * k,const BIGNUM * D,const BIGNUM * n,BN_CTX * ctx)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
bn_strong_lucas_test(int * is_pseudoprime,const BIGNUM * n,const BIGNUM * D,BN_CTX * ctx)159 bn_strong_lucas_test(int *is_pseudoprime, 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_pseudoprime = 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_pseudoprime = 1;
212 goto done;
213 }
214 }
215
216 /*
217 * If we got here, n is definitely composite.
218 */
219
220 *is_pseudoprime = 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
bn_strong_lucas_selfridge(int * is_pseudoprime,const BIGNUM * n,BN_CTX * ctx)238 bn_strong_lucas_selfridge(int *is_pseudoprime, 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_pseudoprime = 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_pseudoprime = 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_pseudoprime, 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 * Fermat criterion in Miller-Rabin test.
305 *
306 * Check whether 1 < base < n - 1 witnesses that n is composite. For prime n:
307 *
308 * * Fermat's little theorem: base^(n-1) = 1 (mod n).
309 * * The only square roots of 1 (mod n) are 1 and -1.
310 *
311 * Calculate base^((n-1)/2) by writing n - 1 = k * 2^s with odd k. Iteratively
312 * compute power = (base^k)^(2^(s-1)) by successive squaring of base^k.
313 *
314 * If power ever reaches -1, base^(n-1) is equal to 1 and n is a pseudoprime
315 * for base. If power reaches 1 before -1 during successive squaring, we have
316 * an unexpected square root of 1 and n is composite. Otherwise base^(n-1) != 1,
317 * and n is composite.
318 */
319
320 static int
bn_fermat(int * is_pseudoprime,const BIGNUM * n,const BIGNUM * n_minus_one,const BIGNUM * k,int s,const BIGNUM * base,BN_CTX * ctx,BN_MONT_CTX * mctx)321 bn_fermat(int *is_pseudoprime, const BIGNUM *n, const BIGNUM *n_minus_one,
322 const BIGNUM *k, int s, const BIGNUM *base, BN_CTX *ctx, BN_MONT_CTX *mctx)
323 {
324 BIGNUM *power;
325 int ret = 0;
326 int i;
327
328 BN_CTX_start(ctx);
329
330 if ((power = BN_CTX_get(ctx)) == NULL)
331 goto err;
332
333 /* Sanity check: ensure that 1 < base < n - 1. */
334 if (BN_cmp(base, BN_value_one()) <= 0 || BN_cmp(base, n_minus_one) >= 0)
335 goto err;
336
337 if (!BN_mod_exp_mont_ct(power, base, k, n, ctx, mctx))
338 goto err;
339
340 if (BN_is_one(power) || BN_cmp(power, n_minus_one) == 0) {
341 *is_pseudoprime = 1;
342 goto done;
343 }
344
345 /* Loop invariant: power is neither 1 nor -1 (mod n). */
346 for (i = 1; i < s; i++) {
347 if (!BN_mod_sqr(power, power, n, ctx))
348 goto err;
349
350 /* n is a pseudoprime for base. */
351 if (BN_cmp(power, n_minus_one) == 0) {
352 *is_pseudoprime = 1;
353 goto done;
354 }
355
356 /* n is composite: there's a square root of unity != 1 or -1. */
357 if (BN_is_one(power)) {
358 *is_pseudoprime = 0;
359 goto done;
360 }
361 }
362
363 /*
364 * If we get here, n is definitely composite: base^(n-1) != 1.
365 */
366
367 *is_pseudoprime = 0;
368
369 done:
370 ret = 1;
371
372 err:
373 BN_CTX_end(ctx);
374
375 return ret;
376 }
377
378 /*
379 * Miller-Rabin primality test for base 2 and for |rounds| of random bases.
380 * On success: is_pseudoprime == 0 implies that n is composite.
381 */
382
383 static int
bn_miller_rabin(int * is_pseudoprime,const BIGNUM * n,BN_CTX * ctx,size_t rounds)384 bn_miller_rabin(int *is_pseudoprime, const BIGNUM *n, BN_CTX *ctx,
385 size_t rounds)
386 {
387 BN_MONT_CTX *mctx = NULL;
388 BIGNUM *base, *k, *n_minus_one;
389 size_t i;
390 int s;
391 int ret = 0;
392
393 BN_CTX_start(ctx);
394
395 if ((base = BN_CTX_get(ctx)) == NULL)
396 goto err;
397 if ((k = BN_CTX_get(ctx)) == NULL)
398 goto err;
399 if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
400 goto err;
401
402 if (BN_is_word(n, 2) || BN_is_word(n, 3)) {
403 *is_pseudoprime = 1;
404 goto done;
405 }
406
407 if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
408 *is_pseudoprime = 0;
409 goto done;
410 }
411
412 if (!BN_sub(n_minus_one, n, BN_value_one()))
413 goto err;
414
415 /*
416 * Factorize n - 1 = k * 2^s.
417 */
418
419 s = 0;
420 while (!BN_is_bit_set(n_minus_one, s))
421 s++;
422 if (!BN_rshift(k, n_minus_one, s))
423 goto err;
424
425 /*
426 * Montgomery setup for n.
427 */
428
429 if ((mctx = BN_MONT_CTX_new()) == NULL)
430 goto err;
431
432 if (!BN_MONT_CTX_set(mctx, n, ctx))
433 goto err;
434
435 /*
436 * Perform a Miller-Rabin test for base 2 as required by BPSW.
437 */
438
439 if (!BN_set_word(base, 2))
440 goto err;
441
442 if (!bn_fermat(is_pseudoprime, n, n_minus_one, k, s, base, ctx, mctx))
443 goto err;
444 if (!*is_pseudoprime)
445 goto done;
446
447 /*
448 * Perform Miller-Rabin tests with random 3 <= base < n - 1 to reduce
449 * risk of false positives in BPSW.
450 */
451
452 for (i = 0; i < rounds; i++) {
453 if (!bn_rand_interval(base, 3, n_minus_one))
454 goto err;
455
456 if (!bn_fermat(is_pseudoprime, n, n_minus_one, k, s, base, ctx,
457 mctx))
458 goto err;
459 if (!*is_pseudoprime)
460 goto done;
461 }
462
463 /*
464 * If we got here, we have a Miller-Rabin pseudoprime.
465 */
466
467 *is_pseudoprime = 1;
468
469 done:
470 ret = 1;
471
472 err:
473 BN_MONT_CTX_free(mctx);
474 BN_CTX_end(ctx);
475
476 return ret;
477 }
478
479 /*
480 * The Baillie-Pomerance-Selfridge-Wagstaff algorithm combines a Miller-Rabin
481 * test for base 2 with a Strong Lucas pseudoprime test.
482 */
483
484 int
bn_is_prime_bpsw(int * is_pseudoprime,const BIGNUM * n,BN_CTX * in_ctx,size_t rounds)485 bn_is_prime_bpsw(int *is_pseudoprime, const BIGNUM *n, BN_CTX *in_ctx,
486 size_t rounds)
487 {
488 BN_CTX *ctx = NULL;
489 BN_ULONG mod;
490 int i;
491 int ret = 0;
492
493 if (BN_is_word(n, 2)) {
494 *is_pseudoprime = 1;
495 goto done;
496 }
497
498 if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
499 *is_pseudoprime = 0;
500 goto done;
501 }
502
503 /* Trial divisions with the first 2048 primes. */
504 for (i = 0; i < NUMPRIMES; i++) {
505 if ((mod = BN_mod_word(n, primes[i])) == (BN_ULONG)-1)
506 goto err;
507 if (mod == 0) {
508 *is_pseudoprime = BN_is_word(n, primes[i]);
509 goto done;
510 }
511 }
512
513 if ((ctx = in_ctx) == NULL)
514 ctx = BN_CTX_new();
515 if (ctx == NULL)
516 goto err;
517
518 if (!bn_miller_rabin(is_pseudoprime, n, ctx, rounds))
519 goto err;
520 if (!*is_pseudoprime)
521 goto done;
522
523 if (!bn_strong_lucas_selfridge(is_pseudoprime, n, ctx))
524 goto err;
525
526 done:
527 ret = 1;
528
529 err:
530 if (ctx != in_ctx)
531 BN_CTX_free(ctx);
532
533 return ret;
534 }
535