Add Miller-Rabin test for random bases to BPSW

The behavior of the BPSW primality test 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 infinitely 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 adds a new implementation of the old MR test that runs before running
the strong Lucas test. I like to imagine that it's slightly cleaner code.
We're effectively at about twice the cost of what we had a year ago. In
addition, it adds some non-determinism in case there actually are false
positives for the BPSW test.

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, that's not worth
much.

This is a first step towards incorporating some of the considerations in
"A performant misuse-resistant API for Primality Testing" by Massimo and
Paterson. Further work will happen in tree. In particular, there are plans
to crank the number of Miller-Rabin tests considerably so as to have a
guaranteed baseline. The manual will be updated shortly.

positive feedback beck
ok jsing

3 files changed, 130 insertions, 33 deletions

diff --git a/src/lib/libcrypto/bn/bn_bpsw.c b/src/lib/libcrypto/bn/bn_bpsw.c
index 9220339f19..a1acbbf1e9 100644
--- a/src/lib/libcrypto/bn/bn_bpsw.c
+++ b/src/lib/libcrypto/bn/bn_bpsw.c
@@ -1,4 +1,4 @@
-/*      $OpenBSD: bn_bpsw.c,v 1.8 2022/11/26 16:08:51 tb Exp $ */
+/*      $OpenBSD: bn_bpsw.c,v 1.9 2023/05/10 12:21:55 tb Exp $ */
 /*
  * Copyright (c) 2022 Martin Grenouilloux <martin.grenouilloux@lse.epita.fr>
  * Copyright (c) 2022 Theo Buehler <tb@openbsd.org>
@@ -301,23 +301,103 @@ bn_strong_lucas_selfridge(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
 }
 
 /*
- * 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 power = (base^k)^(2^(s-1)) by successive squaring of base^k.
+ *
+ * If power ever reaches -1, base^(n-1) is equal to 1 and n is a pseudoprime
+ * for base. If power reaches 1 before -1 during successive squaring, we have
+ * an unexpected square root of 1 and n is composite. Otherwise base^(n-1) != 1,
+ * and 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, const BIGNUM *base, BN_CTX *ctx, BN_MONT_CTX *mctx)
 {
-        BIGNUM *n_minus_one, *k, *x;
-        int i, s;
+        BIGNUM *power;
         int ret = 0;
+        int i;
 
         BN_CTX_start(ctx);
 
-        if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
+        if ((power = BN_CTX_get(ctx)) == NULL)
+                goto err;
+
+        /* 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(power, base, k, n, ctx, mctx))
+                goto err;
+
+        if (BN_is_one(power) || BN_cmp(power, n_minus_one) == 0) {
+                *is_prime = 1;
+                goto done;
+        }
+
+        /* Loop invariant: power is neither 1 nor -1 (mod n). */
+        for (i = 1; i < s; i++) {
+                if (!BN_mod_sqr(power, power, n, ctx))
+                        goto err;
+
+                /* n is a pseudoprime for base. */
+                if (BN_cmp(power, 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(power)) {
+                        *is_prime = 0;
+                        goto done;
+                }
+        }
+
+        /*
+         * If we get here, n is definitely composite: base^(n-1) != 1.
+         */
+
+        *is_prime = 0;
+
+ done:
+        ret = 1;
+
+ err:
+        BN_CTX_end(ctx);
+
+        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 ((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 +424,56 @@ bn_miller_rabin_base_2(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
                 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 random 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 +485,7 @@ bn_miller_rabin_base_2(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
  */
 
 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,13 +517,11 @@ bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx)
         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;
 
diff --git a/src/lib/libcrypto/bn/bn_local.h b/src/lib/libcrypto/bn/bn_local.h
index 24d91af462..78b4157d12 100644
--- a/src/lib/libcrypto/bn/bn_local.h
+++ b/src/lib/libcrypto/bn/bn_local.h
@@ -1,4 +1,4 @@
-/* $OpenBSD: bn_local.h,v 1.21 2023/04/25 17:59:41 tb Exp $ */
+/* $OpenBSD: bn_local.h,v 1.22 2023/05/10 12:21:55 tb Exp $ */
 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
  * All rights reserved.
  *
@@ -324,7 +324,7 @@ int bn_copy(BIGNUM *dst, const BIGNUM *src);
 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 */
diff --git a/src/lib/libcrypto/bn/bn_prime.c b/src/lib/libcrypto/bn/bn_prime.c
index c2fd0fc2e9..b8f0eb69ce 100644
--- a/src/lib/libcrypto/bn/bn_prime.c
+++ b/src/lib/libcrypto/bn/bn_prime.c
@@ -1,4 +1,4 @@
-/* $OpenBSD: bn_prime.c,v 1.31 2023/04/25 19:57:59 tb Exp $ */
+/* $OpenBSD: bn_prime.c,v 1.32 2023/05/10 12:21:55 tb Exp $ */
 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
  * All rights reserved.
  *
@@ -195,12 +195,12 @@ BN_generate_prime_ex(BIGNUM *ret, int bits, int safe, const BIGNUM *add,
                 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 bits, int safe, const BIGNUM *add,
                 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 checks, BN_CTX *ctx_passed,
 {
         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;