From 511819f0f73c5f37d3f7c033f01d29ef88bfa9cb Mon Sep 17 00:00:00 2001 From: Rob Austein Date: Fri, 21 Aug 2015 14:21:10 -0400 Subject: Updated point doubling and addition to use algorithms from the hyperelliptic.org formula database. Compiles, still not tested. --- ecdsa.c | 207 ++++++++++++++++++++++++++++++---------------------------------- 1 file changed, 97 insertions(+), 110 deletions(-) diff --git a/ecdsa.c b/ecdsa.c index 49893de..8e715f4 100644 --- a/ecdsa.c +++ b/ecdsa.c @@ -253,10 +253,16 @@ static inline int point_equal(const ec_point_t * const P, * All of these are operations in the field underlying the specified * curve, and assume that operands are already in Montgomery form. * - * Several of these are written a bit oddly, in an attempt to make - * them run in constant time. Be warned that an optimizing compiler - * may be clever enough to defeat this. In the long run, the real - * solution is probably to perform these field operations in Verilog. + * The ff_add() and ff_sub() are written a bit oddly, in an attempt to + * make them run in constant time. An optimizing compiler may be + * clever enough to defeat this. In the long run, we probably want to + * perform these field operations in Verilog anyway. + * + * We might be able to squeeze a bit more speed out of the point + * arithmetic by making using fp_mul_2d() when multiplying by a power + * of two. Skipping for now as a premature optimization, but if we do + * need this, it'd probably be simplest to add a ff_dbl() function + * which handles overflow in the same way that ff_add() does. */ static inline void ff_add(const ecdsa_curve_t * const curve, @@ -266,9 +272,12 @@ static inline void ff_add(const ecdsa_curve_t * const curve, { fp_int t[2][1]; memset(t, 0, sizeof(t)); + fp_add(unconst_fp_int(a), unconst_fp_int(b), t[0]); fp_sub(t[0], unconst_fp_int(curve->q), t[1]); - fp_copy(t[fp_cmp(t[0], unconst_fp_int(curve->q)) != FP_LT], c); + + fp_copy(t[fp_cmp_d(t[1], 0) != FP_LT], c); + memset(t, 0, sizeof(t)); } @@ -279,21 +288,12 @@ static inline void ff_sub(const ecdsa_curve_t * const curve, { fp_int t[2][1]; memset(t, 0, sizeof(t)); + fp_sub(unconst_fp_int(a), unconst_fp_int(b), t[0]); fp_add(t[0], unconst_fp_int(curve->q), t[1]); - fp_copy(t[fp_cmp_d(c, 0) == FP_LT], c); - memset(t, 0, sizeof(t)); -} -static inline void ff_div2(const ecdsa_curve_t * const curve, - const fp_int * const a, - fp_int *b) -{ - fp_int t[2][1]; - memset(t, 0, sizeof(t)); - fp_copy(unconst_fp_int(a), t[0]); - fp_add(t[0], unconst_fp_int(curve->q), t[1]); - fp_div_2(t[fp_isodd(unconst_fp_int(a))], b); + fp_copy(t[fp_cmp_d(t[0], 0) == FP_LT], c); + memset(t, 0, sizeof(t)); } @@ -314,37 +314,14 @@ static inline void ff_sqr(const ecdsa_curve_t * const curve, fp_montgomery_reduce(b, unconst_fp_int(curve->q), curve->rho); } -#warning Change point arithmetic algorithms? -/* - * The point doubling and addition algorithms we use here are from - * libtomcrypt. The formula database at hyperelliptic.org lists - * faster algorithms satisfying the same preconditions, perhaps we - * should use those instead? - * - * Labels in the following refer to entries on the page: - * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html - * - * The libtomcrypt doubling algorithm looks like a trivial variation - * on dbl-2004-hmv. We might want to use dbl-2001-b instead. - * - * The libtomcrypt addition algorithm doesn't match up exactly with - * any listed algorithm, but I suspect it's a variation on - * add-1998-cmo-2. We might want to use add-2007-bl instead. - * - * There are faster algorithms listed, but all of them appear to - * require whacking one or both points back into affine - * representation, which has its own costs, so, at least for now, it'd - * probably be best to stick with algorithms that don't require this. - */ - /** * Double an EC point. * @param P The point to double * @param R [out] The destination of the double * @param curve The curve parameters structure * - * Algorithm is a minor variation on algorithm 3.21 from Guide to - * Elliptic Curve Cryptography. + * Algorithm is dbl-2001-b from + * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html */ static inline void point_double(const ec_point_t * const P, @@ -353,36 +330,41 @@ static inline void point_double(const ec_point_t * const P, { assert(P != NULL && R != NULL && curve != NULL); - fp_int t1[1]; fp_init(t1); - fp_int t2[1]; fp_init(t2); - - if (P != R) - *R = *P; - - ff_sqr (curve, R->z, t1); /* t1 = Pz ** 2 */ - ff_sub (curve, R->x, t1, t2); /* t2 = Px - Pz ** 2 */ - ff_add (curve, R->x, t1, t1); /* t1 = Px + Pz ** 2 */ - ff_mul (curve, t1, t2, t2); /* t2 = 1 * (Px - Pz ** 2) * (Px + Pz ** 2) */ - ff_add (curve, t2, t2, t1); /* t1 = 2 * (Px - Pz ** 2) * (Px + Pz ** 2) */ - ff_add (curve, t1, t2, t1); /* t1 = 3 * (Px - Pz ** 2) * (Px + Pz ** 2) = A */ - - ff_add (curve, R->y, R->y, R->y); /* Ry = 2 * Py = B */ - ff_mul (curve, R->z, R->y, R->z); /* Rz = B * Pz */ - - ff_sqr (curve, R->y, R->y); /* Ry = B ** 2 = C */ - ff_sqr (curve, R->y, t2); /* t2 = C ** 2 */ - ff_div2 (curve, t2, t2); /* t2 = C ** 2 / 2 */ - ff_mul (curve, R->y, R->x, R->y); /* Ry = C * Px = D */ - - ff_sqr (curve, t1, R->x); /* Rx = A ** 2 */ - ff_sub (curve, R->x, R->y, R->x); /* Rx = A ** 2 - D */ - ff_sub (curve, R->x, R->y, R->x); /* Rx = A ** 2 - 2 * D */ - - ff_sub (curve, R->y, R->x, R->y); /* Ry = D - Rx */ - ff_mul (curve, R->y, t1, R->y); /* Ry = (D - Rx) * A */ - ff_sub (curve, R->y, t2, R->y); /* Ry = (D - Rx) * A - C ** 2 / 2 */ - - fp_zero(t1); fp_zero(t2); + fp_int alpha[1], beta[1], gamma[1], delta[1], t1[1], t2[1]; + + fp_init(alpha); fp_init(beta); fp_init(gamma); fp_init(delta); fp_init(t1); fp_init(t2); + + ff_sqr (curve, P->z, delta); /* delta = Pz ** 2 */ + ff_sqr (curve, P->y, gamma); /* gamma = Py ** 2 */ + ff_mul (curve, P->x, gamma, beta); /* beta = Px * gamma */ + ff_sub (curve, P->x, delta, t1); /* alpha = 3 * (Px - delta) * (Px + delta) */ + ff_add (curve, P->x, delta, t2); + ff_mul (curve, t1, t2, t1); + ff_add (curve, t1, t1, t2); + ff_add (curve, t1, t2, alpha); + + ff_sqr (curve, alpha, t1); /* Rx = (alpha ** 2) - (8 * beta) */ + ff_add (curve, beta, beta, t2); + ff_add (curve, t2, t2, t2); + ff_add (curve, t2, t2, t2); + ff_sub (curve, t1, t2, R->x); + + ff_add (curve, P->y, P->z, t1); /* Rz = ((Py + Pz) ** 2) - gamma - delta */ + ff_sqr (curve, t1, t1); + ff_sub (curve, t1, gamma, t1); + ff_sub (curve, t1, delta, R->z); + + ff_add (curve, beta, beta, t1); /* Ry = (((4 * beta) - Rx) * alpha) - (8 * (gamma ** 2)) */ + ff_add (curve, t1, t1, t1); + ff_sub (curve, t1, R->x, t1); + ff_mul (curve, t1, alpha, t1); + ff_sqr (curve, gamma, t2); + ff_add (curve, t2, t2, t2); + ff_add (curve, t2, t2, t2); + ff_add (curve, t2, t2, t2); + ff_sub (curve, t1, t2, R->y); + + fp_zero(alpha); fp_zero(beta); fp_zero(gamma); fp_zero(delta); fp_zero(t1); fp_zero(t2); } /** @@ -391,7 +373,10 @@ static inline void point_double(const ec_point_t * const P, * @param Q The point to add * @param R [out] The destination of the double * @param curve The curve parameters structure -*/ + * + * Algorithm is add-2007-bl from + * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html + */ static inline void point_add(const ec_point_t * const P, const ec_point_t * const Q, @@ -403,54 +388,56 @@ static inline void point_add(const ec_point_t * const P, if (point_equal(P, Q, curve)) return point_double(P, R, curve); - fp_int t1[1]; fp_init(t1); - fp_int t2[1]; fp_init(t2); + fp_int Z1Z1[1], Z2Z2[1], U1[1], U2[1], S1[1], S2[1], H[1], I[1], J[1], r[1], V[1], t[1]; - if (P != R) - *R = *P; + fp_init(Z1Z1), fp_init(Z2Z2), fp_init(U1), fp_init(U2), fp_init(S1), fp_init(S2); + fp_init(H), fp_init(I), fp_init(J), fp_init(r), fp_init(V), fp_init(t); - /* - * Operations marked {@} are no-ops when Q.z == 1, but probably - * don't save us enough in the long run for optimizing them out to - * be worth even a low-probability risk of a timing channel attack. - */ + ff_sqr (curve, P->z, Z1Z1); /* Z1Z1 = Pz ** 2 */ + + ff_sqr (curve, Q->z, Z2Z2); /* Z2Z1 = Qz ** 2 */ + + ff_mul (curve, P->x, Z2Z2, U1); /* U1 = Px * Z2Z2 */ + + ff_mul (curve, Q->x, Z1Z1, U2); /* U2 = Qx * Z1Z1 */ + + ff_mul (curve, Q->z, Z2Z2, S1); /* S1 = Py * (Qz ** 3) */ + ff_mul (curve, P->y, S1, S1); + + ff_mul (curve, P->z, Z1Z1, S2); /* S2 = Qy * (Pz ** 3) */ + ff_mul (curve, Q->y, S2, S2); - ff_sqr (curve, Q->z, t1); /* t1 = z' ** 2 {@} */ - ff_mul (curve, t1, R->x, R->x); /* x = x * z' ** 2 {@} */ - ff_mul (curve, Q->z, t1, t1); /* t1 = z' ** 3 {@} */ - ff_mul (curve, t1, R->y, R->y); /* y = y * z' ** 3 {@} */ + ff_sub (curve, U2, U1, H); /* H = U2 - U1 */ - ff_sqr (curve, R->z, t1); /* t1 = z * z */ - ff_mul (curve, Q->x, t1, t2); /* t2 = x' * t1 */ - ff_mul (curve, R->z, t1, t1); /* t1 = z * t1 */ - ff_mul (curve, Q->y, t1, t1); /* t1 = y' * t1 */ + ff_add (curve, H, H, I); /* I = (2 * H) ** 2 */ + ff_sqr (curve, I, I); - ff_sub (curve, R->y, t1, R->y); /* y = y - t1 */ - ff_add (curve, t1, t1, t1); /* t1 = 2 * t1 */ - ff_add (curve, t1, R->y, t1); /* t1 = y + t1 */ - ff_sub (curve, R->x, t2, R->x); /* x = x - t2 */ - ff_add (curve, t2, t2, t2); /* t2 = 2 * t2 */ - ff_add (curve, t2, R->x, t2); /* t2 = x + t2 */ + ff_mul (curve, H, I, J); /* J = H * I */ - ff_mul (curve, R->z, Q->z, R->z); /* z = z * z' {@} */ + ff_sub (curve, S2, S1, r); /* r = 2 * (S2 - S1) */ + ff_add (curve, r, r, r); - ff_mul (curve, R->z, R->x, R->z); /* z = z * x */ + ff_mul (curve, U1, I, V); /* V = U1 * I */ - ff_mul (curve, t1, R->x, t1); /* t1 = t1 * x */ - ff_sqr (curve, R->x, R->x); /* x = x * x */ - ff_mul (curve, t2, R->x, t2); /* t2 = t2 * x */ - ff_mul (curve, t1, R->x, t1); /* t1 = t1 * x */ + ff_sqr (curve, r, R->x); /* Rx = (r ** 2) - J - (2 * V) */ + ff_sub (curve, R->x, J, R->x); + ff_sub (curve, R->x, V, R->x); + ff_sub (curve, R->x, V, R->x); - ff_sqr (curve, R->y, R->x); /* x = y * y */ - ff_sub (curve, R->x, t2, R->x); /* x = x - t2 */ + ff_sub (curve, V, R->x, R->y); /* Ry = (r * (V - Rx)) - (2 * S1 * J) */ + ff_mul (curve, r, R->y, R->y); + ff_mul (curve, S1, J, t); + ff_sub (curve, R->y, t, R->y); + ff_sub (curve, R->y, t, R->y); - ff_sub (curve, t2, R->x, t2); /* t2 = t2 - x */ - ff_sub (curve, t2, R->x, t2); /* t2 = t2 - x */ - ff_mul (curve, t2, R->y, t2); /* t2 = t2 * y */ - ff_sub (curve, t2, t1, R->y); /* y = t2 - t1 */ - ff_div2 (curve, R->y, R->y); /* y = y / 2 */ + ff_add (curve, P->z, Q->z, R->z); /* Rz = (((Pz + Qz) ** 2) - Z1Z1 - Z2Z2) * H */ + ff_sqr (curve, R->z, R->z); + ff_sub (curve, R->z, Z1Z1, R->z); + ff_sub (curve, R->z, Z2Z2, R->z); + ff_mul (curve, R->z, H, R->z); - fp_zero(t1); fp_zero(t2); + fp_zero(Z1Z1), fp_zero(Z2Z2), fp_zero(U1), fp_zero(U2), fp_zero(S1), fp_zero(S2); + fp_zero(H), fp_zero(I), fp_zero(J), fp_zero(r), fp_zero(V), fp_zero(t); } /** -- cgit v1.2.3