Here is an updated diff that addresses some of the changes suggested and
requested by jsing.
The main change is that we now pass an is_prime out parameter through
and distinguish between goto done and goto err. All functions now have
boolean return values.
In addition, some functions have been renamed and some comments tweaked.
I also fixed a stupid bug I had in bn_isqrt() that mishandled the BN_CTX
on unlikely early error.
Index: Makefile
===================================================================
RCS file: /cvs/src/lib/libcrypto/Makefile,v
retrieving revision 1.74
diff -u -p -r1.74 Makefile
--- Makefile 8 May 2022 20:59:32 -0000 1.74
+++ Makefile 12 Jul 2022 10:41:09 -0000
@@ -89,6 +89,7 @@ SRCS+= bn_print.c bn_rand.c bn_shift.c b
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
Index: bn/bn_bpsw.c
===================================================================
RCS file: bn/bn_bpsw.c
diff -N bn/bn_bpsw.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ bn/bn_bpsw.c 12 Jul 2022 10:41:09 -0000
@@ -0,0 +1,420 @@
+/* $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_div_by_two_mod_odd_n(BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
+{
+ if (!BN_is_odd(n))
+ return 0;
+
+ if (BN_is_odd(a)) {
+ if (!BN_add(a, a, n))
+ return 0;
+ }
+ if (!BN_rshift1(a, a))
+ return 0;
+ if (!BN_mod_ct(a, a, n, ctx))
+ 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 err;
+
+ /* Calculate D * U^2 before computing U'. */
+ if (!BN_sqr(tmp, U, ctx))
+ goto err;
+ if (!BN_mul(tmp, D, tmp, ctx))
+ goto err;
+
+ /* U' = U * V (mod n). */
+ if (!BN_mod_mul(U, U, V, n, ctx))
+ goto err;
+
+ /* V' = (V^2 + D * U^2) / 2 (mod n). */
+ if (!BN_sqr(V, V, ctx))
+ goto err;
+ if (!BN_add(V, V, tmp))
+ goto err;
+ if (!bn_div_by_two_mod_odd_n(V, n, ctx))
+ goto err;
+
+ if (digit == 1) {
+ /* Calculate D * U' before computing U''. */
+ if (!BN_mul(tmp, D, U, ctx))
+ goto err;
+
+ /* U'' = (U' + V') / 2 (mod n). */
+ if (!BN_add(U, U, V))
+ goto err;
+ if (!bn_div_by_two_mod_odd_n(U, n, ctx))
+ goto err;
+
+ /* V'' = (V' + D * U') / 2 (mod n). */
+ if (!BN_add(V, V, tmp))
+ goto err;
+ if (!bn_div_by_two_mod_odd_n(V, n, ctx))
+ goto err;
+ }
+
+ ret = 1;
+
+ err:
+ 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 err;
+ if (!BN_one(V))
+ goto err;
+
+ /*
+ * 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 err;
+ }
+
+ ret = 1;
+
+ err:
+ 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(int *is_prime, const BIGNUM *n, const BIGNUM *D,
+ BN_CTX *ctx)
+{
+ BIGNUM *k, *U, *V;
+ int r, s;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((k = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((U = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((V = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /*
+ * 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 err;
+ if (!BN_set_bit(k, 0))
+ goto err;
+
+ /*
+ * 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 err;
+
+ if (BN_is_zero(U) || BN_is_zero(V)) {
+ *is_prime = 1;
+ goto done;
+ }
+
+ /*
+ * If any V_{k * 2^r} is zero for 1 <= r < s then n is a strong Lucas
+ * pseudoprime.
+ */
+ for (r = 1; r < s; r++) {
+ if (!bn_lucas_step(U, V, 0, D, n, ctx))
+ goto err;
+
+ if (BN_is_zero(V)) {
+ *is_prime = 1;
+ goto done;
+ }
+ }
+
+ /* If we got here, n is definitely composite. */
+ *is_prime = 0;
+
+ done:
+ ret = 1;
+
+ err:
+ 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(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
+{
+ BIGNUM *D, *two;
+ int is_perfect_square, jacobi_symbol, sign;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ /* If n is a perfect square, it is composite. */
+ if (!bn_is_perfect_square(&is_perfect_square, n, ctx))
+ goto err;
+ if (is_perfect_square) {
+ *is_prime = 0;
+ goto err;
+ }
+
+ /*
+ * Find the first D in the Selfridge sequence 5, -7, 9, -11, 13, ...
+ * such that the Jacobi symbol (D/n) is -1.
+ */
+ if ((D = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((two = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ sign = 1;
+ if (!BN_set_word(D, 5))
+ goto err;
+ if (!BN_set_word(two, 2))
+ goto err;
+
+ while (1) {
+ /* For odd n the Kronecker symbol computes the Jacobi symbol. */
+ if ((jacobi_symbol = BN_kronecker(D, n, ctx)) == -2)
+ goto err;
+
+ /* We found the value for D. */
+ if (jacobi_symbol == -1)
+ break;
+
+ /* n and D have prime factors in common. */
+ if (jacobi_symbol == 0) {
+ *is_prime = 0;
+ goto done;
+ }
+
+ sign = -sign;
+ if (!BN_uadd(D, D, two))
+ goto err;
+ BN_set_negative(D, sign == -1);
+ }
+
+ if (!bn_strong_lucas_test(is_prime, n, D, ctx))
+ goto err;
+
+ done:
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Miller-Rabin primality test for base 2.
+ */
+static int
+bn_miller_rabin_base_2(int *is_prime, const BIGNUM *n, BN_CTX *ctx)
+{
+ BIGNUM *n_minus_one, *k, *x;
+ int i, s;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((n_minus_one = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((k = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((x = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ if (BN_is_word(n, 2) || BN_is_word(n, 3)) {
+ *is_prime = 1;
+ goto done;
+ }
+
+ if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
+ *is_prime = 0;
+ goto done;
+ }
+
+ if (!BN_sub(n_minus_one, n, BN_value_one()))
+ goto err;
+
+ s = 0;
+ while (!BN_is_bit_set(n_minus_one, s))
+ s++;
+ if (!BN_rshift(k, n_minus_one, s))
+ goto err;
+
+ /* If 2^k is 1 or -1 (mod n) then n is a 2-pseudoprime. */
+ if (!BN_set_word(x, 2))
+ goto err;
+ if (!BN_mod_exp_ct(x, x, k, n, ctx))
+ goto err;
+
+ if (BN_is_one(x) || BN_cmp(x, n_minus_one) == 0) {
+ *is_prime = 1;
+ goto done;
+ }
+
+ /*
+ * If 2^{2^i k} == -1 (mod n) for some 1 <= i < s, then n is a
+ * 2-pseudoprime
+ */
+ for (i = 1; i < s; i++) {
+ if (!BN_mod_sqr(x, x, n, ctx))
+ goto err;
+ if (BN_cmp(x, n_minus_one) == 0) {
+ *is_prime = 1;
+ goto done;
+ }
+ }
+
+ /* If we got here, n is definitely composite. */
+ *is_prime = 0;
+
+ done:
+ ret = 1;
+
+ err:
+ 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(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx)
+{
+ BN_CTX *ctx = in_ctx;
+ BN_ULONG mod;
+ int i;
+ int ret = 0;
+
+ if (BN_is_word(n, 2)) {
+ *is_prime = 1;
+ goto done;
+ }
+
+ if (BN_cmp(n, BN_value_one()) <= 0 || !BN_is_odd(n)) {
+ *is_prime = 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 err;
+ if (mod == 0) {
+ *is_prime = BN_is_word(n, primes[i]);
+ goto done;
+ }
+ }
+
+ if (ctx == NULL)
+ ctx = BN_CTX_new();
+ if (ctx == NULL)
+ goto err;
+
+ if (!bn_miller_rabin_base_2(is_prime, n, ctx))
+ 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;
+
+ done:
+ ret = 1;
+
+ err:
+ if (ctx != in_ctx)
+ BN_CTX_free(ctx);
+
+ return ret;
+}
Index: bn/bn_isqrt.c
===================================================================
RCS file: bn/bn_isqrt.c
diff -N bn/bn_isqrt.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ bn/bn_isqrt.c 12 Jul 2022 10:41:09 -0000
@@ -0,0 +1,237 @@
+/* $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 "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 = NULL;
+ BIGNUM *a, *b;
+ int c, d, e, s;
+ int cmp, perfect;
+ 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 ((ctx = in_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_is_zero(n)) {
+ perfect = 1;
+ if (!BN_zero(a))
+ goto err;
+ goto done;
+ }
+
+ 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;
+ }
+
+ perfect = cmp == 0;
+
+ done:
+ if (out_perfect != NULL)
+ *out_perfect = perfect;
+
+ 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);
+
+/*
+ * 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_perfect_square(int *is_perfect_square, const BIGNUM *n, BN_CTX *ctx)
+{
+ BN_ULONG r;
+
+ *is_perfect_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_perfect_square, n, ctx);
+}
Index: bn/bn_kron.c
===================================================================
RCS file: /cvs/src/lib/libcrypto/bn/bn_kron.c,v
retrieving revision 1.9
diff -u -p -r1.9 bn_kron.c
--- bn/bn_kron.c 20 Jun 2022 19:42:58 -0000 1.9
+++ bn/bn_kron.c 12 Jul 2022 10:41:09 -0000
@@ -55,9 +55,6 @@
#include "bn_lcl.h"
-/* The least significant word of a BIGNUM. */
-#define BN_lsw(n) (((n)->top == 0) ? (BN_ULONG) 0 : (n)->d[0])
-
/*
* Kronecker symbol, implemented according to Henri Cohen, "A Course in
* Computational Algebraic Number Theory", Algorithm 1.4.10.
Index: bn/bn_lcl.h
===================================================================
RCS file: /cvs/src/lib/libcrypto/bn/bn_lcl.h,v
retrieving revision 1.31
diff -u -p -r1.31 bn_lcl.h
--- bn/bn_lcl.h 14 Jan 2022 08:01:47 -0000 1.31
+++ bn/bn_lcl.h 12 Jul 2022 10:41:09 -0000
@@ -493,6 +493,9 @@ struct bn_gencb_st {
}
#endif /* !BN_LLONG */
+/* The least significant word of a BIGNUM. */
+#define BN_lsw(n) (((n)->top == 0) ? (BN_ULONG) 0 : (n)->d[0])
+
void bn_mul_normal(BN_ULONG *r, BN_ULONG *a, int na, BN_ULONG *b, int nb);
void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b);
void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b);
@@ -652,6 +655,11 @@ int BN_gcd_ct(BIGNUM *r, const BIGNUM *a
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_perfect_square(int *is_perfect_square, const BIGNUM *n, BN_CTX *ctx);
+
+int bn_is_prime_bpsw(int *is_prime, const BIGNUM *n, BN_CTX *in_ctx);
__END_HIDDEN_DECLS
#endif
Index: bn/bn_prime.c
===================================================================
RCS file: /cvs/src/lib/libcrypto/bn/bn_prime.c,v
retrieving revision 1.19
diff -u -p -r1.19 bn_prime.c
--- bn/bn_prime.c 18 Jun 2022 15:52:35 -0000 1.19
+++ bn/bn_prime.c 12 Jul 2022 10:41:09 -0000
@@ -255,16 +255,28 @@ BN_is_prime_ex(const BIGNUM *a, int chec
return BN_is_prime_fasttest_ex(a, checks, ctx_passed, 0, cb);
}
+#define LIBRESSL_HAS_BPSW
+
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
+ int is_prime;
+
+ /* XXX - tickle BN_GENCB in bn_is_prime_bpsw(). */
+ if (!bn_is_prime_bpsw(&is_prime, a, ctx_passed))
+ return -1;
+
+ return is_prime;
+#endif
if (BN_cmp(a, BN_value_one()) <= 0)
return 0;
