On Fri, Apr 28, 2023 at 10:23:15AM +0200, Theo Buehler wrote:
> 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;
Nit: you call this a pseudo-prime - which I agree with, but the variable you
have added
is "is_prime" - maybe call it "maybe-prime" or "pseudo-prime"
>
> 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.
You are now pre-checking with miller rabin to reduce the likelyhood of
a composite number of any size, but you are not explaining this here.
The point of the miller-rabin pre-test is to cheaply reduce the likelyhood
of a composite number being fed to Ballie-PSW to "very small" so that even
if the future hypothesis finds such numbers we are still protected,
and this should be explained here.
> +.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
