Provide an EC method that uses homogeneous projective coordinates.

This makes use of EC_FIELD_ELEMENT to perform fixed width constant
time operations.

Addition and doubling of points makes use of the formulas from
"Complete addition formulas for prime order elliptic curves"
(https://eprint.iacr.org/2015/1060). These are complete and
operate in constant time.

Further work will continue in tree.

ok tb@

3 files changed, 870 insertions, 2 deletions

diff --git a/src/lib/libcrypto/Makefile b/src/lib/libcrypto/Makefile
index 7564001961..6adfe45e6f 100644
--- a/src/lib/libcrypto/Makefile
+++ b/src/lib/libcrypto/Makefile
@@ -1,4 +1,4 @@
-# $OpenBSD: Makefile,v 1.234 2025/05/25 05:12:05 jsing Exp $
+# $OpenBSD: Makefile,v 1.235 2025/05/25 05:19:26 jsing Exp $
 
 LIB=    crypto
 LIBREBUILD=y
@@ -289,6 +289,7 @@ SRCS+= ec_lib.c
 SRCS+= ec_mult.c
 SRCS+= ec_pmeth.c
 SRCS+= eck_prn.c
+SRCS+= ecp_hp_methods.c
 SRCS+= ecp_methods.c
 SRCS+= ecx_methods.c
 
diff --git a/src/lib/libcrypto/ec/ec_local.h b/src/lib/libcrypto/ec/ec_local.h
index c0ff026fb2..75a3e25247 100644
--- a/src/lib/libcrypto/ec/ec_local.h
+++ b/src/lib/libcrypto/ec/ec_local.h
@@ -1,4 +1,4 @@
-/* $OpenBSD: ec_local.h,v 1.68 2025/05/24 08:25:58 jsing Exp $ */
+/* $OpenBSD: ec_local.h,v 1.69 2025/05/25 05:19:26 jsing Exp $ */
 /*
  * Originally written by Bodo Moeller for the OpenSSL project.
  */
@@ -76,6 +76,7 @@
 #include <openssl/objects.h>
 
 #include "bn_local.h"
+#include "ec_internal.h"
 
 __BEGIN_HIDDEN_DECLS
 
@@ -158,6 +159,10 @@ struct ec_group_st {
 
         /* Montgomery context used by EC_GFp_mont_method. */
         BN_MONT_CTX *mont_ctx;
+
+        EC_FIELD_MODULUS fm;
+        EC_FIELD_ELEMENT fe_a;
+        EC_FIELD_ELEMENT fe_b;
 } /* EC_GROUP */;
 
 struct ec_point_st {
@@ -171,6 +176,10 @@ struct ec_point_st {
         BIGNUM *Y;
         BIGNUM *Z;
         int Z_is_one; /* enable optimized point arithmetics for special case */
+
+        EC_FIELD_ELEMENT fe_x;
+        EC_FIELD_ELEMENT fe_y;
+        EC_FIELD_ELEMENT fe_z;
 } /* EC_POINT */;
 
 const EC_METHOD *EC_GFp_simple_method(void);
diff --git a/src/lib/libcrypto/ec/ecp_hp_methods.c b/src/lib/libcrypto/ec/ecp_hp_methods.c
new file mode 100644
index 0000000000..0c6b095765
--- /dev/null
+++ b/src/lib/libcrypto/ec/ecp_hp_methods.c
@@ -0,0 +1,858 @@
+/*      $OpenBSD: ecp_hp_methods.c,v 1.1 2025/05/25 05:19:26 jsing Exp $        */
+/*
+ * Copyright (c) 2024-2025 Joel Sing <jsing@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 <string.h>
+
+#include <openssl/bn.h>
+#include <openssl/ec.h>
+#include <openssl/err.h>
+
+#include "bn_internal.h"
+#include "ec_local.h"
+#include "ec_internal.h"
+#include "err_local.h"
+
+/*
+ * TODO:
+ *  - Constant time and blinding for scalar multiplication
+ */
+
+static int
+ec_group_set_curve(EC_GROUP *group, const BIGNUM *p, const BIGNUM *a,
+    const BIGNUM *b, BN_CTX *ctx)
+{
+        BIGNUM *t;
+        int ret = 0;
+
+        BN_CTX_start(ctx);
+
+        /* XXX - p must be a prime > 3. */
+
+        if (!bn_copy(group->p, p))
+                goto err;
+        if (!bn_copy(group->a, a))
+                goto err;
+        if (!bn_copy(group->b, b))
+                goto err;
+
+        /* XXX */
+        BN_set_negative(group->p, 0);
+
+        /* XXX */
+        if (!BN_nnmod(group->a, group->a, group->p, ctx))
+                goto err;
+        if (!BN_nnmod(group->b, group->b, group->p, ctx))
+                goto err;
+
+        if ((t = BN_CTX_get(ctx)) == NULL)
+                goto err;
+        if (!BN_set_word(t, 3))
+                goto err;
+        if (!BN_mod_add(t, t, a, group->p, ctx))
+                goto err;
+
+        group->a_is_minus3 = BN_is_zero(t);
+
+        if (!ec_field_modulus_from_bn(&group->fm, group->p, ctx))
+                goto err;
+        if (!ec_field_element_from_bn(&group->fm, group, &group->fe_a, group->a, ctx))
+                goto err;
+        if (!ec_field_element_from_bn(&group->fm, group, &group->fe_b, group->b, ctx))
+                goto err;
+
+        ret = 1;
+
+ err:
+        BN_CTX_end(ctx);
+
+        return ret;
+}
+
+static int
+ec_group_get_curve(const EC_GROUP *group, BIGNUM *p, BIGNUM *a,
+    BIGNUM *b, BN_CTX *ctx)
+{
+        if (p != NULL) {
+                if (!bn_copy(p, group->p))
+                        return 0;
+        }
+        if (a != NULL) {
+                if (!bn_copy(a, group->a))
+                        return 0;
+        }
+        if (b != NULL) {
+                if (!bn_copy(b, group->b))
+                        return 0;
+        }
+        return 1;
+}
+
+static int
+ec_point_is_at_infinity(const EC_GROUP *group, const EC_POINT *point)
+{
+        /* Check if Z is equal to zero. */
+        return ec_field_element_is_zero(&group->fm, &point->fe_z);
+}
+
+static int
+ec_point_set_to_infinity(const EC_GROUP *group, EC_POINT *point)
+{
+        /* Infinity is (x = 0, y = 1, z = 0). */
+
+        memset(&point->fe_x, 0, sizeof(point->fe_x));
+        memset(&point->fe_y, 0, sizeof(point->fe_y));
+        memset(&point->fe_z, 0, sizeof(point->fe_z));
+
+        point->fe_y.w[0] = 1;
+
+        return 1;
+}
+
+static int
+ec_point_set_affine_coordinates(const EC_GROUP *group, EC_POINT *point,
+    const BIGNUM *x, const BIGNUM *y, BN_CTX *ctx)
+{
+        if (x == NULL || y == NULL) {
+                ECerror(ERR_R_PASSED_NULL_PARAMETER);
+                return 0;
+        }
+
+        if (!bn_copy(point->X, x))
+                return 0;
+        if (!bn_copy(point->Y, y))
+                return 0;
+        if (!BN_one(point->Z))
+                return 0;
+
+        /* XXX */
+        if (!BN_nnmod(point->X, point->X, group->p, ctx))
+                return 0;
+        if (!BN_nnmod(point->Y, point->Y, group->p, ctx))
+                return 0;
+
+        if (!ec_field_element_from_bn(&group->fm, group, &point->fe_x, point->X, ctx))
+                return 0;
+        if (!ec_field_element_from_bn(&group->fm, group, &point->fe_y, point->Y, ctx))
+                return 0;
+        if (!ec_field_element_from_bn(&group->fm, group, &point->fe_z, point->Z, ctx))
+                return 0;
+
+        return 1;
+}
+
+static int
+ec_point_get_affine_coordinates(const EC_GROUP *group, const EC_POINT *point,
+    BIGNUM *x, BIGNUM *y, BN_CTX *ctx)
+{
+        BIGNUM *zinv;
+        int ret = 0;
+
+        /*
+         * Convert homogeneous projective coordinates (XZ, YZ, Z) to affine
+         * coordinates (x = X/Z, y = Y/Z).
+         */
+        if (!ec_field_element_to_bn(&group->fm, &point->fe_x, point->X, ctx))
+                return 0;
+        if (!ec_field_element_to_bn(&group->fm, &point->fe_y, point->Y, ctx))
+                return 0;
+        if (!ec_field_element_to_bn(&group->fm, &point->fe_z, point->Z, ctx))
+                return 0;
+
+        BN_CTX_start(ctx);
+
+        if ((zinv = BN_CTX_get(ctx)) == NULL)
+                goto err;
+
+        if (BN_mod_inverse_ct(zinv, point->Z, group->p, ctx) == NULL)
+                goto err;
+
+        if (x != NULL) {
+                if (!BN_mod_mul(x, point->X, zinv, group->p, ctx))
+                        goto err;
+        }
+        if (y != NULL) {
+                if (!BN_mod_mul(y, point->Y, zinv, group->p, ctx))
+                        goto err;
+        }
+
+        ret = 1;
+
+ err:
+        BN_CTX_end(ctx);
+
+        return ret;
+}
+
+static int
+ec_point_add_a1(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a,
+    const EC_POINT *b, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3;
+        EC_FIELD_ELEMENT b3, t0, t1, t2, t3, t4, t5;
+        EC_FIELD_ELEMENT ga, gb;
+
+        /*
+         * Complete, projective point addition for arbitrary prime order short
+         * Weierstrass curves with arbitrary a - see
+         * https://eprint.iacr.org/2015/1060, algorithm 1 and appendix A.1.
+         */
+
+        ec_field_element_copy(&ga, &group->fe_a);
+        ec_field_element_copy(&gb, &group->fe_b);
+
+        ec_field_element_copy(&X1, &a->fe_x);
+        ec_field_element_copy(&Y1, &a->fe_y);
+        ec_field_element_copy(&Z1, &a->fe_z);
+
+        ec_field_element_copy(&X2, &b->fe_x);
+        ec_field_element_copy(&Y2, &b->fe_y);
+        ec_field_element_copy(&Z2, &b->fe_z);
+
+        /* b3 := 3 * b ; */
+        ec_field_element_add(&group->fm, &b3, &gb, &gb);
+        ec_field_element_add(&group->fm, &b3, &b3, &gb);
+
+        /* t0 := X1 * X2 ; t1 := Y1 * Y2 ; t2 := Z1 * Z2 ; */
+        ec_field_element_mul(&group->fm, &t0, &X1, &X2);
+        ec_field_element_mul(&group->fm, &t1, &Y1, &Y2);
+        ec_field_element_mul(&group->fm, &t2, &Z1, &Z2);
+
+        /* t3 := X1 + Y1 ; t4 := X2 + Y2 ; t3 := t3 * t4 ; */
+        ec_field_element_add(&group->fm, &t3, &X1, &Y1);
+        ec_field_element_add(&group->fm, &t4, &X2, &Y2);
+        ec_field_element_mul(&group->fm, &t3, &t3, &t4);
+
+        /* t4 := t0 + t1 ; t3 := t3 - t4 ; t4 := X1 + Z1 ; */
+        ec_field_element_add(&group->fm, &t4, &t0, &t1);
+        ec_field_element_sub(&group->fm, &t3, &t3, &t4);
+        ec_field_element_add(&group->fm, &t4, &X1, &Z1);
+
+        /* t5 := X2 + Z2 ; t4 := t4 * t5 ; t5 := t0 + t2 ; */
+        ec_field_element_add(&group->fm, &t5, &X2, &Z2);
+        ec_field_element_mul(&group->fm, &t4, &t4, &t5);
+        ec_field_element_add(&group->fm, &t5, &t0, &t2);
+
+        /* t4 := t4 - t5 ; t5 := Y1 + Z1 ; X3 := Y2 + Z2 ; */
+        ec_field_element_sub(&group->fm, &t4, &t4, &t5);
+        ec_field_element_add(&group->fm, &t5, &Y1, &Z1);
+        ec_field_element_add(&group->fm, &X3, &Y2, &Z2);
+
+        /* t5 := t5 * X3 ; X3 := t1 + t2 ; t5 := t5 - X3 ; */
+        ec_field_element_mul(&group->fm, &t5, &t5, &X3);
+        ec_field_element_add(&group->fm, &X3, &t1, &t2);
+        ec_field_element_sub(&group->fm, &t5, &t5, &X3);
+
+        /* Z3 := a * t4 ; X3 := b3 * t2 ; Z3 := X3 + Z3 ; */
+        ec_field_element_mul(&group->fm, &Z3, &ga, &t4);
+        ec_field_element_mul(&group->fm, &X3, &b3, &t2);
+        ec_field_element_add(&group->fm, &Z3, &X3, &Z3);
+
+        /* X3 := t1 - Z3 ; Z3 := t1 + Z3 ; Y3 := X3 * Z3 ; */
+        ec_field_element_sub(&group->fm, &X3, &t1, &Z3);
+        ec_field_element_add(&group->fm, &Z3, &t1, &Z3);
+        ec_field_element_mul(&group->fm, &Y3, &X3, &Z3);
+
+        /* t1 := t0 + t0 ; t1 := t1 + t0 ; t2 := a * t2 ; */
+        ec_field_element_add(&group->fm, &t1, &t0, &t0);
+        ec_field_element_add(&group->fm, &t1, &t1, &t0);
+        ec_field_element_mul(&group->fm, &t2, &ga, &t2);
+
+        /* t4 := b3 * t4 ; t1 := t1 + t2 ; t2 := t0 - t2 ; */
+        ec_field_element_mul(&group->fm, &t4, &b3, &t4);
+        ec_field_element_add(&group->fm, &t1, &t1, &t2);
+        ec_field_element_sub(&group->fm, &t2, &t0, &t2);
+
+        /* t2 := a * t2 ; t4 := t4 + t2 ; t0 := t1 * t4 ; */
+        ec_field_element_mul(&group->fm, &t2, &ga, &t2);
+        ec_field_element_add(&group->fm, &t4, &t4, &t2);
+        ec_field_element_mul(&group->fm, &t0, &t1, &t4);
+
+        /* Y3 := Y3 + t0 ; t0 := t5 * t4 ; X3 := t3 * X3 ; */
+        ec_field_element_add(&group->fm, &Y3, &Y3, &t0);
+        ec_field_element_mul(&group->fm, &t0, &t5, &t4);
+        ec_field_element_mul(&group->fm, &X3, &t3, &X3);
+
+        /* X3 := X3 - t0 ; t0 := t3 * t1 ; Z3 := t5 * Z3 ; */
+        ec_field_element_sub(&group->fm, &X3, &X3, &t0);
+        ec_field_element_mul(&group->fm, &t0, &t3, &t1);
+        ec_field_element_mul(&group->fm, &Z3, &t5, &Z3);
+
+        /* Z3 := Z3 + t0 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &t0);
+
+        ec_field_element_copy(&r->fe_x, &X3);
+        ec_field_element_copy(&r->fe_y, &Y3);
+        ec_field_element_copy(&r->fe_z, &Z3);
+
+        return 1;
+}
+
+static int
+ec_point_add_a2(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a,
+    const EC_POINT *b, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3;
+        EC_FIELD_ELEMENT t0, t1, t2, t3, t4;
+        EC_FIELD_ELEMENT gb;
+
+        /*
+         * Complete, projective point addition for arbitrary prime order short
+         * Weierstrass curves with a = -3 - see https://eprint.iacr.org/2015/1060,
+         * algorithm 4 and appendix A.2.
+         */
+
+        ec_field_element_copy(&gb, &group->fe_b);
+
+        ec_field_element_copy(&X1, &a->fe_x);
+        ec_field_element_copy(&Y1, &a->fe_y);
+        ec_field_element_copy(&Z1, &a->fe_z);
+
+        ec_field_element_copy(&X2, &b->fe_x);
+        ec_field_element_copy(&Y2, &b->fe_y);
+        ec_field_element_copy(&Z2, &b->fe_z);
+
+        /* t0 := X1 * X2 ; t1 := Y1 * Y2 ; t2 := Z1 * Z2 ; */
+        ec_field_element_mul(&group->fm, &t0, &X1, &X2);
+        ec_field_element_mul(&group->fm, &t1, &Y1, &Y2);
+        ec_field_element_mul(&group->fm, &t2, &Z1, &Z2);
+
+        /* t3 := X1 + Y1 ; t4 := X2 + Y2 ; t3 := t3 * t4 ; */
+        ec_field_element_add(&group->fm, &t3, &X1, &Y1);
+        ec_field_element_add(&group->fm, &t4, &X2, &Y2);
+        ec_field_element_mul(&group->fm, &t3, &t3, &t4);
+
+        /* t4 := t0 + t1 ; t3 := t3 - t4 ; t4 := Y1 + Z1 ; */
+        ec_field_element_add(&group->fm, &t4, &t0, &t1);
+        ec_field_element_sub(&group->fm, &t3, &t3, &t4);
+        ec_field_element_add(&group->fm, &t4, &Y1, &Z1);
+
+        /* X3 := Y2 + Z2 ; t4 := t4 * X3 ; X3 := t1 + t2 ; */
+        ec_field_element_add(&group->fm, &X3, &Y2, &Z2);
+        ec_field_element_mul(&group->fm, &t4, &t4, &X3);
+        ec_field_element_add(&group->fm, &X3, &t1, &t2);
+
+        /* t4 := t4 - X3 ; X3 := X1 + Z1 ; Y3 := X2 + Z2 ; */
+        ec_field_element_sub(&group->fm, &t4, &t4, &X3);
+        ec_field_element_add(&group->fm, &X3, &X1, &Z1);
+        ec_field_element_add(&group->fm, &Y3, &X2, &Z2);
+
+        /* X3 := X3 * Y3 ; Y3 := t0 + t2 ; Y3 := X3 - Y3 ; */
+        ec_field_element_mul(&group->fm, &X3, &X3, &Y3);
+        ec_field_element_add(&group->fm, &Y3, &t0, &t2);
+        ec_field_element_sub(&group->fm, &Y3, &X3, &Y3);
+
+        /* Z3 := b * t2 ; X3 := Y3 - Z3 ; Z3 := X3 + X3 ; */
+        ec_field_element_mul(&group->fm, &Z3, &gb, &t2);
+        ec_field_element_sub(&group->fm, &X3, &Y3, &Z3);
+        ec_field_element_add(&group->fm, &Z3, &X3, &X3);
+
+        /* X3 := X3 + Z3 ; Z3 := t1 - X3 ; X3 := t1 + X3 ; */
+        ec_field_element_add(&group->fm, &X3, &X3, &Z3);
+        ec_field_element_sub(&group->fm, &Z3, &t1, &X3);
+        ec_field_element_add(&group->fm, &X3, &t1, &X3);
+
+        /* Y3 := b * Y3 ; t1 := t2 + t2 ; t2 := t1 + t2 ; */
+        ec_field_element_mul(&group->fm, &Y3, &gb, &Y3);
+        ec_field_element_add(&group->fm, &t1, &t2, &t2);
+        ec_field_element_add(&group->fm, &t2, &t1, &t2);
+
+        /* Y3 := Y3 - t2 ; Y3 := Y3 - t0 ; t1 := Y3 + Y3 ; */
+        ec_field_element_sub(&group->fm, &Y3, &Y3, &t2);
+        ec_field_element_sub(&group->fm, &Y3, &Y3, &t0);
+        ec_field_element_add(&group->fm, &t1, &Y3, &Y3);
+
+        /* Y3 := t1 + Y3 ; t1 := t0 + t0 ; t0 := t1 + t0 ; */
+        ec_field_element_add(&group->fm, &Y3, &t1, &Y3);
+        ec_field_element_add(&group->fm, &t1, &t0, &t0);
+        ec_field_element_add(&group->fm, &t0, &t1, &t0);
+
+        /* t0 := t0 - t2 ; t1 := t4 * Y3 ; t2 := t0 * Y3 ; */
+        ec_field_element_sub(&group->fm, &t0, &t0, &t2);
+        ec_field_element_mul(&group->fm, &t1, &t4, &Y3);
+        ec_field_element_mul(&group->fm, &t2, &t0, &Y3);
+
+        /* Y3 := X3 * Z3 ; Y3 := Y3 + t2 ; X3 := t3 * X3 ; */
+        ec_field_element_mul(&group->fm, &Y3, &X3, &Z3);
+        ec_field_element_add(&group->fm, &Y3, &Y3, &t2);
+        ec_field_element_mul(&group->fm, &X3, &t3, &X3);
+
+        /* X3 := X3 - t1 ; Z3 := t4 * Z3 ; t1 := t3 * t0 ; */
+        ec_field_element_sub(&group->fm, &X3, &X3, &t1);
+        ec_field_element_mul(&group->fm, &Z3, &t4, &Z3);
+        ec_field_element_mul(&group->fm, &t1, &t3, &t0);
+
+        /* Z3 := Z3 + t1 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &t1);
+
+        ec_field_element_copy(&r->fe_x, &X3);
+        ec_field_element_copy(&r->fe_y, &Y3);
+        ec_field_element_copy(&r->fe_z, &Z3);
+
+        return 1;
+}
+
+static int
+ec_point_add(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a,
+    const EC_POINT *b, BN_CTX *ctx)
+{
+        if (group->a_is_minus3)
+                return ec_point_add_a2(group, r, a, b, ctx);
+
+        return ec_point_add_a1(group, r, a, b, ctx);
+}
+
+static int
+ec_point_dbl_a1(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT X1, Y1, Z1, X3, Y3, Z3;
+        EC_FIELD_ELEMENT b3, t0, t1, t2, t3;
+        EC_FIELD_ELEMENT ga, gb;
+
+        /*
+         * Exception-free point doubling for arbitrary prime order short
+         * Weierstrass curves with arbitrary a - see
+         * https://eprint.iacr.org/2015/1060, algorithm 3 and appendix A.1.
+         */
+
+        ec_field_element_copy(&ga, &group->fe_a);
+        ec_field_element_copy(&gb, &group->fe_b);
+
+        ec_field_element_copy(&X1, &a->fe_x);
+        ec_field_element_copy(&Y1, &a->fe_y);
+        ec_field_element_copy(&Z1, &a->fe_z);
+
+        /* b3 := 3 * b ; */
+        ec_field_element_add(&group->fm, &b3, &gb, &gb);
+        ec_field_element_add(&group->fm, &b3, &b3, &gb);
+
+        /* b3 := 3 * b ; */
+        ec_field_element_add(&group->fm, &b3, &gb, &gb);
+        ec_field_element_add(&group->fm, &b3, &b3, &gb);
+
+        /* t0 := X^2; t1 := Y^2; t2 := Z^2 ; */
+        ec_field_element_sqr(&group->fm, &t0, &X1);
+        ec_field_element_sqr(&group->fm, &t1, &Y1);
+        ec_field_element_sqr(&group->fm, &t2, &Z1);
+
+        /* t3 := X * Y ; t3 := t3 + t3 ; Z3 := X * Z ; */
+        ec_field_element_mul(&group->fm, &t3, &X1, &Y1);
+        ec_field_element_add(&group->fm, &t3, &t3, &t3);
+        ec_field_element_mul(&group->fm, &Z3, &X1, &Z1);
+
+        /* Z3 := Z3 + Z3 ; X3 := a * Z3 ; Y3 := b3 * t2 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &Z3);
+        ec_field_element_mul(&group->fm, &X3, &ga, &Z3);
+        ec_field_element_mul(&group->fm, &Y3, &b3, &t2);
+
+        /* Y3 := X3 + Y3 ; X3 := t1 - Y3 ; Y3 := t1 + Y3 ; */
+        ec_field_element_add(&group->fm, &Y3, &X3, &Y3);
+        ec_field_element_sub(&group->fm, &X3, &t1, &Y3);
+        ec_field_element_add(&group->fm, &Y3, &t1, &Y3);
+
+        /* Y3 := X3 * Y3 ; X3 := t3 * X3 ; Z3 := b3 * Z3 ; */
+        ec_field_element_mul(&group->fm, &Y3, &X3, &Y3);
+        ec_field_element_mul(&group->fm, &X3, &t3, &X3);
+        ec_field_element_mul(&group->fm, &Z3, &b3, &Z3);
+
+        /* t2 := a * t2 ; t3 := t0 - t2 ; t3 := a * t3 ; */
+        ec_field_element_mul(&group->fm, &t2, &ga, &t2);
+        ec_field_element_sub(&group->fm, &t3, &t0, &t2);
+        ec_field_element_mul(&group->fm, &t3, &ga, &t3);
+
+        /* t3 := t3 + Z3 ; Z3 := t0 + t0 ; t0 := Z3 + t0 ; */
+        ec_field_element_add(&group->fm, &t3, &t3, &Z3);
+        ec_field_element_add(&group->fm, &Z3, &t0, &t0);
+        ec_field_element_add(&group->fm, &t0, &Z3, &t0);
+
+        /* t0 := t0 + t2 ; t0 := t0 * t3 ; Y3 := Y3 + t0 ; */
+        ec_field_element_add(&group->fm, &t0, &t0, &t2);
+        ec_field_element_mul(&group->fm, &t0, &t0, &t3);
+        ec_field_element_add(&group->fm, &Y3, &Y3, &t0);
+
+        /* t2 := Y * Z ; t2 := t2 + t2 ; t0 := t2 * t3 ; */
+        ec_field_element_mul(&group->fm, &t2, &Y1, &Z1);
+        ec_field_element_add(&group->fm, &t2, &t2, &t2);
+        ec_field_element_mul(&group->fm, &t0, &t2, &t3);
+
+        /* X3 := X3 - t0 ; Z3 := t2 * t1 ; Z3 := Z3 + Z3 ; */
+        ec_field_element_sub(&group->fm, &X3, &X3, &t0);
+        ec_field_element_mul(&group->fm, &Z3, &t2, &t1);
+        ec_field_element_add(&group->fm, &Z3, &Z3, &Z3);
+
+        /* Z3 := Z3 + Z3 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &Z3);
+
+        ec_field_element_copy(&r->fe_x, &X3);
+        ec_field_element_copy(&r->fe_y, &Y3);
+        ec_field_element_copy(&r->fe_z, &Z3);
+
+        return 1;
+}
+
+static int
+ec_point_dbl_a2(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT X1, Y1, Z1, X3, Y3, Z3;
+        EC_FIELD_ELEMENT t0, t1, t2, t3;
+        EC_FIELD_ELEMENT ga, gb;
+
+        /*
+         * Exception-free point doubling for arbitrary prime order short
+         * Weierstrass curves with a = -3 - see https://eprint.iacr.org/2015/1060,
+         * algorithm 6 and appendix A.2.
+         */
+
+        ec_field_element_copy(&ga, &group->fe_a);
+        ec_field_element_copy(&gb, &group->fe_b);
+
+        ec_field_element_copy(&X1, &a->fe_x);
+        ec_field_element_copy(&Y1, &a->fe_y);
+        ec_field_element_copy(&Z1, &a->fe_z);
+
+        /* t0 := X^2; t1 := Y^2; t2 := Z^2 ; */
+        ec_field_element_sqr(&group->fm, &t0, &X1);
+        ec_field_element_sqr(&group->fm, &t1, &Y1);
+        ec_field_element_sqr(&group->fm, &t2, &Z1);
+
+        /* t3 := X * Y ; t3 := t3 + t3 ; Z3 := X * Z ; */
+        ec_field_element_mul(&group->fm, &t3, &X1, &Y1);
+        ec_field_element_add(&group->fm, &t3, &t3, &t3);
+        ec_field_element_mul(&group->fm, &Z3, &X1, &Z1);
+
+        /* Z3 := Z3 + Z3 ; Y3 := b * t2 ; Y3 := Y3 - Z3 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &Z3);
+        ec_field_element_mul(&group->fm, &Y3, &gb, &t2);
+        ec_field_element_sub(&group->fm, &Y3, &Y3, &Z3);
+
+        /* X3 := Y3 + Y3 ; Y3 := X3 + Y3 ; X3 := t1 - Y3 ; */
+        ec_field_element_add(&group->fm, &X3, &Y3, &Y3);
+        ec_field_element_add(&group->fm, &Y3, &X3, &Y3);
+        ec_field_element_sub(&group->fm, &X3, &t1, &Y3);
+
+        /* Y3 := t1 + Y3 ; Y3 := X3 * Y3 ; X3 := X3 * t3 ; */
+        ec_field_element_add(&group->fm, &Y3, &t1, &Y3);
+        ec_field_element_mul(&group->fm, &Y3, &X3, &Y3);
+        ec_field_element_mul(&group->fm, &X3, &X3, &t3);
+
+        /* t3 := t2 + t2 ; t2 := t2 + t3 ; Z3 := b * Z3 ; */
+        ec_field_element_add(&group->fm, &t3, &t2, &t2);
+        ec_field_element_add(&group->fm, &t2, &t2, &t3);
+        ec_field_element_mul(&group->fm, &Z3, &gb, &Z3);
+
+        /* Z3 := Z3 - t2 ; Z3 := Z3 - t0 ; t3 := Z3 + Z3 ; */
+        ec_field_element_sub(&group->fm, &Z3, &Z3, &t2);
+        ec_field_element_sub(&group->fm, &Z3, &Z3, &t0);
+        ec_field_element_add(&group->fm, &t3, &Z3, &Z3);
+
+        /* Z3 := Z3 + t3 ; t3 := t0 + t0 ; t0 := t3 + t0 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &t3);
+        ec_field_element_add(&group->fm, &t3, &t0, &t0);
+        ec_field_element_add(&group->fm, &t0, &t3, &t0);
+
+        /* t0 := t0 - t2 ; t0 := t0 * Z3 ; Y3 := Y3 + t0 ; */
+        ec_field_element_sub(&group->fm, &t0, &t0, &t2);
+        ec_field_element_mul(&group->fm, &t0, &t0, &Z3);
+        ec_field_element_add(&group->fm, &Y3, &Y3, &t0);
+
+        /* t0 := Y * Z ; t0 := t0 + t0 ; Z3 := t0 * Z3 ; */
+        ec_field_element_mul(&group->fm, &t0, &Y1, &Z1);
+        ec_field_element_add(&group->fm, &t0, &t0, &t0);
+        ec_field_element_mul(&group->fm, &Z3, &t0, &Z3);
+
+        /* X3 := X3 - Z3 ; Z3 := t0 * t1 ; Z3 := Z3 + Z3 ; */
+        ec_field_element_sub(&group->fm, &X3, &X3, &Z3);
+        ec_field_element_mul(&group->fm, &Z3, &t0, &t1);
+        ec_field_element_add(&group->fm, &Z3, &Z3, &Z3);
+
+        /* Z3 := Z3 + Z3 ; */
+        ec_field_element_add(&group->fm, &Z3, &Z3, &Z3);
+
+        ec_field_element_copy(&r->fe_x, &X3);
+        ec_field_element_copy(&r->fe_y, &Y3);
+        ec_field_element_copy(&r->fe_z, &Z3);
+
+        return 1;
+}
+
+static int
+ec_point_dbl(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a, BN_CTX *ctx)
+{
+        if (group->a_is_minus3)
+                return ec_point_dbl_a2(group, r, a, ctx);
+
+        return ec_point_dbl_a1(group, r, a, ctx);
+}
+
+static int
+ec_point_invert(const EC_GROUP *group, EC_POINT *point, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT y;
+        BN_ULONG mask;
+        int i;
+
+        /*
+         * Invert the point by setting Y = p - Y, if Y is non-zero and the point
+         * is not at infinity.
+         */
+
+        mask = ~(0 - (ec_point_is_at_infinity(group, point) |
+            ec_field_element_is_zero(&group->fm, &point->fe_y)));
+
+        /* XXX - masked/conditional subtraction? */
+        ec_field_element_sub(&group->fm, &y, &group->fm.m, &point->fe_y);
+
+        for (i = 0; i < group->fm.n; i++)
+                point->fe_y.w[i] = (point->fe_y.w[i] & ~mask) | (y.w[i] & mask);
+
+        return 1;
+}
+
+static int
+ec_point_is_on_curve(const EC_GROUP *group, const EC_POINT *point, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT sum, axz2, bz3, x3, y2z, z2;
+
+        /*
+         * Curve is defined by a Weierstrass equation y^2 = x^3 + a*x + b.
+         * The given point is in homogeneous projective coordinates
+         * (x = X/Z, y = Y/Z). Substitute and multiply by Z^3 in order to
+         * evaluate as zy^2 = x^3 + axz^2 + bz^3.
+         */
+
+        ec_field_element_sqr(&group->fm, &z2, &point->fe_z);
+
+        ec_field_element_sqr(&group->fm, &y2z, &point->fe_y);
+        ec_field_element_mul(&group->fm, &y2z, &y2z, &point->fe_z);
+
+        ec_field_element_sqr(&group->fm, &x3, &point->fe_x);
+        ec_field_element_mul(&group->fm, &x3, &x3, &point->fe_x);
+
+        ec_field_element_mul(&group->fm, &axz2, &group->fe_a, &point->fe_x);
+        ec_field_element_mul(&group->fm, &axz2, &axz2, &z2);
+
+        ec_field_element_mul(&group->fm, &bz3, &group->fe_b, &point->fe_z);
+        ec_field_element_mul(&group->fm, &bz3, &bz3, &z2);
+
+        ec_field_element_add(&group->fm, &sum, &x3, &axz2);
+        ec_field_element_add(&group->fm, &sum, &sum, &bz3);
+
+        return ec_field_element_equal(&group->fm, &y2z, &sum) |
+            ec_point_is_at_infinity(group, point);
+}
+
+static int
+ec_point_cmp(const EC_GROUP *group, const EC_POINT *a, const EC_POINT *b, BN_CTX *ctx)
+{
+        EC_FIELD_ELEMENT ax, ay, bx, by;
+
+        /*
+         * Compare two points that have homogeneous projection coordinates, that
+         * is (X_a/Z_a, Y_a/Z_a) == (X_b/Z_b, Y_b/Z_b). Return -1 on error, 0 on
+         * equality and 1 on inequality.
+         *
+         * If a and b are both at infinity, Z_a and Z_b will both be zero,
+         * resulting in all values becoming zero, resulting in equality. If a is
+         * at infinity and b is not, then Y_a will be one and Z_b will be
+         * non-zero, hence Y_a * Z_b will be non-zero. Z_a will be zero, hence
+         * Y_b * Z_a will be zero, resulting in inequality. The same applies if
+         * b is at infinity and a is not.
+         */
+
+        ec_field_element_mul(&group->fm, &ax, &a->fe_x, &b->fe_z);
+        ec_field_element_mul(&group->fm, &ay, &a->fe_y, &b->fe_z);
+        ec_field_element_mul(&group->fm, &bx, &b->fe_x, &a->fe_z);
+        ec_field_element_mul(&group->fm, &by, &b->fe_y, &a->fe_z);
+
+        return 1 - (ec_field_element_equal(&group->fm, &ax, &bx) &
+            ec_field_element_equal(&group->fm, &ay, &by));
+}
+
+#if 0
+static int
+ec_points_make_affine(const EC_GROUP *group, size_t num, EC_POINT *points[],
+    BN_CTX *ctx)
+{
+        size_t i;
+
+        /* XXX */
+        for (i = 0; i < num; i++) {
+                if (!EC_POINT_make_affine(group, points[0], ctx))
+                        return 0;
+        }
+
+        return 1;
+}
+#else
+
+static int
+ec_points_make_affine(const EC_GROUP *group, size_t num, EC_POINT *points[],
+    BN_CTX *ctx)
+{
+        BIGNUM **prod_Z = NULL;
+        BIGNUM *tmp, *tmp_Z;
+        size_t i;
+        int ret = 0;
+
+        if (num == 0)
+                return 1;
+
+        BN_CTX_start(ctx);
+
+        if ((tmp = BN_CTX_get(ctx)) == NULL)
+                goto err;
+        if ((tmp_Z = BN_CTX_get(ctx)) == NULL)
+                goto err;
+
+        if ((prod_Z = calloc(num, sizeof *prod_Z)) == NULL)
+                goto err;
+        for (i = 0; i < num; i++) {
+                if ((prod_Z[i] = BN_CTX_get(ctx)) == NULL)
+                        goto err;
+        }
+
+        if (!BN_is_zero(points[0]->Z)) {
+                if (!bn_copy(prod_Z[0], points[0]->Z))
+                        goto err;
+        } else {
+                if (!BN_one(prod_Z[0]))
+                        goto err;
+        }
+
+        for (i = 1; i < num; i++) {
+                if (!BN_is_zero(points[i]->Z)) {
+                        if (!BN_mod_mul(prod_Z[i], prod_Z[i - 1], points[i]->Z,
+                            group->p, ctx))
+                                goto err;
+                } else {
+                        if (!bn_copy(prod_Z[i], prod_Z[i - 1]))
+                                goto err;
+                }
+        }
+
+        if (!BN_mod_inverse_nonct(tmp, prod_Z[num - 1], group->p, ctx)) {
+                ECerror(ERR_R_BN_LIB);
+                goto err;
+        }
+
+        for (i = num - 1; i > 0; i--) {
+                if (BN_is_zero(points[i]->Z))
+                        continue;
+
+                if (!BN_mod_mul(tmp_Z, prod_Z[i - 1], tmp, group->p, ctx))
+                        goto err;
+                if (!BN_mod_mul(tmp, tmp, points[i]->Z, group->p, ctx))
+                        goto err;
+                if (!bn_copy(points[i]->Z, tmp_Z))
+                        goto err;
+        }
+
+        for (i = 0; i < num; i++) {
+                EC_POINT *p = points[i];
+
+                if (BN_is_zero(p->Z))
+                        continue;
+
+                if (!BN_mod_mul(p->X, p->X, p->Z, group->p, ctx))
+                        goto err;
+                if (!BN_mod_mul(p->Y, p->Y, p->Z, group->p, ctx))
+                        goto err;
+
+                if (!BN_one(p->Z))
+                        goto err;
+        }
+
+        ret = 1;
+
+ err:
+        return ret;
+}
+#endif
+
+static int
+ec_mul(const EC_GROUP *group, EC_POINT *r, const BIGNUM *scalar, const EC_POINT *point,
+    BN_CTX *ctx)
+{
+        EC_POINT *rr;
+        int bits, i;
+        int ret = 0;
+
+        /* XXX - need constant time and blinding. */
+
+        if ((rr = EC_POINT_new(group)) == NULL)
+                goto err;
+
+        bits = BN_num_bits(scalar);
+
+        EC_POINT_copy(rr, point);
+
+        for (i = bits - 2; i >= 0; i--) {
+                if (!EC_POINT_dbl(group, rr, rr, ctx))
+                        goto err;
+                if (BN_is_bit_set(scalar, i)) {
+                        if (!EC_POINT_add(group, rr, rr, point, ctx))
+                                goto err;
+                }
+        }
+
+        EC_POINT_copy(r, rr);
+
+        ret = 1;
+
+ err:
+        EC_POINT_free(rr);
+
+        return ret;
+}
+
+static int
+ec_mul_single_ct(const EC_GROUP *group, EC_POINT *r, const BIGNUM *scalar,
+    const EC_POINT *point, BN_CTX *ctx)
+{
+        return ec_mul(group, r, scalar, point, ctx);
+}
+
+static int
+ec_mul_double_nonct(const EC_GROUP *group, EC_POINT *r, const BIGNUM *scalar1,
+    const EC_POINT *point1, const BIGNUM *scalar2, const EC_POINT *point2,
+    BN_CTX *ctx)
+{
+        return ec_wnaf_mul(group, r, scalar1, point1, scalar2, point2, ctx);
+}
+
+static const EC_METHOD ec_GFp_homogeneous_projective_method = {
+        .group_set_curve = ec_group_set_curve,
+        .group_get_curve = ec_group_get_curve,
+        .point_set_to_infinity = ec_point_set_to_infinity,
+        .point_is_at_infinity = ec_point_is_at_infinity,
+        .point_set_affine_coordinates = ec_point_set_affine_coordinates,
+        .point_get_affine_coordinates = ec_point_get_affine_coordinates,
+        .add = ec_point_add,
+        .dbl = ec_point_dbl,
+        .invert = ec_point_invert,
+        .point_is_on_curve = ec_point_is_on_curve,
+        .point_cmp = ec_point_cmp,
+        .points_make_affine = ec_points_make_affine,
+        .mul_single_ct = ec_mul_single_ct,
+        .mul_double_nonct = ec_mul_double_nonct,
+};
+
+const EC_METHOD *
+EC_GFp_homogeneous_projective_method(void)
+{
+        return &ec_GFp_homogeneous_projective_method;
+}
+LCRYPTO_ALIAS(EC_GFp_simple_method);