Hello,
I'm suggesting you a diff against master of the implementation of the
Baillie-PSW algorithm for primality testing.
The code in itself is commented and has information about what is being done
and why.
The reason for this change is that in this paper from 2018
(https://eprint.iacr.org/2018/749.pdf),
researchers pointed out the weaknesses of main algorithms for primality testing
where one could generate
pseudoprimes to those tests resulting in fooling them and having composite
numbers being taken as primes.
There is one that has been studied and that has no known pseudoprime and it is
called the Baillie-PSW algorithm.
It is this one I and Theo Buehler have implemented as it is safer.
Me and Theo Buehler have been actively working on it for the past month to
finally present it to you.
In hope this finds you well,
Regards,
Martin Grenouilloux.
EPITA - LSE
diff --git a/lib/libcrypto/Makefile b/lib/libcrypto/Makefile
index d6432cdc518..7e3a07c5fba 100644
--- a/lib/libcrypto/Makefile
+++ b/lib/libcrypto/Makefile
@@ -89,6 +89,7 @@ SRCS+= bn_print.c bn_rand.c bn_shift.c bn_word.c bn_blind.c
SRCS+= bn_kron.c bn_sqrt.c bn_gcd.c bn_prime.c bn_err.c bn_sqr.c
SRCS+= bn_recp.c bn_mont.c bn_mpi.c bn_exp2.c bn_gf2m.c bn_nist.c
SRCS+= bn_depr.c bn_const.c bn_x931p.c
+SRCS+= bn_bpsw.c bn_isqrt.c
# buffer/
SRCS+= buffer.c buf_err.c buf_str.c
diff --git a/lib/libcrypto/bn/bn_bpsw.c b/lib/libcrypto/bn/bn_bpsw.c
new file mode 100644
index 00000000000..d198899b2a4
--- /dev/null
+++ b/lib/libcrypto/bn/bn_bpsw.c
@@ -0,0 +1,401 @@
+/* $OpenBSD$ */
+/*
+ * Copyright (c) 2022 Martin Grenouilloux <martin.grenouill...@lse.epita.fr>
+ * Copyright (c) 2022 Theo Buehler <t...@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <openssl/bn.h>
+
+#include "bn_lcl.h"
+#include "bn_prime.h"
+
+/*
+ * For an odd n compute a / 2 (mod n). If a is even, we can do a plain
+ * division, otherwise calculate (a + n) / 2. Then reduce (mod n).
+ */
+static int
+bn_division_by_two_mod_n(BIGNUM *r, BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
+{
+ if (!BN_is_odd(n))
+ return 0;
+
+ if (!BN_mod_ct(r, a, n, ctx))
+ return 0;
+
+ if (BN_is_odd(r)) {
+ if (!BN_add(r, r, n))
+ return 0;
+ }
+
+ if (!BN_rshift1(r, r))
+ return 0;
+
+ return 1;
+}
+
+/*
+ * Given the next binary digit of k and the current Lucas terms U and V, this
+ * helper computes the next terms in the Lucas sequence defined as follows:
+ *
+ * U' = U * V (mod n)
+ * V' = (V^2 + D * U^2) / 2 (mod n)
+ *
+ * If digit == 0, bn_lucas_step() returns U' and V'. If digit == 1, it returns
+ *
+ * U'' = (U' + V') / 2 (mod n)
+ * V'' = (V' + D * U') / 2 (mod n)
+ *
+ * Compare with FIPS 186-4, Appendix C.3.3, step 6.
+ */
+static int
+bn_lucas_step(BIGNUM *U, BIGNUM *V, int digit, const BIGNUM *D,
+ const BIGNUM *n, BN_CTX *ctx)
+{
+ BIGNUM *tmp;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((tmp = BN_CTX_get(ctx)) == NULL)
+ goto done;
+
+ /* Store D * U^2 before computing U'. */
+ if (!BN_sqr(tmp, U, ctx))
+ goto done;
+ if (!BN_mul(tmp, D, tmp, ctx))
+ goto done;
+
+ /* U' = U * V (mod n). */
+ if (!BN_mod_mul(U, U, V, n, ctx))
+ goto done;
+
+ /* V' = (V^2 + D * U^2) / 2 (mod n). */
+ if (!BN_sqr(V, V, ctx))
+ goto done;
+ if (!BN_add(V, V, tmp))
+ goto done;
+ if (!bn_division_by_two_mod_n(V, V, n, ctx))
+ goto done;
+
+ if (digit == 1) {
+ /* Store D * U' before computing U''. */
+ if (!BN_mul(tmp, D, U, ctx))
+ goto done;
+
+ /* U'' = (U' + V') / 2 (mod n). */
+ if (!BN_add(U, U, V))
+ goto done;
+ if (!bn_division_by_two_mod_n(U, U, n, ctx))
+ goto done;
+
+ /* V'' = (V' + D * U') / 2 (mod n). */
+ if (!BN_add(V, V, tmp))
+ goto done;
+ if (!bn_division_by_two_mod_n(V, V, n, ctx))
+ goto done;
+ }
+
+ ret = 1;
+
+ done:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Compute the Lucas terms U_k, V_k, see FIPS 186-4, Appendix C.3.3, steps 4-6.
+ */
+static int
+bn_lucas(BIGNUM *U, BIGNUM *V, const BIGNUM *k, const BIGNUM *D,
+ const BIGNUM *n, BN_CTX *ctx)
+{
+ int digit, i;
+ int ret = 0;
+
+ if (!BN_one(U))
+ goto done;
+ if (!BN_one(V))
+ goto done;
+
+ /*
+ * Iterate over the digits of k from MSB to LSB. Start at digit 2
+ * since the first digit is dealt with by setting U = 1 and V = 1.
+ */
+ for (i = BN_num_bits(k) - 2; i >= 0; i--) {
+ digit = BN_is_bit_set(k, i);
+
+ if (!bn_lucas_step(U, V, digit, D, n, ctx))
+ goto done;
+ }
+
+ ret = 1;
+
+ done:
+ return ret;
+}
+
+/*
+ * This is a stronger variant of the Lucas test in FIPS 186-4, Appendix C.3.3.
+ * Every strong Lucas pseudoprime n is also a Lucas pseudoprime since
+ * U_{n+1} == 0 follows from U_k == 0 or V_{k * 2^r} == 0 for 0 <= r < s.
+ */
+static int
+bn_strong_lucas_test(const BIGNUM *n, const BIGNUM *D, BN_CTX *ctx)
+{
+ BIGNUM *k, *U, *V;
+ int r, s;
+ int ret = -1;
+
+ BN_CTX_start(ctx);
+
+ if ((k = BN_CTX_get(ctx)) == NULL)
+ goto done;
+ if ((U = BN_CTX_get(ctx)) == NULL)
+ goto done;
+ if ((V = BN_CTX_get(ctx)) == NULL)
+ goto done;
+
+ /*
+ * Factorize n + 1 = k * 2^s with odd k: shift away the s trailing ones
+ * of n and set the lowest bit of the resulting number k.
+ */
+ s = 0;
+ while (BN_is_bit_set(n, s))
+ s++;
+ if (!BN_rshift(k, n, s))
+ goto done;
+ if (!BN_set_bit(k, 0))
+ goto done;
+
+ /*
+ * Calculate the Lucas terms U_k and V_k. If either of them is zero,
+ * then n is a strong Lucas pseudoprime.
+ */
+ if (!bn_lucas(U, V, k, D, n, ctx))
+ goto done;
+
+ if (BN_is_zero(U) || BN_is_zero(V)) {
+ ret = 1;
+ goto done;
+ }
+
+ /*
+ * Check if any V_{k * d^r} is zero for 1 <= r < s.
+ */
+ for (r = 1; r < s; r++) {
+ if (!bn_lucas_step(U, V, 0, D, n, ctx))
+ goto done;
+
+ if (BN_is_zero(V)) {
+ ret = 1;
+ goto done;
+ }
+ }
+
+ /* If we got here, n is definitely composite. */
+ ret = 0;
+
+ done:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Test n for primality using the strong Lucas test with Selfridge's
+ * parameters. Returns 1 if n is prime or a strong Lucas-Selfridge
+ * pseudoprime. Returns 0 if n is definitely composite.
+ */
+static int
+bn_strong_lucas_selfridge(const BIGNUM *n, BN_CTX *ctx)
+{
+ BIGNUM *D, *two;
+ int jacobi_symbol, perfect_square, sign;
+ int ret = -1;
+
+ BN_CTX_start(ctx);
+
+ /* If n is a perfect square, it is composite. */
+ if (!bn_is_square(&perfect_square, n, ctx))
+ goto done;
+ if (perfect_square) {
+ ret = 0;
+ goto done;
+ }
+
+ /*
+ * Find the first element D in the sequence 5, -7, 9, -11, 13, ...
+ * such that Jacobi(D, n) = -1 (Selfridge's algorithm).
+ */
+ if ((D = BN_CTX_get(ctx)) == NULL)
+ goto done;
+ if ((two = BN_CTX_get(ctx)) == NULL)
+ goto done;
+
+ sign = 1;
+ if (!BN_set_word(D, 5))
+ goto done;
+ if (!BN_set_word(two, 2))
+ goto done;
+
+ while (1) {
+ /* For odd n the Kronecker symbol computes the Jacobi symbol. */
+ if ((jacobi_symbol = BN_kronecker(D, n, ctx)) == -2)
+ goto done;
+
+ /* We found the value for D. */
+ if (jacobi_symbol == -1)
+ break;
+
+ /* n and D have prime factors in common. */
+ if (jacobi_symbol == 0) {
+ ret = 0;
+ goto done;
+ }
+
+ /* Subtract or add 2 to follow the sequence described above. */
+ sign = -sign;
+ if (!BN_uadd(D, D, two))
+ goto done;
+ BN_set_negative(D, sign == -1);
+ }
+
+ ret = bn_strong_lucas_test(n, D, ctx);
+
+ done:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Miller-Rabin primality test for base 2.
+ */
+static int
+bn_miller_rabin_base_2(const BIGNUM *n, BN_CTX *ctx)
+{
+ BIGNUM *n_minus_one, *k, *x;
+ int i, s;
+ int ret = -1;
+
+ BN_CTX_start(ctx);
+
+ if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
+ goto done;
+ if ((k = BN_CTX_get(ctx)) == NULL)
+ goto done;
+ if ((x = BN_CTX_get(ctx)) == NULL)
+ goto done;
+
+ if (BN_is_word(n, 2) || BN_is_word(n, 3)) {
+ ret = 1;
+ goto done;
+ }
+
+ if (BN_cmp(n, BN_value_one()) == 0 || !BN_is_odd(n)) {
+ ret = 0;
+ goto done;
+ }
+
+ if (!BN_sub(n_minus_one, n, BN_value_one()))
+ goto done;
+
+ s = 0;
+ while (!BN_is_bit_set(n_minus_one, s))
+ s++;
+ if (!BN_rshift(k, n_minus_one, s))
+ goto done;
+
+ /* If 2^k is 1 or -1 (mod n) then n is a 2-pseudoprime. */
+ if (!BN_set_word(x, 2))
+ goto done;
+ if (!BN_mod_exp_ct(x, x, k, n, ctx))
+ goto done;
+
+ if (BN_is_one(x) || BN_cmp(x, n_minus_one) == 0) {
+ ret = 1;
+ goto done;
+ }
+
+ for (i = 1; i < s; i++) {
+ if (!BN_mod_sqr(x, x, n, ctx))
+ goto done;
+ if (BN_cmp(x, n_minus_one) == 0) {
+ ret = 1;
+ goto done;
+ }
+ }
+
+ /* If we got here, n is definitely composite. */
+ ret = 0;
+
+ done:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * The Baillie-Pomerance-Selfridge-Wagstaff algorithm combines a Miller-Rabin
+ * test for base 2 with a Strong Lucas pseudoprime test.
+ */
+
+int
+bn_is_prime_bpsw(const BIGNUM *n, BN_CTX *in_ctx)
+{
+ BN_CTX *ctx = in_ctx;
+ BN_ULONG mod;
+ int i;
+ int ret = -1;
+
+ if (BN_is_word(n, 2)) {
+ ret = 1;
+ goto done;
+ }
+
+ if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
+ ret = 0;
+ goto done;
+ }
+
+ /* Trial divisions with the first 2048 primes. */
+ for (i = 0; i < NUMPRIMES; i++) {
+ if ((mod = BN_mod_word(n, primes[i])) == (BN_ULONG)-1)
+ goto done;
+ if (mod == 0) {
+ ret = BN_is_word(n, primes[i]);
+ goto done;
+ }
+ }
+
+ if (ctx == NULL)
+ ctx = BN_CTX_new();
+ if (ctx == NULL)
+ goto done;
+
+ if ((ret = bn_miller_rabin_base_2(n, ctx)) <= 0)
+ goto done;
+
+ /* XXX - Miller-Rabin for random bases? - see FIPS 186-4, Table C.1. */
+
+ ret = bn_strong_lucas_selfridge(n, ctx);
+
+ done:
+ if (ctx != in_ctx)
+ BN_CTX_free(ctx);
+
+ return ret;
+}
diff --git a/lib/libcrypto/bn/bn_isqrt.c b/lib/libcrypto/bn/bn_isqrt.c
new file mode 100644
index 00000000000..ea776636642
--- /dev/null
+++ b/lib/libcrypto/bn/bn_isqrt.c
@@ -0,0 +1,239 @@
+/* $OpenBSD$ */
+/*
+ * Copyright (c) 2022 Theo Buehler <t...@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <stddef.h>
+#include <stdint.h>
+
+#include <openssl/bn.h>
+#include <openssl/err.h>
+#include <openssl/ossl_typ.h>
+
+#include "bn_lcl.h"
+
+#define CTASSERT(x) extern char _ctassert[(x) ? 1 : -1 ] \
+ __attribute__((__unused__))
+
+/*
+ * Calculate integer square root of |n| using a variant of Newton's method.
+ *
+ * Returns the integer square root of |n| in the caller-provided |out_sqrt|;
+ * |*out_perfect| is set to 1 if and only if |n| is a perfect square.
+ * One of |out_sqrt| and |out_perfect| can be NULL; |in_ctx| can be NULL.
+ *
+ * Returns 0 on error, 1 on success.
+ *
+ * Adapted from pure Python describing cpython's math.isqrt(), without bothering
+ * with any of the optimizations in the C code. A correctness proof is here:
+ * https://github.com/mdickinson/snippets/blob/master/proofs/isqrt/src/isqrt.lean
+ * The comments in the Python code also give a rather detailed proof.
+ */
+
+int
+bn_isqrt(BIGNUM *out_sqrt, int *out_perfect, const BIGNUM *n, BN_CTX *in_ctx)
+{
+ BN_CTX *ctx = in_ctx;
+ BIGNUM *a, *b;
+ int c, d, e, s;
+ int cmp;
+ int ret = 0;
+
+ if (out_perfect == NULL && out_sqrt == NULL) {
+ BNerror(ERR_R_PASSED_NULL_PARAMETER);
+ goto err;
+ }
+
+ if (BN_is_negative(n)) {
+ BNerror(BN_R_INVALID_RANGE);
+ goto err;
+ }
+
+ if (BN_is_zero(n)) {
+ if (out_perfect != NULL)
+ *out_perfect = 1;
+ if (out_sqrt != NULL) {
+ if (!BN_zero(out_sqrt))
+ goto err;
+ }
+ return 1;
+ }
+
+ if (ctx == NULL)
+ ctx = BN_CTX_new();
+ if (ctx == NULL)
+ goto err;
+
+ BN_CTX_start(ctx);
+
+ if ((a = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((b = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ if (!BN_one(a))
+ goto err;
+
+ c = (BN_num_bits(n) - 1) / 2;
+ d = 0;
+
+ /* Calculate s = floor(log(c)). */
+ if (!BN_set_word(b, c))
+ goto err;
+ s = BN_num_bits(b) - 1;
+
+ /*
+ * By definition, the loop below is run <= floor(log(log(n))) times.
+ * Comments in the cpython code establish the loop invariant that
+ *
+ * (a - 1)^2 < n / 4^(c - d) < (a + 1)^2
+ *
+ * holds true in every iteration. Once this is proved via induction,
+ * correctness of the algorithm is easy.
+ *
+ * Roughly speaking, A = (a << (d - e)) is used for one Newton step
+ * "a = (A >> 1) + (m >> 1) / A" approximating m = (n >> 2 * (c - d)).
+ */
+
+ for (; s >= 0; s--) {
+ e = d;
+ d = c >> s;
+
+ if (!BN_rshift(b, n, 2 * c - d - e + 1))
+ goto err;
+
+ if (!BN_div_ct(b, NULL, b, a, ctx))
+ goto err;
+
+ if (!BN_lshift(a, a, d - e - 1))
+ goto err;
+
+ if (!BN_add(a, a, b))
+ goto err;
+ }
+
+ /*
+ * The loop invariant implies that either a or a - 1 is isqrt(n).
+ * Figure out which one it is. The invariant also implies that for
+ * a perfect square n, a must be the square root.
+ */
+
+ if (!BN_sqr(b, a, ctx))
+ goto err;
+
+ /* If a^2 > n, we must have isqrt(n) == a - 1. */
+ if ((cmp = BN_cmp(b, n)) > 0) {
+ if (!BN_sub_word(a, 1))
+ goto err;
+ }
+
+ if (out_perfect != NULL)
+ *out_perfect = cmp == 0;
+
+ if (out_sqrt != NULL) {
+ if (!BN_copy(out_sqrt, a))
+ goto err;
+ }
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+ if (ctx != in_ctx)
+ BN_CTX_free(ctx);
+
+ return ret;
+}
+
+/*
+ * is_square_mod_N[r % N] indicates whether r % N has a square root modulo N.
+ * The tables are generated in regress/lib/libcrypto/bn/bn_isqrt.c.
+ */
+
+static const uint8_t is_square_mod_11[] = {
+ 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0,
+};
+CTASSERT(sizeof(is_square_mod_11) == 11);
+
+static const uint8_t is_square_mod_63[] = {
+ 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0,
+ 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0,
+ 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
+ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
+};
+CTASSERT(sizeof(is_square_mod_63) == 63);
+
+static const uint8_t is_square_mod_64[] = {
+ 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
+ 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
+ 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
+ 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
+};
+CTASSERT(sizeof(is_square_mod_64) == 64);
+
+static const uint8_t is_square_mod_65[] = {
+ 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0,
+ 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0,
+ 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
+ 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0,
+ 1,
+};
+CTASSERT(sizeof(is_square_mod_65) == 65);
+
+#define BN_lsw(n) (((n)->top == 0) ? (BN_ULONG) 0 : (n)->d[0])
+
+/*
+ * Determine whether n is a perfect square or not.
+ *
+ * Returns 1 on success and 0 on error. In case of success, |*is_square| is
+ * set to 1 if and only if |n| is a perfect square.
+ */
+
+int
+bn_is_square(int *is_square, const BIGNUM *n, BN_CTX *ctx)
+{
+ BN_ULONG r;
+
+ *is_square = 0;
+
+ if (BN_is_negative(n))
+ return 1;
+
+ /*
+ * Before performing an expensive bn_isqrt() operation, weed out many
+ * obvious non-squares. See H. Cohen, "A course in computational
+ * algebraic number theory", Algorithm 1.7.3.
+ *
+ * The idea is that a square remains a square when reduced modulo any
+ * number. The moduli are chosen in such a way that a non-square has
+ * probability < 1% of passing the four table lookups.
+ */
+
+ /* n % 64 */
+ r = BN_lsw(n) & 0x3f;
+
+ if (!is_square_mod_64[r % 64])
+ return 1;
+
+ if ((r = BN_mod_word(n, 11 * 63 * 65)) == (BN_ULONG)-1)
+ return 0;
+
+ if (!is_square_mod_63[r % 63] ||
+ !is_square_mod_65[r % 65] ||
+ !is_square_mod_11[r % 11])
+ return 1;
+
+ return bn_isqrt(NULL, is_square, n, ctx);
+}
diff --git a/lib/libcrypto/bn/bn_lcl.h b/lib/libcrypto/bn/bn_lcl.h
index d8c9e20f402..872e15ea0ea 100644
--- a/lib/libcrypto/bn/bn_lcl.h
+++ b/lib/libcrypto/bn/bn_lcl.h
@@ -653,5 +653,10 @@ int BN_gcd_nonct(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx);
int BN_swap_ct(BN_ULONG swap, BIGNUM *a, BIGNUM *b, size_t nwords);
+int bn_isqrt(BIGNUM *out_sqrt, int *out_perfect, const BIGNUM *n, BN_CTX *ctx);
+int bn_is_square(int *is_square, const BIGNUM *n, BN_CTX *ctx);
+
+int bn_is_prime_bpsw(const BIGNUM *n, BN_CTX *in_ctx);
+
__END_HIDDEN_DECLS
#endif
diff --git a/lib/libcrypto/bn/bn_prime.c b/lib/libcrypto/bn/bn_prime.c
index e78c5686ab5..3eee43bf1bd 100644
--- a/lib/libcrypto/bn/bn_prime.c
+++ b/lib/libcrypto/bn/bn_prime.c
@@ -259,12 +259,16 @@ int
BN_is_prime_fasttest_ex(const BIGNUM *a, int checks, BN_CTX *ctx_passed,
int do_trial_division, BN_GENCB *cb)
{
- int i, j, ret = -1;
- int k;
BN_CTX *ctx = NULL;
BIGNUM *A1, *A1_odd, *check; /* taken from ctx */
BN_MONT_CTX *mont = NULL;
const BIGNUM *A = NULL;
+ int i, j, k;
+ int ret = -1;
+
+#ifdef LIBRESSL_HAS_BPSW
+ return bn_is_prime_bpsw(a, ctx_passed);
+#endif
if (BN_cmp(a, BN_value_one()) <= 0)
return 0;
diff --git a/lib/libcrypto/opensslfeatures.h b/lib/libcrypto/opensslfeatures.h
index 49a5f15b597..56af8d0713e 100644
--- a/lib/libcrypto/opensslfeatures.h
+++ b/lib/libcrypto/opensslfeatures.h
@@ -5,6 +5,7 @@
*/
#define LIBRESSL_HAS_TLS1_3
#define LIBRESSL_HAS_DTLS1_2
+#define LIBRESSL_HAS_BPSW
#define OPENSSL_THREADS
