/*	$OpenBSD: ecp_hp_methods.c,v 1.3 2025/06/01 03:23:33 tb 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:
	BN_CTX_end(ctx);
	free(prod_Z);

	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;
}