EC GFp ladder

This commit leverages the Montgomery ladder scaffold introduced in #6690
(alongside a specialized Lopez-Dahab ladder for binary curves) to
provide a specialized differential addition-and-double implementation to
speedup prime curves, while keeping all the features of
`ec_scalar_mul_ladder` against SCA attacks.

The arithmetic in ladder_pre, ladder_step and ladder_post is auto
generated with tooling, from the following formulae:

- `ladder_pre`: Formula 3 for doubling from Izu-Takagi "A fast parallel
  elliptic curve multiplication resistant against side channel attacks",
  as described at
  https://hyperelliptic.org/EFD/g1p/auto-shortw-xz.html#doubling-dbl-2002-it-2
- `ladder_step`: differential addition-and-doubling Eq. (8) and (10)
  from Izu-Takagi "A fast parallel elliptic curve multiplication
  resistant against side channel attacks", as described at
  https://hyperelliptic.org/EFD/g1p/auto-shortw-xz.html#ladder-ladd-2002-it-3
- `ladder_post`: y-coordinate recovery using Eq. (8) from Brier-Joye
  "Weierstrass Elliptic Curves and Side-Channel Attacks", modified to
  work in projective coordinates.
Co-authored-by: NNicola Tuveri <nic.tuv@gmail.com>
Reviewed-by: NAndy Polyakov <appro@openssl.org>
Reviewed-by: NRich Salz <rsalz@openssl.org>
(Merged from https://github.com/openssl/openssl/pull/6772)

CHANGES

@@ -9,6 +9,15 @@
 
  Changes between 1.1.0h and 1.1.1 [xx XXX xxxx]
 
+  *) Use the new ec_scalar_mul_ladder scaffold to implement a specialized ladder
+     step for prime curves. The new implementation is based on formulae from
+     differential addition-and-doubling in homogeneous projective coordinates
+     from Izu-Takagi "A fast parallel elliptic curve multiplication resistant
+     against side channel attacks" and Brier-Joye "Weierstrass Elliptic Curves
+     and Side-Channel Attacks" Eq. (8) for y-coordinate recovery, modified
+     to work in projective coordinates.
+     [Billy Bob Brumley, Nicola Tuveri]
+
   *) Change generating and checking of primes so that the error rate of not
      being prime depends on the intended use based on the size of the input.
      For larger primes this will result in more rounds of Miller-Rabin.
@@ -30,7 +39,7 @@
      [Andy Polyakov]
 
   *) Use the new ec_scalar_mul_ladder scaffold to implement a specialized ladder
-     step for binary curves. The new implementation is based on formulas from
+     step for binary curves. The new implementation is based on formulae from
      differential addition-and-doubling in mixed Lopez-Dahab projective
      coordinates, modified to independently blind the operands.
      [Billy Bob Brumley, Sohaib ul Hassan, Nicola Tuveri]

crypto/ec/ec_lcl.h

@@ -301,7 +301,6 @@ struct ec_point_st {
                                  * special case */
 };
 
-
 static ossl_inline int ec_point_is_compat(const EC_POINT *point,
                                           const EC_GROUP *group)
 {
@@ -314,7 +313,6 @@ static ossl_inline int ec_point_is_compat(const EC_POINT *point,
     return 1;
 }
 
-
 NISTP224_PRE_COMP *EC_nistp224_pre_comp_dup(NISTP224_PRE_COMP *);
 NISTP256_PRE_COMP *EC_nistp256_pre_comp_dup(NISTP256_PRE_COMP *);
 NISTP521_PRE_COMP *EC_nistp521_pre_comp_dup(NISTP521_PRE_COMP *);
@@ -394,7 +392,16 @@ int ec_GFp_simple_field_mul(const EC_GROUP *, BIGNUM *r, const BIGNUM *a,
 int ec_GFp_simple_field_sqr(const EC_GROUP *, BIGNUM *r, const BIGNUM *a,
                             BN_CTX *);
 int ec_GFp_simple_blind_coordinates(const EC_GROUP *group, EC_POINT *p,
-                                   BN_CTX *ctx);
+                                    BN_CTX *ctx);
+int ec_GFp_simple_ladder_pre(const EC_GROUP *group,
+                             EC_POINT *r, EC_POINT *s,
+                             EC_POINT *p, BN_CTX *ctx);
+int ec_GFp_simple_ladder_step(const EC_GROUP *group,
+                              EC_POINT *r, EC_POINT *s,
+                              EC_POINT *p, BN_CTX *ctx);
+int ec_GFp_simple_ladder_post(const EC_GROUP *group,
+                              EC_POINT *r, EC_POINT *s,
+                              EC_POINT *p, BN_CTX *ctx);
 
 /* method functions in ecp_mont.c */
 int ec_GFp_mont_group_init(EC_GROUP *);

crypto/ec/ecp_mont.c

@@ -64,9 +64,9 @@ const EC_METHOD *EC_GFp_mont_method(void)
         ecdh_simple_compute_key,
         0, /* field_inverse_mod_ord */
         ec_GFp_simple_blind_coordinates,
-        0, /* ladder_pre */
-        0, /* ladder_step */
-        0  /* ladder_post */
+        ec_GFp_simple_ladder_pre,
+        ec_GFp_simple_ladder_step,
+        ec_GFp_simple_ladder_post
     };
 
     return &ret;

