--- /dev/null
+/* $OpenBSD: bn_mod_sqrt.c,v 1.1 2023/04/11 10:08:44 tb Exp $ */
+
+/*
+ * Copyright (c) 2022 Theo Buehler <tb@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/err.h>
+
+#include "bn_local.h"
+
+/*
+ * Tonelli-Shanks according to H. Cohen "A Course in Computational Algebraic
+ * Number Theory", Section 1.5.1, Springer GTM volume 138, Berlin, 1996.
+ *
+ * Under the assumption that p is prime and a is a quadratic residue, we know:
+ *
+ * a^[(p-1)/2] = 1 (mod p). (*)
+ *
+ * To find a square root of a (mod p), we handle three cases of increasing
+ * complexity. In the first two cases, we can compute a square root using an
+ * explicit formula, thus avoiding the probabilistic nature of Tonelli-Shanks.
+ *
+ * 1. p = 3 (mod 4).
+ *
+ * Set n = (p+1)/4. Then 2n = 1 + (p-1)/2 and (*) shows that x = a^n (mod p)
+ * is a square root of a: x^2 = a^(2n) = a * a^[(p-1)/2] = a (mod p).
+ *
+ * 2. p = 5 (mod 8).
+ *
+ * This uses a simplification due to Atkin. By Theorem 1.4.7 and 1.4.9, the
+ * Kronecker symbol (2/p) evaluates to (-1)^[(p^2-1)/8]. From p = 5 (mod 8)
+ * we get (p^2-1)/8 = 1 (mod 2), so (2/p) = -1, and thus
+ *
+ * 2^[(p-1)/2] = -1 (mod p). (**)
+ *
+ * Set b = (2a)^[(p-5)/8]. With (p-1)/2 = 2 + (p-5)/2, (*) and (**) show
+ *
+ * i = 2 a b^2 is a square root of -1 (mod p).
+ *
+ * Indeed, i^2 = 2^2 a^2 b^4 = 2^[(p-1)/2] a^[(p-1)/2] = -1 (mod p). Because
+ * of (i-1)^2 = -2i (mod p) and i (-i) = 1 (mod p), a square root of a is
+ *
+ * x = a b (i-1)
+ *
+ * as x^2 = a^2 b^2 (-2i) = a (2 a b^2) (-i) = a (mod p).
+ *
+ * 3. p = 1 (mod 8).
+ *
+ * This is the Tonelli-Shanks algorithm. For a prime p, the multiplicative
+ * group of GF(p) is cyclic of order p - 1 = 2^s q, with odd q. Denote its
+ * 2-Sylow subgroup by S. It is cyclic of order 2^s. The squares in S have
+ * order dividing 2^(s-1). They are the even powers of any generator z of S.
+ * If a is a quadratic residue, 1 = a^[(p-1)/2] = (a^q)^[2^(s-1)], so b = a^q
+ * is a square in S. Therefore there is an integer k such that b z^(2k) = 1.
+ * Set x = a^[(q+1)/2] z^k, and find x^2 = a (mod p).
+ *
+ * The problem is thus reduced to finding a generator z of the 2-Sylow
+ * subgroup S of GF(p)* and finding k. An iterative constructions avoids
+ * the need for an explicit k, a generator is found by a randomized search.
+ *
+ * While we do not actually know that p is a prime number, we can still apply
+ * the formulas in cases 1 and 2 and verify that we have indeed found a square
+ * root of p. Similarly, in case 3, we can try to find a quadratic non-residue,
+ * which will fail for example if p is a square. The iterative construction
+ * may or may not find a candidate square root which we can then validate.
+ */
+
+/*
+ * Handle the cases where p is 2, p isn't odd or p is one. Since BN_mod_sqrt()
+ * can run on untrusted data, a primality check is too expensive. Also treat
+ * the obvious cases where a is 0 or 1.
+ */
+
+static int
+bn_mod_sqrt_trivial_cases(int *done, BIGNUM *out_sqrt, const BIGNUM *a,
+ const BIGNUM *p, BN_CTX *ctx)
+{
+ *done = 1;
+
+ if (BN_abs_is_word(p, 2))
+ return BN_set_word(out_sqrt, BN_is_odd(a));
+
+ if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ return 0;
+ }
+
+ if (BN_is_zero(a) || BN_is_one(a))
+ return BN_set_word(out_sqrt, BN_is_one(a));
+
+ *done = 0;
+
+ return 1;
+}
+
+/*
+ * Case 1. We know that (a/p) = 1 and that p = 3 (mod 4).
+ */
+
+static int
+bn_mod_sqrt_p_is_3_mod_4(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
+ BN_CTX *ctx)
+{
+ BIGNUM *n;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((n = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /* Calculate n = (|p| + 1) / 4. */
+ if (!BN_uadd(n, p, BN_value_one()))
+ goto err;
+ if (!BN_rshift(n, n, 2))
+ goto err;
+
+ /* By case 1 above, out_sqrt = a^n is a square root of a (mod p). */
+ if (!BN_mod_exp_ct(out_sqrt, a, n, p, ctx))
+ goto err;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Case 2. We know that (a/p) = 1 and that p = 5 (mod 8).
+ */
+
+static int
+bn_mod_sqrt_p_is_5_mod_8(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
+ BN_CTX *ctx)
+{
+ BIGNUM *b, *i, *n, *tmp;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((b = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((i = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((n = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((tmp = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /* Calculate n = (|p| - 5) / 8. Since p = 5 (mod 8), simply shift. */
+ if (!BN_rshift(n, p, 3))
+ goto err;
+ BN_set_negative(n, 0);
+
+ /* Compute tmp = 2a (mod p) for later use. */
+ if (!BN_mod_lshift1(tmp, a, p, ctx))
+ goto err;
+
+ /* Calculate b = (2a)^n (mod p). */
+ if (!BN_mod_exp_ct(b, tmp, n, p, ctx))
+ goto err;
+
+ /* Calculate i = 2 a b^2 (mod p). */
+ if (!BN_mod_sqr(i, b, p, ctx))
+ goto err;
+ if (!BN_mod_mul(i, tmp, i, p, ctx))
+ goto err;
+
+ /* A square root is out_sqrt = a b (i-1) (mod p). */
+ if (!BN_sub_word(i, 1))
+ goto err;
+ if (!BN_mod_mul(out_sqrt, a, b, p, ctx))
+ goto err;
+ if (!BN_mod_mul(out_sqrt, out_sqrt, i, p, ctx))
+ goto err;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Case 3. We know that (a/p) = 1 and that p = 1 (mod 8).
+ */
+
+/*
+ * Simple helper. To find a generator of the 2-Sylow subgroup of GF(p)*, we
+ * need to find a quadratic non-residue of p, i.e., n such that (n/p) = -1.
+ */
+
+static int
+bn_mod_sqrt_n_is_non_residue(int *is_non_residue, const BIGNUM *n,
+ const BIGNUM *p, BN_CTX *ctx)
+{
+ switch (BN_kronecker(n, p, ctx)) {
+ case -1:
+ *is_non_residue = 1;
+ return 1;
+ case 1:
+ *is_non_residue = 0;
+ return 1;
+ case 0:
+ /* n divides p, so ... */
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ return 0;
+ default:
+ return 0;
+ }
+}
+
+/*
+ * The following is the only non-deterministic part preparing Tonelli-Shanks.
+ *
+ * If we find n such that (n/p) = -1, then n^q (mod p) is a generator of the
+ * 2-Sylow subgroup of GF(p)*. To find such n, first try some small numbers,
+ * then random ones.
+ */
+
+static int
+bn_mod_sqrt_find_sylow_generator(BIGNUM *out_generator, const BIGNUM *p,
+ const BIGNUM *q, BN_CTX *ctx)
+{
+ BIGNUM *n, *p_abs, *thirty_two;
+ int i, is_non_residue;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((n = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((thirty_two = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((p_abs = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ for (i = 2; i < 32; i++) {
+ if (!BN_set_word(n, i))
+ goto err;
+ if (!bn_mod_sqrt_n_is_non_residue(&is_non_residue, n, p, ctx))
+ goto err;
+ if (is_non_residue)
+ goto found;
+ }
+
+ if (!BN_set_word(thirty_two, 32))
+ goto err;
+ if (!bn_copy(p_abs, p))
+ goto err;
+ BN_set_negative(p_abs, 0);
+
+ for (i = 0; i < 128; i++) {
+ if (!bn_rand_interval(n, thirty_two, p_abs))
+ goto err;
+ if (!bn_mod_sqrt_n_is_non_residue(&is_non_residue, n, p, ctx))
+ goto err;
+ if (is_non_residue)
+ goto found;
+ }
+
+ /*
+ * The probability to get here is < 2^(-128) for prime p. For squares
+ * it is easy: for p = 1369 = 37^2 this happens in ~3% of runs.
+ */
+
+ BNerror(BN_R_TOO_MANY_ITERATIONS);
+ goto err;
+
+ found:
+ /*
+ * If p is prime, n^q generates the 2-Sylow subgroup S of GF(p)*.
+ */
+
+ if (!BN_mod_exp_ct(out_generator, n, q, p, ctx))
+ goto err;
+
+ /* Sanity: p is not necessarily prime, so we could have found 0 or 1. */
+ if (BN_is_zero(out_generator) || BN_is_one(out_generator)) {
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ goto err;
+ }
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Initialization step for Tonelli-Shanks.
+ *
+ * In the end, b = a^q (mod p) and x = a^[(q+1)/2] (mod p). Cohen optimizes this
+ * to minimize taking powers of a. This is a bit confusing and distracting, so
+ * factor this into a separate function.
+ */
+
+static int
+bn_mod_sqrt_tonelli_shanks_initialize(BIGNUM *b, BIGNUM *x, const BIGNUM *a,
+ const BIGNUM *p, const BIGNUM *q, BN_CTX *ctx)
+{
+ BIGNUM *k;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((k = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /* k = (q-1)/2. Since q is odd, we can shift. */
+ if (!BN_rshift1(k, q))
+ goto err;
+
+ /* x = a^[(q-1)/2] (mod p). */
+ if (!BN_mod_exp_ct(x, a, k, p, ctx))
+ goto err;
+
+ /* b = ax^2 = a^q (mod p). */
+ if (!BN_mod_sqr(b, x, p, ctx))
+ goto err;
+ if (!BN_mod_mul(b, a, b, p, ctx))
+ goto err;
+
+ /* x = ax = a^[(q+1)/2] (mod p). */
+ if (!BN_mod_mul(x, a, x, p, ctx))
+ goto err;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Find smallest exponent m such that b^(2^m) = 1 (mod p). Assuming that a
+ * is a quadratic residue and p is a prime, we know that 1 <= m < r.
+ */
+
+static int
+bn_mod_sqrt_tonelli_shanks_find_exponent(int *out_exponent, const BIGNUM *b,
+ const BIGNUM *p, int r, BN_CTX *ctx)
+{
+ BIGNUM *x;
+ int m;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((x = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /*
+ * If r <= 1, the Tonelli-Shanks iteration should have terminated as
+ * r == 1 implies b == 1.
+ */
+ if (r <= 1) {
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ goto err;
+ }
+
+ /*
+ * Sanity check to ensure taking squares actually does something:
+ * If b is 1, the Tonelli-Shanks iteration should have terminated.
+ * If b is 0, something's very wrong, in particular p can't be prime.
+ */
+ if (BN_is_zero(b) || BN_is_one(b)) {
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ goto err;
+ }
+
+ if (!bn_copy(x, b))
+ goto err;
+
+ for (m = 1; m < r; m++) {
+ if (!BN_mod_sqr(x, x, p, ctx))
+ goto err;
+ if (BN_is_one(x))
+ break;
+ }
+
+ if (m >= r) {
+ /* This means a is not a quadratic residue. As (a/p) = 1, ... */
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ goto err;
+ }
+
+ *out_exponent = m;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * The update step. With the minimal m such that b^(2^m) = 1 (mod m),
+ * set t = y^[2^(r-m-1)] (mod p) and update x = xt, y = t^2, b = by.
+ * This preserves the loop invariants a b = x^2, y^[2^(r-1)] = -1 and
+ * b^[2^(r-1)] = 1.
+ */
+
+static int
+bn_mod_sqrt_tonelli_shanks_update(BIGNUM *b, BIGNUM *x, BIGNUM *y,
+ const BIGNUM *p, int m, int r, BN_CTX *ctx)
+{
+ BIGNUM *t;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((t = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /* t = y^[2^(r-m-1)] (mod p). */
+ if (!BN_set_bit(t, r - m - 1))
+ goto err;
+ if (!BN_mod_exp_ct(t, y, t, p, ctx))
+ goto err;
+
+ /* x = xt (mod p). */
+ if (!BN_mod_mul(x, x, t, p, ctx))
+ goto err;
+
+ /* y = t^2 = y^[2^(r-m)] (mod p). */
+ if (!BN_mod_sqr(y, t, p, ctx))
+ goto err;
+
+ /* b = by (mod p). */
+ if (!BN_mod_mul(b, b, y, p, ctx))
+ goto err;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+static int
+bn_mod_sqrt_p_is_1_mod_8(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
+ BN_CTX *ctx)
+{
+ BIGNUM *b, *q, *x, *y;
+ int e, m, r;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((b = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((q = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((x = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((y = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ /*
+ * Factor p - 1 = 2^e q with odd q. Since p = 1 (mod 8), we know e >= 3.
+ */
+
+ e = 1;
+ while (!BN_is_bit_set(p, e))
+ e++;
+ if (!BN_rshift(q, p, e))
+ goto err;
+
+ if (!bn_mod_sqrt_find_sylow_generator(y, p, q, ctx))
+ goto err;
+
+ /*
+ * Set b = a^q (mod p) and x = a^[(q+1)/2] (mod p).
+ */
+ if (!bn_mod_sqrt_tonelli_shanks_initialize(b, x, a, p, q, ctx))
+ goto err;
+
+ /*
+ * The Tonelli-Shanks iteration. Starting with r = e, the following loop
+ * invariants hold at the start of the loop.
+ *
+ * a b = x^2 (mod p)
+ * y^[2^(r-1)] = -1 (mod p)
+ * b^[2^(r-1)] = 1 (mod p)
+ *
+ * In particular, if b = 1 (mod p), x is a square root of a.
+ *
+ * Since p - 1 = 2^e q, we have 2^(e-1) q = (p - 1) / 2, so in the first
+ * iteration this follows from (a/p) = 1, (n/p) = -1, y = n^q, b = a^q.
+ *
+ * In subsequent iterations, t = y^[2^(r-m-1)], where m is the smallest
+ * m such that b^(2^m) = 1. With x = xt (mod p) and b = bt^2 (mod p) the
+ * first invariant is preserved, the second and third follow from
+ * y = t^2 (mod p) and r = m as well as the choice of m.
+ *
+ * Finally, r is strictly decreasing in each iteration. If p is prime,
+ * let S be the 2-Sylow subgroup of GF(p)*. We can prove the algorithm
+ * stops: Let S_r be the subgroup of S consisting of elements of order
+ * dividing 2^r. Then S_r = <y> and b is in S_(r-1). The S_r form a
+ * descending filtration of S and when r = 1, then b = 1.
+ */
+
+ for (r = e; r >= 1; r = m) {
+ /*
+ * Termination condition. If b == 1 then x is a square root.
+ */
+ if (BN_is_one(b))
+ goto done;
+
+ /* Find smallest exponent 1 <= m < r such that b^(2^m) == 1. */
+ if (!bn_mod_sqrt_tonelli_shanks_find_exponent(&m, b, p, r, ctx))
+ goto err;
+
+ /*
+ * With t = y^[2^(r-m-1)], update x = xt, y = t^2, b = by.
+ */
+ if (!bn_mod_sqrt_tonelli_shanks_update(b, x, y, p, m, r, ctx))
+ goto err;
+
+ /*
+ * Sanity check to make sure we don't loop indefinitely.
+ * bn_mod_sqrt_tonelli_shanks_find_exponent() ensures m < r.
+ */
+ if (r <= m)
+ goto err;
+ }
+
+ /*
+ * If p is prime, we should not get here.
+ */
+
+ BNerror(BN_R_NOT_A_SQUARE);
+ goto err;
+
+ done:
+ if (!bn_copy(out_sqrt, x))
+ goto err;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Choose the smaller of sqrt and |p| - sqrt.
+ */
+
+static int
+bn_mod_sqrt_normalize(BIGNUM *sqrt, const BIGNUM *p, BN_CTX *ctx)
+{
+ BIGNUM *x;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((x = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ if (!BN_lshift1(x, sqrt))
+ goto err;
+
+ if (BN_ucmp(x, p) > 0) {
+ if (!BN_usub(sqrt, p, sqrt))
+ goto err;
+ }
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+/*
+ * Verify that a = (sqrt_a)^2 (mod p). Requires that a is reduced (mod p).
+ */
+
+static int
+bn_mod_sqrt_verify(const BIGNUM *a, const BIGNUM *sqrt_a, const BIGNUM *p,
+ BN_CTX *ctx)
+{
+ BIGNUM *x;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((x = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ if (!BN_mod_sqr(x, sqrt_a, p, ctx))
+ goto err;
+
+ if (BN_cmp(x, a) != 0) {
+ BNerror(BN_R_NOT_A_SQUARE);
+ goto err;
+ }
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+static int
+bn_mod_sqrt_internal(BIGNUM *out_sqrt, const BIGNUM *a, const BIGNUM *p,
+ BN_CTX *ctx)
+{
+ BIGNUM *a_mod_p, *sqrt;
+ BN_ULONG lsw;
+ int done;
+ int kronecker;
+ int ret = 0;
+
+ BN_CTX_start(ctx);
+
+ if ((a_mod_p = BN_CTX_get(ctx)) == NULL)
+ goto err;
+ if ((sqrt = BN_CTX_get(ctx)) == NULL)
+ goto err;
+
+ if (!BN_nnmod(a_mod_p, a, p, ctx))
+ goto err;
+
+ if (!bn_mod_sqrt_trivial_cases(&done, sqrt, a_mod_p, p, ctx))
+ goto err;
+ if (done)
+ goto verify;
+
+ /*
+ * Make sure that the Kronecker symbol (a/p) == 1. In case p is prime
+ * this is equivalent to a having a square root (mod p). The cost of
+ * BN_kronecker() is O(log^2(n)). This is small compared to the cost
+ * O(log^4(n)) of Tonelli-Shanks.
+ */
+
+ if ((kronecker = BN_kronecker(a_mod_p, p, ctx)) == -2)
+ goto err;
+ if (kronecker <= 0) {
+ /* This error is only accurate if p is known to be a prime. */
+ BNerror(BN_R_NOT_A_SQUARE);
+ goto err;
+ }
+
+ lsw = BN_lsw(p);
+
+ if (lsw % 4 == 3) {
+ if (!bn_mod_sqrt_p_is_3_mod_4(sqrt, a_mod_p, p, ctx))
+ goto err;
+ } else if (lsw % 8 == 5) {
+ if (!bn_mod_sqrt_p_is_5_mod_8(sqrt, a_mod_p, p, ctx))
+ goto err;
+ } else if (lsw % 8 == 1) {
+ if (!bn_mod_sqrt_p_is_1_mod_8(sqrt, a_mod_p, p, ctx))
+ goto err;
+ } else {
+ /* Impossible to hit since the trivial cases ensure p is odd. */
+ BNerror(BN_R_P_IS_NOT_PRIME);
+ goto err;
+ }
+
+ if (!bn_mod_sqrt_normalize(sqrt, p, ctx))
+ goto err;
+
+ verify:
+ if (!bn_mod_sqrt_verify(a_mod_p, sqrt, p, ctx))
+ goto err;
+
+ if (!bn_copy(out_sqrt, sqrt))
+ goto err;
+
+ ret = 1;
+
+ err:
+ BN_CTX_end(ctx);
+
+ return ret;
+}
+
+BIGNUM *
+BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
+{
+ BIGNUM *out_sqrt;
+
+ if ((out_sqrt = in) == NULL)
+ out_sqrt = BN_new();
+ if (out_sqrt == NULL)
+ goto err;
+
+ if (!bn_mod_sqrt_internal(out_sqrt, a, p, ctx))
+ goto err;
+
+ return out_sqrt;
+
+ err:
+ if (out_sqrt != in)
+ BN_free(out_sqrt);
+
+ return NULL;
+}
+++ /dev/null
-/* $OpenBSD: bn_sqrt.c,v 1.16 2023/03/27 10:25:02 tb Exp $ */
-/* Written by Lenka Fibikova <fibikova@exp-math.uni-essen.de>
- * and Bodo Moeller for the OpenSSL project. */
-/* ====================================================================
- * Copyright (c) 1998-2000 The OpenSSL Project. All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in
- * the documentation and/or other materials provided with the
- * distribution.
- *
- * 3. All advertising materials mentioning features or use of this
- * software must display the following acknowledgment:
- * "This product includes software developed by the OpenSSL Project
- * for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
- *
- * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
- * endorse or promote products derived from this software without
- * prior written permission. For written permission, please contact
- * openssl-core@openssl.org.
- *
- * 5. Products derived from this software may not be called "OpenSSL"
- * nor may "OpenSSL" appear in their names without prior written
- * permission of the OpenSSL Project.
- *
- * 6. Redistributions of any form whatsoever must retain the following
- * acknowledgment:
- * "This product includes software developed by the OpenSSL Project
- * for use in the OpenSSL Toolkit (http://www.openssl.org/)"
- *
- * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
- * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
- * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE OpenSSL PROJECT OR
- * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
- * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
- * OF THE POSSIBILITY OF SUCH DAMAGE.
- * ====================================================================
- *
- * This product includes cryptographic software written by Eric Young
- * (eay@cryptsoft.com). This product includes software written by Tim
- * Hudson (tjh@cryptsoft.com).
- *
- */
-
-#include <openssl/err.h>
-
-#include "bn_local.h"
-
-/*
- * Returns 'ret' such that ret^2 == a (mod p), if it exists, using the
- * Tonelli-Shanks algorithm following Henri Cohen, "A Course in Computational
- * Algebraic Number Theory", algorithm 1.5.1, Springer, Berlin, 1996.
- *
- * Note: 'p' must be prime!
- */
-
-BIGNUM *
-BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
-{
- BIGNUM *ret = in;
- int err = 1;
- int r;
- BIGNUM *A, *b, *q, *t, *x, *y;
- int e, i, j;
-
- if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
- if (BN_abs_is_word(p, 2)) {
- if (ret == NULL)
- ret = BN_new();
- if (ret == NULL)
- goto end;
- if (!BN_set_word(ret, BN_is_bit_set(a, 0))) {
- if (ret != in)
- BN_free(ret);
- return NULL;
- }
- return ret;
- }
-
- BNerror(BN_R_P_IS_NOT_PRIME);
- return (NULL);
- }
-
- if (BN_is_zero(a) || BN_is_one(a)) {
- if (ret == NULL)
- ret = BN_new();
- if (ret == NULL)
- goto end;
- if (!BN_set_word(ret, BN_is_one(a))) {
- if (ret != in)
- BN_free(ret);
- return NULL;
- }
- return ret;
- }
-
- BN_CTX_start(ctx);
- if ((A = BN_CTX_get(ctx)) == NULL)
- goto end;
- if ((b = BN_CTX_get(ctx)) == NULL)
- goto end;
- if ((q = BN_CTX_get(ctx)) == NULL)
- goto end;
- if ((t = BN_CTX_get(ctx)) == NULL)
- goto end;
- if ((x = BN_CTX_get(ctx)) == NULL)
- goto end;
- if ((y = BN_CTX_get(ctx)) == NULL)
- goto end;
-
- if (ret == NULL)
- ret = BN_new();
- if (ret == NULL)
- goto end;
-
- /* A = a mod p */
- if (!BN_nnmod(A, a, p, ctx))
- goto end;
-
- /* now write |p| - 1 as 2^e*q where q is odd */
- e = 1;
- while (!BN_is_bit_set(p, e))
- e++;
- /* we'll set q later (if needed) */
-
- if (e == 1) {
- /* The easy case: (|p|-1)/2 is odd, so 2 has an inverse
- * modulo (|p|-1)/2, and square roots can be computed
- * directly by modular exponentiation.
- * We have
- * 2 * (|p|+1)/4 == 1 (mod (|p|-1)/2),
- * so we can use exponent (|p|+1)/4, i.e. (|p|-3)/4 + 1.
- */
- if (!BN_rshift(q, p, 2))
- goto end;
- q->neg = 0;
- if (!BN_add_word(q, 1))
- goto end;
- if (!BN_mod_exp_ct(ret, A, q, p, ctx))
- goto end;
- err = 0;
- goto vrfy;
- }
-
- if (e == 2) {
- /* |p| == 5 (mod 8)
- *
- * In this case 2 is always a non-square since
- * Legendre(2,p) = (-1)^((p^2-1)/8) for any odd prime.
- * So if a really is a square, then 2*a is a non-square.
- * Thus for
- * b := (2*a)^((|p|-5)/8),
- * i := (2*a)*b^2
- * we have
- * i^2 = (2*a)^((1 + (|p|-5)/4)*2)
- * = (2*a)^((p-1)/2)
- * = -1;
- * so if we set
- * x := a*b*(i-1),
- * then
- * x^2 = a^2 * b^2 * (i^2 - 2*i + 1)
- * = a^2 * b^2 * (-2*i)
- * = a*(-i)*(2*a*b^2)
- * = a*(-i)*i
- * = a.
- *
- * (This is due to A.O.L. Atkin,
- * <URL: http://listserv.nodak.edu/scripts/wa.exe?A2=ind9211&L=nmbrthry&O=T&P=562>,
- * November 1992.)
- */
-
- /* t := 2*a */
- if (!BN_mod_lshift1_quick(t, A, p))
- goto end;
-
- /* b := (2*a)^((|p|-5)/8) */
- if (!BN_rshift(q, p, 3))
- goto end;
- q->neg = 0;
- if (!BN_mod_exp_ct(b, t, q, p, ctx))
- goto end;
-
- /* y := b^2 */
- if (!BN_mod_sqr(y, b, p, ctx))
- goto end;
-
- /* t := (2*a)*b^2 - 1*/
- if (!BN_mod_mul(t, t, y, p, ctx))
- goto end;
- if (!BN_sub_word(t, 1))
- goto end;
-
- /* x = a*b*t */
- if (!BN_mod_mul(x, A, b, p, ctx))
- goto end;
- if (!BN_mod_mul(x, x, t, p, ctx))
- goto end;
-
- if (!bn_copy(ret, x))
- goto end;
- err = 0;
- goto vrfy;
- }
-
- /* e > 2, so we really have to use the Tonelli/Shanks algorithm.
- * First, find some y that is not a square. */
- if (!bn_copy(q, p)) /* use 'q' as temp */
- goto end;
- q->neg = 0;
- i = 2;
- do {
- /* For efficiency, try small numbers first;
- * if this fails, try random numbers.
- */
- if (i < 22) {
- if (!BN_set_word(y, i))
- goto end;
- } else {
- if (!BN_pseudo_rand(y, BN_num_bits(p), 0, 0))
- goto end;
- if (BN_ucmp(y, p) >= 0) {
- if (p->neg) {
- if (!BN_add(y, y, p))
- goto end;
- } else {
- if (!BN_sub(y, y, p))
- goto end;
- }
- }
- /* now 0 <= y < |p| */
- if (BN_is_zero(y))
- if (!BN_set_word(y, i))
- goto end;
- }
-
- r = BN_kronecker(y, q, ctx); /* here 'q' is |p| */
- if (r < -1)
- goto end;
- if (r == 0) {
- /* m divides p */
- BNerror(BN_R_P_IS_NOT_PRIME);
- goto end;
- }
- } while (r == 1 && ++i < 82);
-
- if (r != -1) {
- /* Many rounds and still no non-square -- this is more likely
- * a bug than just bad luck.
- * Even if p is not prime, we should have found some y
- * such that r == -1.
- */
- BNerror(BN_R_TOO_MANY_ITERATIONS);
- goto end;
- }
-
- /* Here's our actual 'q': */
- if (!BN_rshift(q, q, e))
- goto end;
-
- /* Now that we have some non-square, we can find an element
- * of order 2^e by computing its q'th power. */
- if (!BN_mod_exp_ct(y, y, q, p, ctx))
- goto end;
- if (BN_is_one(y)) {
- BNerror(BN_R_P_IS_NOT_PRIME);
- goto end;
- }
-
- /* Now we know that (if p is indeed prime) there is an integer
- * k, 0 <= k < 2^e, such that
- *
- * a^q * y^k == 1 (mod p).
- *
- * As a^q is a square and y is not, k must be even.
- * q+1 is even, too, so there is an element
- *
- * X := a^((q+1)/2) * y^(k/2),
- *
- * and it satisfies
- *
- * X^2 = a^q * a * y^k
- * = a,
- *
- * so it is the square root that we are looking for.
- */
-
- /* t := (q-1)/2 (note that q is odd) */
- if (!BN_rshift1(t, q))
- goto end;
-
- /* x := a^((q-1)/2) */
- if (BN_is_zero(t)) { /* special case: p = 2^e + 1 */
- if (!BN_nnmod(t, A, p, ctx))
- goto end;
- if (BN_is_zero(t)) {
- /* special case: a == 0 (mod p) */
- BN_zero(ret);
- err = 0;
- goto end;
- } else if (!BN_one(x))
- goto end;
- } else {
- if (!BN_mod_exp_ct(x, A, t, p, ctx))
- goto end;
- if (BN_is_zero(x)) {
- /* special case: a == 0 (mod p) */
- BN_zero(ret);
- err = 0;
- goto end;
- }
- }
-
- /* b := a*x^2 (= a^q) */
- if (!BN_mod_sqr(b, x, p, ctx))
- goto end;
- if (!BN_mod_mul(b, b, A, p, ctx))
- goto end;
-
- /* x := a*x (= a^((q+1)/2)) */
- if (!BN_mod_mul(x, x, A, p, ctx))
- goto end;
-
- while (1) {
- /* Now b is a^q * y^k for some even k (0 <= k < 2^E
- * where E refers to the original value of e, which we
- * don't keep in a variable), and x is a^((q+1)/2) * y^(k/2).
- *
- * We have a*b = x^2,
- * y^2^(e-1) = -1,
- * b^2^(e-1) = 1.
- */
-
- if (BN_is_one(b)) {
- if (!bn_copy(ret, x))
- goto end;
- err = 0;
- goto vrfy;
- }
-
- /* Find the smallest i with 0 < i < e such that b^(2^i) = 1. */
- for (i = 1; i < e; i++) {
- if (i == 1) {
- if (!BN_mod_sqr(t, b, p, ctx))
- goto end;
- } else {
- if (!BN_mod_sqr(t, t, p, ctx))
- goto end;
- }
- if (BN_is_one(t))
- break;
- }
- if (i >= e) {
- BNerror(BN_R_NOT_A_SQUARE);
- goto end;
- }
-
- /* t := y^2^(e - i - 1) */
- if (!bn_copy(t, y))
- goto end;
- for (j = e - i - 1; j > 0; j--) {
- if (!BN_mod_sqr(t, t, p, ctx))
- goto end;
- }
- if (!BN_mod_mul(y, t, t, p, ctx))
- goto end;
- if (!BN_mod_mul(x, x, t, p, ctx))
- goto end;
- if (!BN_mod_mul(b, b, y, p, ctx))
- goto end;
- e = i;
- }
-
-vrfy:
- if (!err) {
- /* verify the result -- the input might have been not a square
- * (test added in 0.9.8) */
-
- if (!BN_mod_sqr(x, ret, p, ctx))
- err = 1;
-
- if (!err && 0 != BN_cmp(x, A)) {
- BNerror(BN_R_NOT_A_SQUARE);
- err = 1;
- }
- }
-
-end:
- if (err) {
- if (ret != NULL && ret != in) {
- BN_free(ret);
- }
- ret = NULL;
- }
- BN_CTX_end(ctx);
- return ret;
-}
projects / openbsd / commitdiff