From 07cffb6a9623966ff329ea6953f11f66ab6431ed Mon Sep 17 00:00:00 2001 From: Mark Wooding Date: Sat, 3 Aug 2013 19:25:44 +0100 Subject: [PATCH] sshbn.c: Some number-theoretic utilities. The `kronecker' function computes the Kronecker symbol (a|n), which is an extension of the Jacobi symbol (itself an extension of the Legendre symbol) with three useful properties. * It's fairly fast to compute -- a little slower than a GCD using Euclid's algorithm. * It's well-defined on all inputs (unlike the Jacobi and Legendre symbols), so it's safe to call. * When n is prime, the result tells you whether a is a quadratic residue mod n. The code doesn't describe the algorithm in detail, but it basically follows from Gauss's work on quadratic residues in `Disquisitiones Arithmeticae'. The relevant theory is part of most introductory number theory texts; I like Shoup, `A Computational Introduction to Number Theory and Algebra', . The `modsqrt' function calculates square roots in a prime field. It takes a bit longer than a `modpow' call. The algorithm is basically the same as described in Menezes, van Oorschott and Vanstone, `Handbook of Applied Cryptography', ; but the explanation is based on Crandall and Pomerance, `Primes: A Computational Introduction'. The `kronecker' function is currently untested (except via `modsqrt'); `modsqrt' itself is tested via `testdata/bignum-fixed.txt' which will gather more non-randomized tests with time. Signed-off-by: Mark Wooding --- ssh.h | 2 + sshbn.c | 233 ++++++++++++++++++++++++++++++++++++++++++++++ testdata/bignum-fixed.txt | 10 ++ 3 files changed, 245 insertions(+) create mode 100644 testdata/bignum-fixed.txt diff --git a/ssh.h b/ssh.h index c604554b..1ac4e5e8 100644 --- a/ssh.h +++ b/ssh.h @@ -457,6 +457,8 @@ Bignum bigmuladd(Bignum a, Bignum b, Bignum addend); Bignum bigdiv(Bignum a, Bignum b); Bignum bigmod(Bignum a, Bignum b); Bignum modinv(Bignum number, Bignum modulus); +int kronecker(Bignum a, Bignum n); +Bignum modsqrt(Bignum a, Bignum p); Bignum bignum_bitmask(Bignum number); Bignum bignum_rshift(Bignum number, int shift); int bignum_cmp(Bignum a, Bignum b); diff --git a/sshbn.c b/sshbn.c index 2d81ee4c..0313071a 100644 --- a/sshbn.c +++ b/sshbn.c @@ -1674,6 +1674,203 @@ Bignum modinv(Bignum number, Bignum modulus) } /* + * Extract the largest power of 2 dividing x, storing it in p2, and returning + * the product of the remaining factors. + */ +static Bignum extract_p2(Bignum x, unsigned *p2) +{ + unsigned i, j, k, n; + Bignum y; + + /* If x is zero then the following won't work. And if x is odd then + * there's nothing very useful to do. + */ + if (!x[0] || (x[1] & 1)) { + *p2 = 0; + return copybn(x); + } + + /* Find the power of two. */ + for (i = 0; !x[i + 1]; i++); + for (j = 0; !((x[i + 1] >> j) & 1); j++); + *p2 = i*BIGNUM_INT_BITS + j; + + /* Work out how big the copy should be. */ + n = x[0] - i - 1; + if (x[x[0]] >> j) n++; + + /* Copy and shift down. */ + y = newbn(n); + for (k = 1; k <= n; k++) { + y[k] = x[k + i] >> j; + if (j && k < x[0]) y[k] |= x[k + i + 1] << (BIGNUM_INT_BITS - j); + } + + /* Done. */ + return y; +} + +/* + * Kronecker symbol (a|n). The result is always in { -1, 0, +1 }, and is + * zero if and only if a and n have a nontrivial common factor. Most + * usefully, if n is prime, this is the Legendre symbol, taking the value +1 + * if a is a quadratic residue mod n, and -1 otherwise; i.e., (a|p) == + * a^{(p-1)/2} (mod p). + */ +int kronecker(Bignum a, Bignum n) +{ + unsigned s, nn; + int r = +1; + Bignum t; + + /* Special case for n = 0. This is the same convention PARI uses, + * except that we can't represent negative numbers. + */ + if (bignum_cmp(n, Zero) == 0) { + if (bignum_cmp(a, One) == 0) return +1; + else return 0; + } + + /* Write n = 2^s t, with t odd. If s > 0 and a is even, then the answer + * is zero; otherwise throw in a factor of (-1)^s if a == 3 or 5 (mod 8). + * + * At this point, we have a copy of n, and must remember to free it when + * we're done. It's convenient to take a copy of a at the same time. + */ + a = copybn(a); + n = extract_p2(n, &s); + + if (s && (!a[0] || !(a[1] & 1))) { r = 0; goto done; } + else if ((s & 1) && ((a[1] & 7) == 3 || (a[1] & 7) == 5)) r = -r; + + /* If n is (now) a unit then we're done. */ + if (bignum_cmp(n, One) == 0) goto done; + + /* Reduce a modulo n before we go any further. */ + if (bignum_cmp(a, n) >= 0) { t = bigmod(a, n); freebn(a); a = t; } + + /* Main loop. */ + for (;;) { + if (bignum_cmp(a, Zero) == 0) { r = 0; goto done; } + + /* Strip out and handle powers of two from a. */ + t = extract_p2(a, &s); freebn(a); a = t; + nn = n[1] & 7; + if ((s & 1) && (nn == 3 || nn == 5)) r = -r; + if (bignum_cmp(a, One) == 0) break; + + /* Swap, applying quadratic reciprocity. */ + if ((nn & 3) == 3 && (a[1] & 3) == 3) r = -r; + t = bigmod(n, a); freebn(n); n = a; a = t; + } + + /* Tidy up: we're done. */ +done: + freebn(a); freebn(n); + return r; +} + +/* + * Modular square root. We must have p prime: extracting square roots modulo + * composites is equivalent to factoring (but we don't check: you'll just get + * the wrong answer). Returns NULL if x is not a quadratic residue mod p. + */ +Bignum modsqrt(Bignum x, Bignum p) +{ + Bignum xinv, b, c, r, t, z, X, mone; + unsigned i, j, s; + + /* If x is not a quadratic residue then we will not go to space today. */ + if (kronecker(x, p) != +1) return NULL; + + /* We need a quadratic nonresidue from somewhere. Exactly half of all + * units mod p are quadratic residues, but no efficient deterministic + * algorithm for finding one is known. So pick at random: we don't + * expect this to take long. + */ + z = newbn(p[0]); + do { + for (i = 1; i <= p[0]; i++) z[i] = rand(); + z[0] = p[0]; bn_restore_invariant(z); + } while (kronecker(z, p) != -1); + b = bigmod(z, p); freebn(z); + + /* We need to compute a few things before we really get started. */ + xinv = modinv(x, p); /* x^{-1} mod p */ + mone = bigsub(p, One); /* p - 1 == -1 (mod p) */ + t = extract_p2(mone, &s); /* 2^s t = p - 1 */ + c = modpow(b, t, p); /* b^t (mod p) */ + z = bigadd(t, One); freebn(t); t = z; /* (t + 1) */ + shift_right(t + 1, t[0], 1); if (!t[t[0]]) t[0]--; + r = modpow(x, t, p); /* x^{(t+1)/2} (mod p) */ + freebn(b); freebn(mone); freebn(t); + + /* OK, so how does this work anyway? + * + * We know that x^t is somewhere in the order-2^s subgroup of GF(p)^*; + * and g = c^{-1} is a generator for this subgroup (since we know that + * g^{2^{s-1}} = b^{(p-1)/2} = (b|p) = -1); so x^t = g^m for some m. In + * fact, we know that m is even because x is a square. Suppose we can + * determine m; then we know that x^t/g^m = 1, so x^{t+1}/c^m = x -- but + * both t + 1 and m are even, so x^{(t+1)/2}/g^{m/2} is a square root of + * x. + * + * Conveniently, finding the discrete log of an element X in a group of + * order 2^s is easy. Write X = g^m = g^{m_0+2k'}; then X^{2^{s-1}} = + * g^{m_0 2^{s-1}} c^{m' 2^s} = g^{m_0 2^{s-1}} is either -1 or +1, + * telling us that m_0 is 1 or 0 respectively. Then X/g^{m_0} = + * (g^2)^{m'} has order 2^{s-1} so we can continue inductively. What we + * end up with at the end is X/g^m. + * + * There are a few wrinkles. As we proceed through the induction, the + * generator for the subgroup will be c^{-2}, since we know that m is + * even. While we want the discrete log of X = x^t, we're actually going + * to keep track of r, which will eventually be x^{(t+1)/2}/g^{m/2} = + * x^{(t+1)/2} c^m, recovering X/g^m = r^2/x as we go. We don't actually + * form the discrete log explicitly, because the final result will + * actually be the square root we want. + */ + for (i = 1; i < s; i++) { + + /* Determine X. We could optimize this, only recomputing it when + * it's been invalidated, but that's fiddlier and this isn't + * performance critical. + */ + z = modmul(r, r, p); + X = modmul(z, xinv, p); + freebn(z); + + /* Determine X^{2^{s-1-i}}. */ + for (j = i + 1; j < s; j++) + z = modmul(X, X, p), freebn(X), X = z; + + /* Maybe accumulate a factor of c. */ + if (bignum_cmp(X, One) != 0) + z = modmul(r, c, p), freebn(r), r = z; + + /* Move on to the next smaller subgroup. */ + z = modmul(c, c, p), freebn(c), c = z; + freebn(X); + } + + /* Of course, there are two square roots of x. For predictability's sake + * we'll always return the one in [1..(p - 1)/2]. The other is, of + * course, p - r. + */ + z = bigsub(p, r); + if (bignum_cmp(r, z) < 0) + freebn(z); + else { + freebn(r); + r = z; + } + + /* We're done. */ + freebn(xinv); freebn(c); + return r; +} + +/* * Render a bignum into decimal. Return a malloced string holding * the decimal representation. */ @@ -1899,6 +2096,42 @@ int main(int argc, char **argv) freebn(modulus); freebn(expected); freebn(answer); + } else if (!strcmp(buf, "modsqrt")) { + Bignum x, p, expected, answer; + + if (ptrnum != 3) { + printf("%d: modsqrt with %d parameters, expected 3\n", line, ptrnum); + exit(1); + } + + x = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]); + p = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]); + expected = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]); + answer = modsqrt(x, p); + if (!answer) + answer = copybn(Zero); + + if (bignum_cmp(expected, answer) == 0) { + passes++; + } else { + char *xs = bignum_decimal(x); + char *ps = bignum_decimal(p); + char *qs = bignum_decimal(answer); + char *ws = bignum_decimal(expected); + + printf("%d: fail: sqrt(%s) mod %s gave %s expected %s\n", + line, xs, ps, qs, ws); + fails++; + + sfree(xs); + sfree(ps); + sfree(qs); + sfree(ws); + } + freebn(p); + freebn(x); + freebn(expected); + freebn(answer); } else { printf("%d: unrecognised test keyword: '%s'\n", line, buf); exit(1); diff --git a/testdata/bignum-fixed.txt b/testdata/bignum-fixed.txt new file mode 100644 index 00000000..baf15606 --- /dev/null +++ b/testdata/bignum-fixed.txt @@ -0,0 +1,10 @@ +modsqrt 1 3 1 +modsqrt 4 5 2 +modsqrt 87a9e12efc9fd544 b9d9d9fb4440f7c7 2d4e1588719986aa +modsqrt 712ec43313c96fe4 8d68905434b020eb 1cc7a93174d21d9c +modsqrt 2d146b877226105d 7c91fd10cc197bed 10f02ba313d1c57d +modsqrt 2fd8dff12ab86c17 8cf04fa2e760f1a5 2709aec642ddb824 +modsqrt 142c3b15dcac5e49 2af5fd8dd69a48b1 117e3c98c1fd7d6e +modsqrt 2eb75cd52c109ab0 3420d17b340dd381 87815684affc279 +modsqrt 35b331ec9bae9085 4ff9adae43a38f4f 26b2be1e41a0645b +modsqrt 92bdc828e44bf30 180089dd35d2dd17 951108874f6c30c -- 2.11.0