crypto/ec/ecp_nist.c

@@ -66,9 +66,9 @@ const EC_METHOD *EC_GFp_nist_method(void)
         ecdh_simple_compute_key,
         0, /* field_inverse_mod_ord */
         ec_GFp_simple_blind_coordinates,
-        0, /* ladder_pre */
-        0, /* ladder_step */
-        0  /* ladder_post */
+        ec_GFp_simple_ladder_pre,
+        ec_GFp_simple_ladder_step,
+        ec_GFp_simple_ladder_post
     };
 
     return &ret;

crypto/ec/ecp_smpl.c

@@ -65,9 +65,9 @@ const EC_METHOD *EC_GFp_simple_method(void)
         ecdh_simple_compute_key,
         0, /* field_inverse_mod_ord */
         ec_GFp_simple_blind_coordinates,
-        0, /* ladder_pre */
-        0, /* ladder_step */
-        0 /* ladder_post */
+        ec_GFp_simple_ladder_pre,
+        ec_GFp_simple_ladder_step,
+        ec_GFp_simple_ladder_post
     };
 
     return &ret;
@@ -1418,6 +1418,224 @@ int ec_GFp_simple_blind_coordinates(const EC_GROUP *group, EC_POINT *p,
     ret = 1;
 
  err:
-     BN_CTX_end(ctx);
-     return ret;
+    BN_CTX_end(ctx);
+    return ret;
+}
+
+/*-
+ * Set s := p, r := 2p.
+ *
+ * For doubling we use Formula 3 from Izu-Takagi "A fast parallel elliptic curve
+ * multiplication resistant against side channel attacks" appendix, as described
+ * at
+ * https://hyperelliptic.org/EFD/g1p/auto-shortw-xz.html#doubling-dbl-2002-it-2
+ *
+ * The input point p will be in randomized Jacobian projective coords:
+ *      x = X/Z**2, y=Y/Z**3
+ *
+ * The output points p, s, and r are converted to standard (homogeneous)
+ * projective coords:
+ *      x = X/Z, y=Y/Z
+ */
+int ec_GFp_simple_ladder_pre(const EC_GROUP *group,
+                             EC_POINT *r, EC_POINT *s,
+                             EC_POINT *p, BN_CTX *ctx)
+{
+    BIGNUM *t1, *t2, *t3, *t4, *t5, *t6 = NULL;
+
+    t1 = r->Z;
+    t2 = r->Y;
+    t3 = s->X;
+    t4 = r->X;
+    t5 = s->Y;
+    t6 = s->Z;
+
+    /* convert p: (X,Y,Z) -> (XZ,Y,Z**3) */
+    if (!group->meth->field_mul(group, p->X, p->X, p->Z, ctx)
+        || !group->meth->field_sqr(group, t1, p->Z, ctx)
+        || !group->meth->field_mul(group, p->Z, p->Z, t1, ctx)
+        /* r := 2p */
+        || !group->meth->field_sqr(group, t2, p->X, ctx)
+        || !group->meth->field_sqr(group, t3, p->Z, ctx)
+        || !group->meth->field_mul(group, t4, t3, group->a, ctx)
+        || !BN_mod_sub_quick(t5, t2, t4, group->field)
+        || !BN_mod_add_quick(t2, t2, t4, group->field)
+        || !group->meth->field_sqr(group, t5, t5, ctx)
+        || !group->meth->field_mul(group, t6, t3, group->b, ctx)
+        || !group->meth->field_mul(group, t1, p->X, p->Z, ctx)
+        || !group->meth->field_mul(group, t4, t1, t6, ctx)
+        || !BN_mod_lshift_quick(t4, t4, 3, group->field)
+        /* r->X coord output */
+        || !BN_mod_sub_quick(r->X, t5, t4, group->field)
+        || !group->meth->field_mul(group, t1, t1, t2, ctx)
+        || !group->meth->field_mul(group, t2, t3, t6, ctx)
+        || !BN_mod_add_quick(t1, t1, t2, group->field)
+        /* r->Z coord output */
+        || !BN_mod_lshift_quick(r->Z, t1, 2, group->field)
+        || !EC_POINT_copy(s, p))
+        return 0;
+
+    r->Z_is_one = 0;
+    s->Z_is_one = 0;
+    p->Z_is_one = 0;
+
+    return 1;
+}
+
+/*-
+ * Differential addition-and-doubling using  Eq. (8) and (10) from Izu-Takagi
+ * "A fast parallel elliptic curve multiplication resistant against side channel
+ * attacks", as described at
+ * https://hyperelliptic.org/EFD/g1p/auto-shortw-xz.html#ladder-ladd-2002-it-3
+ */
+int ec_GFp_simple_ladder_step(const EC_GROUP *group,
+                              EC_POINT *r, EC_POINT *s,
+                              EC_POINT *p, BN_CTX *ctx)
+{
+    int ret = 0;
+    BIGNUM *t0, *t1, *t2, *t3, *t4, *t5, *t6, *t7 = NULL;
+
+    BN_CTX_start(ctx);
+    t0 = BN_CTX_get(ctx);
+    t1 = BN_CTX_get(ctx);
+    t2 = BN_CTX_get(ctx);
+    t3 = BN_CTX_get(ctx);
+    t4 = BN_CTX_get(ctx);
+    t5 = BN_CTX_get(ctx);
+    t6 = BN_CTX_get(ctx);
+    t7 = BN_CTX_get(ctx);
+
+    if (t7 == NULL
+        || !group->meth->field_mul(group, t0, r->X, s->X, ctx)
+        || !group->meth->field_mul(group, t1, r->Z, s->Z, ctx)
+        || !group->meth->field_mul(group, t2, r->X, s->Z, ctx)
+        || !group->meth->field_mul(group, t3, r->Z, s->X, ctx)
+        || !group->meth->field_mul(group, t4, group->a, t1, ctx)
+        || !BN_mod_sub_quick(t4, t0, t4, group->field)
+        || !BN_mod_add_quick(t5, t3, t2, group->field)
+        || !group->meth->field_sqr(group, t4, t4, ctx)
+        || !group->meth->field_mul(group, t5, t1, t5, ctx)
+        || !BN_mod_lshift_quick(t0, group->b, 2, group->field)
+        || !group->meth->field_mul(group, t5, t0, t5, ctx)
+        || !BN_mod_sub_quick(t5, t4, t5, group->field)
+        /* s->X coord output */
+        || !group->meth->field_mul(group, s->X, t5, p->Z, ctx)
+        || !BN_mod_sub_quick(t3, t2, t3, group->field)
+        || !group->meth->field_sqr(group, t3, t3, ctx)
+        /* s->Z coord output */
+        || !group->meth->field_mul(group, s->Z, t3, p->X, ctx)
+        || !group->meth->field_sqr(group, t2, r->X, ctx)
+        || !group->meth->field_sqr(group, t4, r->Z, ctx)
+        || !group->meth->field_mul(group, t1, t4, group->a, ctx)
+        || !BN_mod_add_quick(t6, r->X, r->Z, group->field)
+        || !group->meth->field_sqr(group, t6, t6, ctx)
+        || !BN_mod_sub_quick(t6, t6, t2, group->field)
+        || !BN_mod_sub_quick(t6, t6, t4, group->field)
+        || !BN_mod_sub_quick(t7, t2, t1, group->field)
+        || !group->meth->field_sqr(group, t7, t7, ctx)
+        || !group->meth->field_mul(group, t5, t4, t6, ctx)
+        || !group->meth->field_mul(group, t5, t0, t5, ctx)
+        /* r->X coord output */
+        || !BN_mod_sub_quick(r->X, t7, t5, group->field)
+        || !BN_mod_add_quick(t2, t2, t1, group->field)
+        || !group->meth->field_sqr(group, t5, t4, ctx)
+        || !group->meth->field_mul(group, t5, t5, t0, ctx)
+        || !group->meth->field_mul(group, t6, t6, t2, ctx)
+        || !BN_mod_lshift1_quick(t6, t6, group->field)
+        /* r->Z coord output */
+        || !BN_mod_add_quick(r->Z, t5, t6, group->field))
+        goto err;
+
+    ret = 1;
+
+ err:
+    BN_CTX_end(ctx);
+    return ret;
+}
+
+/*-
+ * Recovers the y-coordinate of r using Eq. (8) from Brier-Joye, "Weierstrass
+ * Elliptic Curves and Side-Channel Attacks", modified to work in projective
+ * coordinates and return r in Jacobian projective coordinates.
+ *
+ * X4 = two*Y1*X2*Z3*Z2*Z1;
+ * Y4 = two*b*Z3*SQR(Z2*Z1) + Z3*(a*Z2*Z1+X1*X2)*(X1*Z2+X2*Z1) - X3*SQR(X1*Z2-X2*Z1);
+ * Z4 = two*Y1*Z3*SQR(Z2)*Z1;
+ *
+ * Z4 != 0 because:
+ *  - Z1==0 implies p is at infinity, which would have caused an early exit in
+ *    the caller;
+ *  - Z2==0 implies r is at infinity (handled by the BN_is_zero(r->Z) branch);
+ *  - Z3==0 implies s is at infinity (handled by the BN_is_zero(s->Z) branch);
+ *  - Y1==0 implies p has order 2, so either r or s are infinity and handled by
+ *    one of the BN_is_zero(...) branches.
+ */
+int ec_GFp_simple_ladder_post(const EC_GROUP *group,
+                              EC_POINT *r, EC_POINT *s,
+                              EC_POINT *p, BN_CTX *ctx)
+{
+    int ret = 0;
+    BIGNUM *t0, *t1, *t2, *t3, *t4, *t5, *t6 = NULL;
+
+    if (BN_is_zero(r->Z))
+        return EC_POINT_set_to_infinity(group, r);
+
+    if (BN_is_zero(s->Z)) {
+        /* (X,Y,Z) -> (XZ,YZ**2,Z) */
+        if (!group->meth->field_mul(group, r->X, p->X, p->Z, ctx)
+            || !group->meth->field_sqr(group, r->Z, p->Z, ctx)
+            || !group->meth->field_mul(group, r->Y, p->Y, r->Z, ctx)
+            || !BN_copy(r->Z, p->Z)
+            || !EC_POINT_invert(group, r, ctx))
+            return 0;
+        return 1;
+    }
+
+    BN_CTX_start(ctx);
+    t0 = BN_CTX_get(ctx);
+    t1 = BN_CTX_get(ctx);
+    t2 = BN_CTX_get(ctx);
+    t3 = BN_CTX_get(ctx);
+    t4 = BN_CTX_get(ctx);
+    t5 = BN_CTX_get(ctx);
+    t6 = BN_CTX_get(ctx);
+
+    if (t6 == NULL
+        || !BN_mod_lshift1_quick(t0, p->Y, group->field)
+        || !group->meth->field_mul(group, t1, r->X, p->Z, ctx)
+        || !group->meth->field_mul(group, t2, r->Z, s->Z, ctx)
+        || !group->meth->field_mul(group, t2, t1, t2, ctx)
+        || !group->meth->field_mul(group, t3, t2, t0, ctx)
+        || !group->meth->field_mul(group, t2, r->Z, p->Z, ctx)
+        || !group->meth->field_sqr(group, t4, t2, ctx)
+        || !BN_mod_lshift1_quick(t5, group->b, group->field)
+        || !group->meth->field_mul(group, t4, t4, t5, ctx)
+        || !group->meth->field_mul(group, t6, t2, group->a, ctx)
+        || !group->meth->field_mul(group, t5, r->X, p->X, ctx)
+        || !BN_mod_add_quick(t5, t6, t5, group->field)
+        || !group->meth->field_mul(group, t6, r->Z, p->X, ctx)
+        || !BN_mod_add_quick(t2, t6, t1, group->field)
+        || !group->meth->field_mul(group, t5, t5, t2, ctx)
+        || !BN_mod_sub_quick(t6, t6, t1, group->field)
+        || !group->meth->field_sqr(group, t6, t6, ctx)
+        || !group->meth->field_mul(group, t6, t6, s->X, ctx)
+        || !BN_mod_add_quick(t4, t5, t4, group->field)
+        || !group->meth->field_mul(group, t4, t4, s->Z, ctx)
+        || !BN_mod_sub_quick(t4, t4, t6, group->field)
+        || !group->meth->field_sqr(group, t5, r->Z, ctx)
+        || !group->meth->field_mul(group, r->Z, p->Z, s->Z, ctx)
+        || !group->meth->field_mul(group, r->Z, t5, r->Z, ctx)
+        || !group->meth->field_mul(group, r->Z, r->Z, t0, ctx)
+        /* t3 := X, t4 := Y */
+        /* (X,Y,Z) -> (XZ,YZ**2,Z) */
+        || !group->meth->field_mul(group, r->X, t3, r->Z, ctx)
+        || !group->meth->field_sqr(group, t3, r->Z, ctx)
+        || !group->meth->field_mul(group, r->Y, t4, t3, ctx))
+        goto err;
+
+    ret = 1;
+
+ err:
+    BN_CTX_end(ctx);
+    return ret;
 }