The behavior of BPSW for numbers > 2^64 is not very well understood.
While there is no known composite that passes the test, there are
heuristics that indicate that there are likely many. Therefore it seems
appropriate to harden the test. Having a settable number of MR rounds
before doing a version of BPSW is also the approach taken by Go's
primality check in math/big.
This is a first step towards incorporating some of the considerations in
"A performant misuse-resistant API for Primality Testing" by Massimo and
Paterson. It should be noted that it is based on this paper (and due to
FIPS, I suppose) that OpenSSL 3 made their primality check so slow that
the bn_prime.c regress test took over 10 minutes on an M3000 before the
libcrypto bump (it took less than a minute with ours).
This basically adds back an equivalent of the old inadequate prime number
test before running the strong Lucas test, with slightly cleaner code.
While I don't think performance matters a lot here, we're then effectively
at about twice the cost of what we had a year ago. I also think that having
some non-determinism in the algorithm is a plus.
The implementation is straightforward. It could easily be tweaked to use
the additional gcds in the "enhanced" MR test of FIPS 186-5, but as long
as we are only going to throw away the additional info, this doesn't add
all that much.
I'd like to land this and then do further work in tree. We probably want
to crank the number of MR rounds done by default considerably. I will
also need to undo some of the clean-up done in BN_generate_prime_ex()
after BPSW landed.
Index: bn/bn_bpsw.c
===================================================================
RCS file: /cvs/src/lib/libcrypto/bn/bn_bpsw.c,v
retrieving revision 1.8
diff -u -p -r1.8 bn_bpsw.c
--- bn/bn_bpsw.c 26 Nov 2022 16:08:51 -0000 1.8
+++ bn/bn_bpsw.c 28 Apr 2023 07:56:01 -0000
@@ -301,23 +301,94 @@ bn_strong_lucas_selfridge(int *is_prime,
}
/*
- * Miller-Rabin primality test for base 2.
+ * Fermat criterion in Miller-Rabin test.
+ *
+ * Check whether 1 < base < n - 1 witnesses that n is composite. For prime n:
+ *
+ * * Fermat's little theorem: base^(n-1) = 1 (mod n).
+ * * The only square roots of 1 (mod n) are 1 and -1.
+ *
+ * Calculate base^((n-1)/2) by writing n - 1 = k * 2^s with odd k: iteratively
+ * compute (base^k)^(2^(s-1)) by successive squaring of base^k.
+ *
+ * If -1 is ever reached, base^(n-1) works out to 1 and n is a pseudoprime
+ * for base. If 1 is reached before -1, we have an unexpected square root of
+ * unity and n is composite. Otherwise base^(n-1) != 1, so n is composite.
*/
static int
-bn_miller_rabin_base_2(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
+bn_fermat(int *is_prime, const BIGNUM *n, const BIGNUM *n_minus_one,
+ const BIGNUM *k, int s, BIGNUM *base, BN_CTX *ctx, BN_MONT_CTX *mctx)
{
- BIGNUM *n_minus_one, *k, *x;
- int i, s;
+ int ret = 0;
+ int i;
+
+ /* Sanity check: ensure that 1 < base < n - 1. */
+ if (BN_cmp(base, BN_value_one()) <= 0 || BN_cmp(base, n_minus_one) >= 0)
+ goto err;
+
+ if (!BN_mod_exp_mont_ct(base, base, k, n, ctx, mctx))
+ goto err;
+
+ if (BN_is_one(base) || BN_cmp(base, n_minus_one) == 0) {
+ *is_prime = 1;
+ goto done;
+ }
+
+ /* Loop invariant: base is neither 1 nor -1 (mod n). */
+ for (i = 1; i < s; i++) {
+ if (!BN_mod_sqr(base, base, n, ctx))
+ goto err;
+
+ /* n is a pseudoprime for base. */
+ if (BN_cmp(base, n_minus_one) == 0) {
+ *is_prime = 1;
+ goto done;
+ }
+
+ /* n is composite: there's a square root of unity != 1 or -1. */
+ if (BN_is_one(base)) {
+ *is_prime = 0;
+ goto done;
+ }
+ }
+
+ /*
+ * If we get here, n is definitely composite: base^(n-1) != 1.
+ */
+
+ *is_prime = 0;
+
+ done:
+ ret = 1;
+
+ err:
+ return ret;
+}
+
+/*
+ * Miller-Rabin primality test for base 2 and for |rounds| of random bases.
+ * On success: is_prime == 0 implies n is composite - the converse is false.
+ */
+
+static int
+bn_miller_rabin(int *is_prime, const BIGNUM *n, BN_CTX *ctx, size_t rounds)
+{
+ BN_MONT_CTX *mctx = NULL;
+ BIGNUM *base, *k, *n_minus_one, *three;
+ size_t i;
+ int s;
int ret = 0;
BN_CTX_start(ctx);
- if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
+ if ((base = BN_CTX_get(ctx)) == NULL)
goto err;
if ((k = BN_CTX_get(ctx)) == NULL)
goto err;
- if ((x = BN_CTX_get(ctx)) == NULL)
+ if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((three = BN_CTX_get(ctx)) == NULL)
goto err;
if (BN_is_word(n, 2) || BN_is_word(n, 3)) {
@@ -344,43 +415,56 @@ bn_miller_rabin_base_2(int *is_prime, co
goto err;
/*
- * If 2^k is 1 or -1 (mod n) then n is a 2-pseudoprime.
+ * Montgomery setup for n.
*/
- if (!BN_set_word(x, 2))
+ if ((mctx = BN_MONT_CTX_new()) == NULL)
goto err;
- if (!BN_mod_exp_ct(x, x, k, n, ctx))
+
+ if (!BN_MONT_CTX_set(mctx, n, ctx))
goto err;
- if (BN_is_one(x) || BN_cmp(x, n_minus_one) == 0) {
- *is_prime = 1;
+ /*
+ * Perform a Miller-Rabin test for base 2 as required by BPSW.
+ */
+
+ if (!BN_set_word(base, 2))
+ goto err;
+
+ if (!bn_fermat(is_prime, n, n_minus_one, k, s, base, ctx, mctx))
+ goto err;
+ if (!*is_prime)
goto done;
- }
/*
- * If 2^{2^i k} == -1 (mod n) for some 1 <= i < s, then n is a
- * 2-pseudoprime.
+ * Perform Miller-Rabin tests with 3 <= base < n - 1 to reduce risk of
+ * false positives in BPSW.
*/
- for (i = 1; i < s; i++) {
- if (!BN_mod_sqr(x, x, n, ctx))
+ if (!BN_set_word(three, 3))
+ goto err;
+
+ for (i = 0; i < rounds; i++) {
+ if (!bn_rand_interval(base, three, n_minus_one))
goto err;
- if (BN_cmp(x, n_minus_one) == 0) {
- *is_prime = 1;
+
+ if (!bn_fermat(is_prime, n, n_minus_one, k, s, base, ctx, mctx))
+ goto err;
+ if (!*is_prime)
goto done;
- }
}
/*
- * If we got here, n is definitely composite.
+ * If we got here, we have a Miller-Rabin pseudoprime.
*/
- *is_prime = 0;
+ *is_prime = 1;
done:
ret = 1;
err:
+ BN_MONT_CTX_free(mctx);
BN_CTX_end(ctx);
return ret;
@@ -392,7 +476,7 @@ bn_miller_rabin_base_2(int *is_prime, co
*/
int
-bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx)
+bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx, size_t rounds)
{
BN_CTX *ctx = NULL;
BN_ULONG mod;
@@ -424,12 +508,10 @@ bn_is_prime_bpsw(int *is_prime, const BI
if (ctx == NULL)
goto err;
- if (!bn_miller_rabin_base_2(is_prime, n, ctx))
+ if (!bn_miller_rabin(is_prime, n, ctx, rounds))
goto err;
if (!*is_prime)
goto done;
-
- /* XXX - Miller-Rabin for random bases? See FIPS 186-4, Table C.1. */
if (!bn_strong_lucas_selfridge(is_prime, n, ctx))
goto err;
Index: bn/bn_local.h
===================================================================
RCS file: /cvs/src/lib/libcrypto/bn/bn_local.h,v
retrieving revision 1.21
diff -u -p -r1.21 bn_local.h
--- bn/bn_local.h 25 Apr 2023 17:59:41 -0000 1.21
+++ bn/bn_local.h 28 Apr 2023 07:56:01 -0000
@@ -324,7 +324,7 @@ int bn_copy(BIGNUM *dst, const BIGNUM *s
int bn_isqrt(BIGNUM *out_sqrt, int *out_perfect, const BIGNUM *n, BN_CTX *ctx);
int bn_is_perfect_square(int *out_perfect, const BIGNUM *n, BN_CTX *ctx);
-int bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx);
+int bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *ctx, size_t
rounds);
__END_HIDDEN_DECLS
#endif /* !HEADER_BN_LOCAL_H */
Index: bn/bn_prime.c
===================================================================
RCS file: /cvs/src/lib/libcrypto/bn/bn_prime.c,v
retrieving revision 1.31
diff -u -p -r1.31 bn_prime.c
--- bn/bn_prime.c 25 Apr 2023 19:57:59 -0000 1.31
+++ bn/bn_prime.c 28 Apr 2023 07:56:01 -0000
@@ -195,12 +195,12 @@ BN_generate_prime_ex(BIGNUM *ret, int bi
goto err;
if (!safe) {
- if (!bn_is_prime_bpsw(&is_prime, ret, ctx))
+ if (!bn_is_prime_bpsw(&is_prime, ret, ctx, 1))
goto err;
if (!is_prime)
goto loop;
} else {
- if (!bn_is_prime_bpsw(&is_prime, ret, ctx))
+ if (!bn_is_prime_bpsw(&is_prime, ret, ctx, 1))
goto err;
if (!is_prime)
goto loop;
@@ -213,7 +213,7 @@ BN_generate_prime_ex(BIGNUM *ret, int bi
if (!BN_rshift1(p, ret))
goto err;
- if (!bn_is_prime_bpsw(&is_prime, p, ctx))
+ if (!bn_is_prime_bpsw(&is_prime, p, ctx, 1))
goto err;
if (!is_prime)
goto loop;
@@ -243,8 +243,14 @@ BN_is_prime_fasttest_ex(const BIGNUM *a,
{
int is_prime;
+ if (checks < 0)
+ return -1;
+
+ if (checks == BN_prime_checks)
+ checks = BN_prime_checks_for_size(BN_num_bits(a));
+
/* XXX - tickle BN_GENCB in bn_is_prime_bpsw(). */
- if (!bn_is_prime_bpsw(&is_prime, a, ctx_passed))
+ if (!bn_is_prime_bpsw(&is_prime, a, ctx_passed, checks))
return -1;
return is_prime;
Index: man/BN_generate_prime.3
===================================================================
RCS file: /cvs/src/lib/libcrypto/man/BN_generate_prime.3,v
retrieving revision 1.20
diff -u -p -r1.20 BN_generate_prime.3
--- man/BN_generate_prime.3 24 Nov 2022 19:06:38 -0000 1.20
+++ man/BN_generate_prime.3 28 Apr 2023 07:56:01 -0000
@@ -84,7 +84,7 @@
.Nm BN_is_prime ,
.Nm BN_is_prime_fasttest
.\" Nm BN_prime_checks_for_size is intentionally undocumented
-.\" because it is no longer used by LibreSSL.
+.\" because it should not be used outside of libcrypto.
.Nd generate primes and test for primality
.Sh SYNOPSIS
.In openssl/bn.h
@@ -178,18 +178,33 @@ test whether the number
.Fa a
is prime.
In LibreSSL, both functions behave identically,
-use the Baillie-Pomerance-Selfridge-Wagstaff algorithm,
-and ignore the
+use the Baillie-Pomerance-Selfridge-Wagstaff algorithm
+combined with
.Fa checks
-and
+rounds of the Miller-Rabin test.
+They ignore the
.Fa do_trial_division
-arguments.
+argument.
.Pp
It is unknown whether any composite number exists that the
Baillie-PSW algorithm misclassifies as a prime.
Some suspect that there may be infinitely many such numbers,
but not a single one is currently known.
It is known that no such number exists below 2\(ha64.
+.Pp
+The number of Miller-Rabin rounds performed by
+.Fn BN_is_prime_fasttest_ex
+and
+.Fn BN_is_prime_ex
+is determined by
+.Fa checks :
+If it is positive, it is used as the number of rounds.
+If
+.Fa checks
+is zero or the special value
+.Dv BN_prime_checks ,
+a suitable number of rounds is calculated from the bit length of
+.Fa a .
.Pp
If
.Dv NULL
