ipxe/src/crypto/weierstrass.c

/*
 * Copyright (C) 2024 Michael Brown <mbrown@fensystems.co.uk>.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 *
 * You can also choose to distribute this program under the terms of
 * the Unmodified Binary Distribution Licence (as given in the file
 * COPYING.UBDL), provided that you have satisfied its requirements.
 */

FILE_LICENCE ( GPL2_OR_LATER_OR_UBDL );

/** @file
 *
 * Weierstrass elliptic curves
 *
 * The implementation is based upon Algorithm 1 from "Complete
 * addition formulas for prime order elliptic curves" (Joost Renes,
 * Craig Costello, and Lejla Batina), available from
 *
 *   https://www.microsoft.com/en-us/research/wp-content/uploads/2016/06/complete-2.pdf
 *
 * The steps within the algorithm have been reordered and temporary
 * variables shuffled to reduce stack usage, and calculations are
 * carried out modulo small multiples of the field prime in order to
 * elide reductions after intermediate addition and subtraction
 * operations.
 *
 * The algorithm is encoded using a bytecode representation, since
 * this substantially reduces the code size compared to direct
 * implementation of the big integer operations.
 */

#include <errno.h>
#include <ipxe/weierstrass.h>

/** Big integer register names */
enum weierstrass_register {

	/*
	 * Read-only registers
	 */

	/* Curve constant "a" (for multiply), zero (for add/subtract) */
	WEIERSTRASS_a = 0,
	/* Curve constant "3b" */
	WEIERSTRASS_3b,
	/* Augend (x,y,z) co-ordinates */
	WEIERSTRASS_x1,
	WEIERSTRASS_y1,
	WEIERSTRASS_z1,
	/* Addend (x,y,z) co-ordinates */
	WEIERSTRASS_x2,
	WEIERSTRASS_y2,
	WEIERSTRASS_z2,

	/*
	 * Read-write registers
	 */

	/* Temporary working registers */
	WEIERSTRASS_Wt,
	WEIERSTRASS_Wxy,
	WEIERSTRASS_Wyz,
	WEIERSTRASS_Wzx,
	/* Low half of multiplication product */
	WEIERSTRASS_Wp,
	/* Result (x,y,z) co-ordinates */
	WEIERSTRASS_x3,
	WEIERSTRASS_y3,
	WEIERSTRASS_z3,

	/* Number of registers */
	WEIERSTRASS_NUM_REGISTERS = 16
};

/** Zero register (for add/subtract operations */
#define WEIERSTRASS_zero WEIERSTRASS_a

/** Construct big integer register index */
#define WEIERSTRASS_REGISTER( name ) _C2 ( WEIERSTRASS_, name )

/** Bytecode operation codes */
enum weierstrass_opcode {
	/** Subtract big integers (and add nothing)*/
	WEIERSTRASS_OP_SUB_0N = 0,
	/** Subtract big integers (and add 2N) */
	WEIERSTRASS_OP_SUB_2N = WEIERSTRASS_2N,
	/** Subtract big integers (and add 4N) */
	WEIERSTRASS_OP_SUB_4N = WEIERSTRASS_4N,
	/** Add big integers */
	WEIERSTRASS_OP_ADD,
	/** Multiply big integers (and perform Montgomery reduction) */
	WEIERSTRASS_OP_MUL,
};

/**
 * Define a bytecode operation
 *
 * @v opcode		Operation code
 * @v dest		Destination big integer register name
 * @v left		Left source big integer register name
 * @v right		Right source big integer register name
 */
#define WEIERSTRASS_OP( opcode, dest, left, right )	\
	( ( (opcode) << 12 ) |				\
	  ( WEIERSTRASS_REGISTER ( dest ) << 8 ) |	\
	  ( WEIERSTRASS_REGISTER ( left ) << 4 ) |	\
	  ( WEIERSTRASS_REGISTER ( right ) << 0 ) )

/** Extract bytecode operation code */
#define WEIERSTRASS_OPCODE( op ) ( ( (op) >> 12 ) & 0xf )

/** Extract destination big integer register */
#define WEIERSTRASS_DEST( op ) ( ( (op) >> 8 ) & 0xf )

/** Extract left source big integer register */
#define WEIERSTRASS_LEFT( op ) ( ( (op) >> 4 ) & 0xf )

/** Extract right source big integer register */
#define WEIERSTRASS_RIGHT( op ) ( ( (op) >> 0 ) & 0xf )

/** Define a three-argument addition operation */
#define WEIERSTRASS_ADD3( dest, augend, addend ) \
	WEIERSTRASS_OP ( WEIERSTRASS_OP_ADD, dest, augend, addend )

/** Define a two-argument addition operation */
#define WEIERSTRASS_ADD2( augend, addend ) \
	WEIERSTRASS_ADD3 ( augend, augend, addend )

/** Define a move operation */
#define WEIERSTRASS_MOV( dest, source ) \
	WEIERSTRASS_ADD3( dest, source, zero )

/** Define a three-argument subtraction operation */
#define WEIERSTRASS_SUB3( dest, minuend, subtrahend, multiple ) \
	WEIERSTRASS_OP ( _C2 ( WEIERSTRASS_OP_SUB_, multiple ),	  \
			 dest, minuend, subtrahend )

/** Define a two-argument subtraction operation */
#define WEIERSTRASS_SUB2( minuend, subtrahend, multiple ) \
	WEIERSTRASS_SUB3 ( minuend, minuend, subtrahend, multiple )

/** Define a stop operation */
#define WEIERSTRASS_STOP WEIERSTRASS_SUB2 ( zero, zero, 0N )

/** Define a three-argument multiplication operation */
#define WEIERSTRASS_MUL3( dest, multiplicand, multiplier ) \
	WEIERSTRASS_OP ( WEIERSTRASS_OP_MUL, dest, multiplicand, multiplier )

/** Define a two-argument multiplication operation */
#define WEIERSTRASS_MUL2( multiplicand, multiplier ) \
	WEIERSTRASS_MUL3 ( multiplicand, multiplicand, multiplier )

/**
 * Initialise curve
 *
 * @v curve		Weierstrass curve
 */
static void weierstrass_init ( struct weierstrass_curve *curve ) {
	unsigned int size = curve->size;
	bigint_t ( size ) __attribute__ (( may_alias )) *prime =
		( ( void * ) curve->prime[0] );
	bigint_t ( size ) __attribute__ (( may_alias )) *fermat =
		( ( void * ) curve->fermat );
	bigint_t ( size ) __attribute__ (( may_alias )) *square =
		( ( void * ) curve->square );
	bigint_t ( size ) __attribute__ (( may_alias )) *one =
		( ( void * ) curve->one );
	bigint_t ( size ) __attribute__ (( may_alias )) *a =
		( ( void * ) curve->a );
	bigint_t ( size ) __attribute__ (( may_alias )) *b3 =
		( ( void * ) curve->b3 );
	bigint_t ( size ) __attribute__ (( may_alias )) *mont =
		( ( void * ) curve->mont[0] );
	bigint_t ( size ) __attribute__ (( may_alias )) *temp =
		( ( void * ) curve->prime[1] );
	bigint_t ( size * 2 ) __attribute__ (( may_alias )) *product =
		( ( void * ) temp );
	bigint_t ( size ) __attribute__ (( may_alias )) *two =
		( ( void * ) temp );
	static const uint8_t one_raw[] = { 1 };
	static const uint8_t two_raw[] = { 2 };
	size_t len = curve->len;
	unsigned int i;

	/* Initialise field prime */
	bigint_init ( prime, curve->prime_raw, len );
	DBGC ( curve, "WEIERSTRASS %s   N = %s\n",
	       curve->name, bigint_ntoa ( prime ) );

	/* Calculate Montgomery constant R^2 mod N */
	bigint_reduce ( prime, square );
	DBGC ( curve, "WEIERSTRASS %s R^2 = %s mod N\n",
	       curve->name, bigint_ntoa ( square ) );

	/* Calculate constant "3b" */
	bigint_init ( b3, curve->b_raw, len );
	DBGC ( curve, "WEIERSTRASS %s   b = %s\n",
	       curve->name, bigint_ntoa ( b3 ) );
	bigint_copy ( b3, a );
	bigint_add ( b3, b3 );
	bigint_add ( a, b3 );

	/* Initialise "a" */
	bigint_init ( a, curve->a_raw, len );
	DBGC ( curve, "WEIERSTRASS %s   a = %s\n",
	       curve->name, bigint_ntoa ( a ) );

	/* Initialise "1" */
	bigint_init ( one, one_raw, sizeof ( one_raw ) );

	/* Convert relevant constants to Montgomery form
	 *
	 * We rely on the fact that the prime multiples have not yet
	 * been calculated, and so can be used as a temporary buffer.
	 */
	for ( i = 0 ; i < WEIERSTRASS_NUM_MONT ; i++ ) {
		static const char *names[] = { "  ", " a", "3b" };
		bigint_multiply ( &mont[i], square, product );
		bigint_montgomery ( prime, product, &mont[i] );
		DBGC ( curve, "WEIERSTRASS %s %sR = %s mod N\n",
		       curve->name, names[i], bigint_ntoa ( &mont[i] ) );
	}

	/* Calculate constant "N-2"
	 *
	 * We rely on the fact that the prime multiples have not yet
	 * been calculated, and so can be used as a temporary buffer.
	 */
	bigint_copy ( prime, fermat );
	bigint_init ( two, two_raw, sizeof ( two_raw ) );
	bigint_subtract ( two, fermat );
	DBGC ( curve, "WEIERSTRASS %s N-2 = %s\n",
	       curve->name, bigint_ntoa ( fermat ) );

	/* Calculate multiples of field prime */
	for ( i = 1 ; i < WEIERSTRASS_NUM_MULTIPLES ; i++ ) {
		bigint_copy ( &prime[ i - 1 ], &prime[i] );
		bigint_add ( &prime[i], &prime[i] );
		DBGC ( curve, "WEIERSTRASS %s  %dN = %s\n",
		       curve->name, ( 1 << i ), bigint_ntoa ( &prime[i] ) );
	}
}

/**
 * Execute bytecode instruction
 *
 * @v curve		Weierstrass curve
 * @v regs		Registers
 * @v size		Big integer size
 * @v op		Operation
 */
static void weierstrass_exec ( const struct weierstrass_curve *curve,
			       void **regs, unsigned int size,
			       unsigned int op ) {
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*prime = ( ( const void * ) curve->prime[0] );
	bigint_t ( size * 2 ) __attribute__ (( may_alias ))
		*product = regs[WEIERSTRASS_Wp];
	bigint_t ( size ) __attribute__ (( may_alias )) *dest;
	const bigint_t ( size ) __attribute__ (( may_alias )) *left;
	const bigint_t ( size ) __attribute__ (( may_alias )) *right;
	const bigint_t ( size ) __attribute__ (( may_alias )) *addend;
	const bigint_t ( size ) __attribute__ (( may_alias )) *subtrahend;
	unsigned int op_code;
	unsigned int op_dest;
	unsigned int op_left;
	unsigned int op_right;

	/* Decode instruction */
	op_code = WEIERSTRASS_OPCODE ( op );
	op_dest = WEIERSTRASS_DEST ( op );
	op_left = WEIERSTRASS_LEFT ( op );
	op_right = WEIERSTRASS_RIGHT ( op );
	dest = regs[op_dest];
	left = regs[op_left];
	right = regs[op_right];

	/* Check destination is a writable register */
	assert ( op_dest >= WEIERSTRASS_Wt );

	/* Handle multiplications */
	if ( op_code == WEIERSTRASS_OP_MUL ) {
		assert ( op_left != WEIERSTRASS_Wp );
		assert ( op_right != WEIERSTRASS_Wp );
		bigint_multiply ( left, right, product );
		bigint_montgomery_relaxed ( prime, product, dest );
		DBGCP ( curve, "WEIERSTRASS %s R%d := R%d x R%d = %s\n",
			curve->name, op_dest, op_left, op_right,
			bigint_ntoa ( dest ) );
		return;
	}

	/* Copy left source, if required */
	if ( op_dest != op_left )
		bigint_copy ( left, dest );

	/* Do nothing more if addend/subtrahend is zero */
	if ( ! op_right ) {
		DBGCP ( curve, "WEIERSTRASS %s R%d := R%d = %s\n",
			curve->name, op_dest, op_left, bigint_ntoa ( dest ) );
		return;
	}

	/* Determine addend and subtrahend */
	addend = NULL;
	subtrahend = NULL;
	if ( op_code == WEIERSTRASS_OP_ADD ) {
		DBGCP ( curve, "WEIERSTRASS %s R%d := R%d + R%d = ",
			curve->name, op_dest, op_left, op_right );
		addend = ( ( const void * ) right );
	} else {
		subtrahend = ( ( const void * ) right );
		if ( op_code > WEIERSTRASS_OP_SUB_0N ) {
			DBGCP ( curve, "WEIERSTRASS %s R%d := R%d - R%d + "
				"%dN = ", curve->name, op_dest, op_left,
				op_right, ( 1 << op_code ) );
			addend = ( ( const void * ) curve->prime[op_code] );
		} else {
			DBGCP ( curve, "WEIERSTRASS %s R%d := R%d - R%d = ",
				curve->name, op_dest, op_left, op_right );
		}
	}

	/* Perform addition and subtraction */
	if ( addend )
		bigint_add ( addend, dest );
	if ( subtrahend )
		bigint_subtract ( subtrahend, dest );
	DBGCP ( curve, "%s\n", bigint_ntoa ( dest ) );
}

/**
 * Add points on curve
 *
 * @v curve		Weierstrass curve
 * @v augend0		Element 0 of point (x1,y1,z1) to be added
 * @v addend0		Element 0 of point (x2,y2,z2) to be added
 * @v result0		Element 0 of point (x3,y3,z3) to hold result
 *
 * Points are represented in projective coordinates, with all values
 * in Montgomery form and in the range [0,4N) where N is the field
 * prime.
 *
 * The augend may have the same value as the addend (i.e. this routine
 * may be used to perform point doubling as well as point addition),
 * and either or both may be the point at infinity.
 *
 * The result may overlap either input, since the inputs are fully
 * consumed before the result is written.
 */
static void weierstrass_add_raw ( const struct weierstrass_curve *curve,
				  const bigint_element_t *augend0,
				  const bigint_element_t *addend0,
				  bigint_element_t *result0 ) {
	unsigned int size = curve->size;
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*prime = ( ( const void * ) curve->prime[0] );
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*a = ( ( const void * ) curve->a );
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*b3 = ( ( const void * ) curve->b3 );
	const weierstrass_t ( size ) __attribute__ (( may_alias ))
		*augend = ( ( const void * ) augend0 );
	const weierstrass_t ( size ) __attribute__ (( may_alias ))
		*addend = ( ( const void * ) addend0 );
	weierstrass_t ( size ) __attribute__ (( may_alias ))
		*result = ( ( void * ) result0 );
	struct {
		bigint_t ( size ) Wt;
		bigint_t ( size ) Wxy;
		bigint_t ( size ) Wyz;
		bigint_t ( size ) Wzx;
		bigint_t ( size * 2 ) Wp;
	} temp;
	void *regs[WEIERSTRASS_NUM_REGISTERS];
	unsigned int schedule;
	const uint16_t *op;
	unsigned int i;

	/* On entry, we assume that x1, x2, y1, y2, z1, z2 are all in
	 * the range [0,4N).  Additions will extend the range.
	 * Subtractions will extend the range (and require an addition
	 * of a suitable multiple of the modulus to ensure that the
	 * result is a positive value).  Relaxed Montgomery
	 * multiplications will reduce the range to [0,2N).  The
	 * outputs x3, y3, z3 will be in the range [0,4N) and
	 * therefore usable as subsequent inputs.
	 */
	static const uint16_t ops[] = {
		/* [Wxy] Qxy = (x1+y1)*(x2+y2)			    (mod 2N) */
		WEIERSTRASS_ADD3 ( Wt, x1, y1 ),
		WEIERSTRASS_ADD3 ( Wxy, x2, y2 ),
		WEIERSTRASS_MUL2 ( Wxy, Wt ),
		/* [Wyz] Qyz = (y1+z1)*(y2+z2)			    (mod 2N) */
		WEIERSTRASS_ADD3 ( Wt, y1, z1 ),
		WEIERSTRASS_ADD3 ( Wyz, y2, z2 ),
		WEIERSTRASS_MUL2 ( Wyz, Wt ),
		/* [Wzx] Qzx = (z1+x1)*(z2+x2)			    (mod 2N) */
		WEIERSTRASS_ADD3 ( Wt, z1, x1 ),
		WEIERSTRASS_ADD3 ( Wzx, z2, x2 ),
		WEIERSTRASS_MUL2 ( Wzx, Wt ),
		/* [x3] Px = x1*x2				    (mod 2N) */
		WEIERSTRASS_MUL3 ( x3, x1, x2 ),
		/* [y3] Py = y1*y2				    (mod 2N) */
		WEIERSTRASS_MUL3 ( y3, y1, y2 ),
		/* [z3] Pz = z1*z2				    (mod 2N) */
		WEIERSTRASS_MUL3 ( z3, z1, z2 ),
		/* [Wxy] Rxy = Qxy - Px - Py			    (mod 6N)
		 *	     = (x1+y1)*(x2+y2) - x1*x2 - y1*y2	    (mod 6N)
		 *	     = x1*y2 + x2*y1			    (mod 6N)
		 */
		WEIERSTRASS_SUB2 ( Wxy, x3, 0N ),
		WEIERSTRASS_SUB2 ( Wxy, y3, 4N ),
		/* [Wyz] Ryz = Qyz - Py - Pz			    (mod 6N)
		 *	     = (y1+z1)*(y2+z2) - y1*y2 - z1*z2	    (mod 6N)
		 *	     = y1*z2 + y2*z1			    (mod 6N)
		 */
		WEIERSTRASS_SUB2 ( Wyz, y3, 0N ),
		WEIERSTRASS_SUB2 ( Wyz, z3, 4N ),
		/* [Wzx] Rzx = Qzx - Pz - Px			    (mod 6N)
		 *	     = (z1+x1)*(z2+x2) - z1*z2 - x1*x2	    (mod 6N)
		 *	     = x1*z2 + x2*z1			    (mod 6N)
		 */
		WEIERSTRASS_SUB2 ( Wzx, z3, 0N ),
		WEIERSTRASS_SUB2 ( Wzx, x3, 4N ),
		/* [Wt] aRzx = a * Rzx				    (mod 2N)
		 *	     = a * (x1*z2 + x2*z1)		    (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( Wt, a, Wzx ),
		/* [Wp] 3bPz = 3b * Pz				    (mod 2N)
		 *	     = 3b*z1*z2				    (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( Wp, 3b, z3 ),
		/* [Wp] Sy = aRzx + 3bPz			    (mod 4N)
		 *	   = a*(x1*z2 + x2*z1) + 3b*z1*z2	    (mod 4N)
		 */
		WEIERSTRASS_ADD2 ( Wp, Wt ),
		/* [Wt] Syz = Py + Sy				    (mod 6N)
		 *	    = y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2  (mod 6N)
		 */
		WEIERSTRASS_ADD3 ( Wt, y3, Wp ),
		/* [y3] Sxy = Py - Sy				    (mod 6N)
		 *	    = y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2  (mod 6N)
		 */
		WEIERSTRASS_SUB2 ( y3, Wp, 4N ),
		/* [z3] aPz = a * Pz				    (mod 2N)
		 *	    = a * z1*z2				    (mod 2N)
		 */
		WEIERSTRASS_MUL2 ( z3, a ),
		/* [Wzx] 3bRzx = 3b * Rzx			    (mod 2N)
		 *	       = 3b * (x1*z2 + x2*z1)		    (mod 2N)
		 */
		WEIERSTRASS_MUL2 ( Wzx, 3b ),
		/* [x3] aPzx' = Px - aPz			    (mod 4N)
		 *	      = x1*x2 - a*z1*z2			    (mod 4N)
		 */
		WEIERSTRASS_SUB2 ( x3, z3, 2N ),
		/* [Wp] Szx = a * aPzx'				    (mod 2N)
		 *	    = a * (x1*x2 - a*z1*z2)		    (mod 2N)
		 * 	    = a*x1*x2 - (a^2)*z1*z2		    (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( Wp, a, x3 ),
		/* [x3] Px = aPzx' + aPz			    (mod 6N)
		 *	   = x1*x2 - a*z1*z2 + a*z1*z2		    (mod 6N)
		 *	   = x1*x2				    (mod 6N)
		 */
		WEIERSTRASS_ADD2 ( x3, z3 ),
		/* [Wzx] Tzx = 3bRzx + Szx			    (mod 4N)
		 *	     = a*x1*x2 + 3b*(x1*z2 + x2*z1) -
		 *	       (a^2)*z1*z2			    (mod 4N)
		 */
		WEIERSTRASS_ADD2 ( Wzx, Wp ),
		/* [z3] aPzx = Px + aPz				    (mod 8N)
		 *	     = x1*x2 + a*z1*z2			    (mod 8N)
		 */
		WEIERSTRASS_ADD2 ( z3, x3 ),
		/* [x3] 2Px = Px + Px				    (mod 12N)
		 *          = 2*x1*x2				    (mod 12N)
		 */
		WEIERSTRASS_ADD2 ( x3, x3 ),
		/* [x3] Tyz = 2Px + aPzx			    (mod 20N)
		 *          = 2*x1*x2 + x1*x2 + a*z1*z2		    (mod 20N)
		 *	    = 3*x1*x2 + a*z1*z2			    (mod 20N)
		 */
		WEIERSTRASS_ADD2 ( x3, z3 ),
		/* [z3] Syz = Syz				    (mod 6N)
		 *	    = y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2  (mod 6N)
		 */
		WEIERSTRASS_MOV ( z3, Wt ),
		/* [Wt] Tyz = Tyz				    (mod 20N)
		 *	    = 3*x1*x2 + a*z1*z2			    (mod 20N)
		 */
		WEIERSTRASS_MOV ( Wt, x3 ),
		/* [x3] Ux = Rxy * Sxy				    (mod 2N)
		 *	   = (x1*y2 + x2*y1) *
		 *	     (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2) (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( x3, Wxy, y3 ),
		/* [y3] Uy = Syz * Sxy				    (mod 2N)
		 *	   = (y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2) *
		 *	     (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2) (mod 2N)
		 */
		WEIERSTRASS_MUL2 ( y3, z3 ),
		/* [z3] Uz = Ryz * Syz				    (mod 2N)
		 *	   = (y1*z2 + y2*z1) *
		 *	     (y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2) (mod 2N)
		 */
		WEIERSTRASS_MUL2 ( z3, Wyz ),
		/* [Wp] Vx = Ryz * Tzx				    (mod 2N)
		 *	   = (y1*z2 + y2*z1) *
		 *	     (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2)
		 *						    (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( Wp, Wyz, Wzx ),
		/* [x3] x3 = Ux - Vx				    (mod 4N)
		 *	   = ((x1*y2 + x2*y1) *
		 *	      (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2)) -
		 *	     ((y1*z2 + y2*z1) *
		 *	      (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2))
		 *						    (mod 4N)
		 */
		WEIERSTRASS_SUB2 ( x3, Wp, 2N ),
		/* [Wp] Vy = Tyz * Tzx				    (mod 2N)
		 *	   = (3*x1*x2 + a*z1*z2) *
		 *	     (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2)
		 *						    (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( Wp, Wt, Wzx ),
		/* [y3] y3 = Vy + Uy				    (mod 4N)
		 *         = ((3*x1*x2 + a*z1*z2) *
		 *	      (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2)) +
		 *	     ((y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2) *
		 *	      (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2))
		 *						    (mod 4N)
		 */
		WEIERSTRASS_ADD2 ( y3, Wp ),
		/* [Wp] Vz = Rxy * Tyz				    (mod 2N)
		 *	   = (x1*y2 + x2*y1) * (3*x1*x2 + a*z1*z2)  (mod 2N)
		 */
		WEIERSTRASS_MUL3 ( Wp, Wxy, Wt ),
		/* [z3] z3 = Uz + Vz				    (mod 4N)
		 *	   = ((y1*z2 + y2*z1) *
		 *	      (y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2)) +
		 *	     ((x1*y2 + x2*y1) * (3*x1*x2 + a*z1*z2))
		 *						    (mod 4N)
		 */
		WEIERSTRASS_ADD2 ( z3, Wp ),
		/* Stop */
		WEIERSTRASS_STOP
	};

	/* Initialise register list */
	regs[WEIERSTRASS_a] = ( ( void * ) a );
	regs[WEIERSTRASS_3b] = ( ( void * ) b3 );
	regs[WEIERSTRASS_x1] = ( ( void * ) &augend->x );
	regs[WEIERSTRASS_x2] = ( ( void * ) &addend->x );
	regs[WEIERSTRASS_x3] = ( ( void * ) &result->x );
	regs[WEIERSTRASS_Wt] = &temp;
	schedule = ( ( ( 1 << WEIERSTRASS_NUM_REGISTERS ) - 1 )
		     - ( 1 << WEIERSTRASS_a )
		     - ( 1 << WEIERSTRASS_3b )
		     - ( 1 << WEIERSTRASS_x1 )
		     - ( 1 << WEIERSTRASS_x2 )
		     - ( 1 << WEIERSTRASS_x3 )
		     - ( 1 << WEIERSTRASS_Wt ) );
	for ( i = 0 ; schedule ; i++, schedule >>= 1 ) {
		if ( schedule & 1 )
			regs[i] = ( regs[ i - 1 ] + sizeof ( *prime ) );
	}
	DBGC2 ( curve, "WEIERSTRASS %s augend (%s,",
		curve->name, bigint_ntoa ( &augend->x ) );
	DBGC2 ( curve, "%s,", bigint_ntoa ( &augend->y ) );
	DBGC2 ( curve, "%s)\n", bigint_ntoa ( &augend->z ) );
	DBGC2 ( curve, "WEIERSTRASS %s addend (%s,",
		curve->name, bigint_ntoa ( &addend->x ) );
	DBGC2 ( curve, "%s,", bigint_ntoa ( &addend->y ) );
	DBGC2 ( curve, "%s)\n", bigint_ntoa ( &addend->z ) );

	/* Sanity checks */
	assert ( regs[WEIERSTRASS_a] == a );
	assert ( regs[WEIERSTRASS_3b] == b3 );
	assert ( regs[WEIERSTRASS_x1] == &augend->x );
	assert ( regs[WEIERSTRASS_y1] == &augend->y );
	assert ( regs[WEIERSTRASS_z1] == &augend->z );
	assert ( regs[WEIERSTRASS_x2] == &addend->x );
	assert ( regs[WEIERSTRASS_y2] == &addend->y );
	assert ( regs[WEIERSTRASS_z2] == &addend->z );
	assert ( regs[WEIERSTRASS_x3] == &result->x );
	assert ( regs[WEIERSTRASS_y3] == &result->y );
	assert ( regs[WEIERSTRASS_z3] == &result->z );
	assert ( regs[WEIERSTRASS_Wt] == &temp.Wt );
	assert ( regs[WEIERSTRASS_Wxy] == &temp.Wxy );
	assert ( regs[WEIERSTRASS_Wyz] == &temp.Wyz );
	assert ( regs[WEIERSTRASS_Wzx] == &temp.Wzx );
	assert ( regs[WEIERSTRASS_Wp] == &temp.Wp );

	/* Execute bytecode instruction sequence */
	for ( op = ops ; *op != WEIERSTRASS_STOP ; op++ )
		weierstrass_exec ( curve, regs, size, *op );
	DBGC2 ( curve, "WEIERSTRASS %s result (%s,",
		curve->name, bigint_ntoa ( &result->x ) );
	DBGC2 ( curve, "%s,", bigint_ntoa ( &result->y ) );
	DBGC2 ( curve, "%s)\n", bigint_ntoa ( &result->z ) );
}

/**
 * Add points on curve
 *
 * @v curve		Weierstrass curve
 * @v augend		Point (x1,y1,z1) to be added
 * @v addend		Point (x2,y2,z2) to be added
 * @v result0		Point (x3,y3,z3) to hold result
 */
#define weierstrass_add( curve, augend, addend, result ) do {		\
	weierstrass_add_raw ( (curve), (augend)->all.element,		\
			      (addend)->all.element,			\
			      (result)->all.element );			\
	} while ( 0 )

/**
 * Add points on curve as part of a Montgomery ladder
 *
 * @v operand		Element 0 of first input operand (may overlap result)
 * @v result		Element 0 of second input operand and result
 * @v size		Number of elements in operands and result
 * @v ctx		Operation context
 * @v tmp		Temporary working space (not used)
 */
static void weierstrass_add_ladder ( const bigint_element_t *operand0,
				     bigint_element_t *result0,
				     unsigned int size, const void *ctx,
				     void *tmp __unused ) {
	const struct weierstrass_curve *curve = ctx;
	const weierstrass_t ( curve->size ) __attribute__ (( may_alias ))
		*operand = ( ( const void * ) operand0 );
	weierstrass_t ( curve->size ) __attribute__ (( may_alias ))
		*result = ( ( void * ) result0 );

	/* Add curve points */
	assert ( size == bigint_size ( &operand->all ) );
	assert ( size == bigint_size ( &result->all ) );
	weierstrass_add ( curve, operand, result, result );
}

/**
 * Verify point is on curve
 *
 * @v curve		Weierstrass curve
 * @v point0		Element 0 of point (x,y,z) to be verified
 * @ret rc		Return status code
 *
 * As with point addition, points are represented in projective
 * coordinates, with all values in Montgomery form and in the range
 * [0,4N) where N is the field prime.
 */
static int weierstrass_verify_raw ( const struct weierstrass_curve *curve,
				    const bigint_element_t *point0 ) {
	unsigned int size = curve->size;
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*prime = ( ( const void * ) curve->prime[0] );
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*a = ( ( const void * ) curve->a );
	const bigint_t ( size ) __attribute__ (( may_alias ))
		*b3 = ( ( const void * ) curve->b3 );
	const weierstrass_t ( size ) __attribute__ (( may_alias ))
		*point = ( ( const void * ) point0 );
	struct {
		bigint_t ( size ) Wt;
		bigint_t ( size * 2 ) Wp;
	} temp;
	void *regs[WEIERSTRASS_NUM_REGISTERS];
	const uint16_t *op;

	/* Calculate 3*(x^3 + a*x + b - y^2) */
	static const uint16_t ops[] = {
		/* [Wt] Tx = x^2				    (mod 2N) */
		WEIERSTRASS_MUL3 ( Wt, x1, x1 ),
		/* [Wt] Txa = Tx + a				    (mod 3N)
		 *	    = x^2 + a				    (mod 3N)
		 */
		WEIERSTRASS_MOV ( Wp, a ),
		WEIERSTRASS_ADD2 ( Wt, Wp ),
		/* [Wt] Txax = Txa * x				    (mod 2N)
		 *	     = (x^2 + a)*x			    (mod 2N)
		 *	     = x^3 + a*x			    (mod 2N)
		 */
		WEIERSTRASS_MUL2 ( Wt, x1 ),
		/* [Wp] Ty = y^2				    (mod 2N) */
		WEIERSTRASS_MUL3 ( Wp, y1, y1 ),
		/* [Wt] Txaxy = Txax - Ty			    (mod 4N)
		 *	      = x^3 + a*x - y^2			    (mod 4N)
		 */
		WEIERSTRASS_SUB2 ( Wt, Wp, 2N ),
		/* [Wp] 2Txaxy = Txaxy + Txaxy			    (mod 8N)
		 *	       = 2*(x^3 + a*x - y^2)		    (mod 8N)
		 */
		WEIERSTRASS_ADD3 ( Wp, Wt, Wt ),
		/* [Wt] 3Txaxy = 2Txaxy + Txaxy			    (mod 12N)
		 *	       = 3*(x^3 + a*x - y^2)		    (mod 12N)
		 */
		WEIERSTRASS_ADD2 ( Wt, Wp ),
		/* [Wt] 3Txaxyb = 3Txaxy + 3b			    (mod 13N)
		 *		= 3*(x^3 + a*x - y^2) + 3b	    (mod 13N)
		 *		= 3*(x^3 + a*x + b - y^2)	    (mod 13N)
		 */
		WEIERSTRASS_ADD2 ( Wt, 3b ),
		/* Stop */
		WEIERSTRASS_STOP
	};

	/* Initialise register list */
	regs[WEIERSTRASS_a] = ( ( void * ) a );
	regs[WEIERSTRASS_3b] = ( ( void * ) b3 );
	regs[WEIERSTRASS_x1] = ( ( void * ) &point->x );
	regs[WEIERSTRASS_y1] = ( ( void * ) &point->y );
	regs[WEIERSTRASS_Wt] = &temp.Wt;
	regs[WEIERSTRASS_Wp] = &temp.Wp;

	/* Execute bytecode instruction sequence */
	for ( op = ops ; *op != WEIERSTRASS_STOP ; op++ )
		weierstrass_exec ( curve, regs, size, *op );

	/* Check that result is zero (modulo the field prime) */
	bigint_grow ( &temp.Wt, &temp.Wp );
	bigint_montgomery ( prime, &temp.Wp, &temp.Wt );
	if ( ! bigint_is_zero ( &temp.Wt ) ) {
		DBGC ( curve, "WEIERSTRASS %s base point is not on curve\n",
		       curve->name );
		return -EINVAL;
	}

	return 0;
}

/**
 * Verify point is on curve
 *
 * @v curve		Weierstrass curve
 * @v point		Point (x,y,z) to be verified
 * @ret rc		Return status code
 */
#define weierstrass_verify( curve, point ) ( {				\
	weierstrass_verify_raw ( (curve), (point)->all.element );	\
	} )

/**
 * Multiply curve point by scalar
 *
 * @v curve		Weierstrass curve
 * @v base		Base point (or NULL to use generator)
 * @v scalar		Scalar multiple
 * @v result		Result point to fill in
 * @ret rc		Return status code
 */
int weierstrass_multiply ( struct weierstrass_curve *curve, const void *base,
			   const void *scalar, void *result ) {
	unsigned int size = curve->size;
	size_t len = curve->len;
	const bigint_t ( size ) __attribute__ (( may_alias )) *prime =
		( ( const void * ) curve->prime[0] );
	const bigint_t ( size ) __attribute__ (( may_alias )) *prime2 =
		( ( const void * ) curve->prime[WEIERSTRASS_2N] );
	const bigint_t ( size ) __attribute__ (( may_alias )) *fermat =
		( ( const void * ) curve->fermat );
	const bigint_t ( size ) __attribute__ (( may_alias )) *square =
		( ( const void * ) curve->square );
	const bigint_t ( size ) __attribute__ (( may_alias )) *one =
		( ( const void * ) curve->one );
	struct {
		union {
			weierstrass_t ( size ) result;
			bigint_t ( size * 2 ) product_in;
		};
		union {
			weierstrass_t ( size ) multiple;
			bigint_t ( size * 2 ) product_out;
		};
		bigint_t ( bigint_required_size ( len ) ) scalar;
	} temp;
	size_t offset;
	unsigned int i;
	int rc;

	/* Initialise curve, if not already done
	 *
	 * The least significant element of the field prime must be
	 * odd, and so the least significant element of the
	 * (initialised) first multiple of the field prime must be
	 * non-zero.
	 */
	if ( ! prime2->element[0] )
		weierstrass_init ( curve );

	/* Use generator if applicable */
	if ( ! base )
		base = curve->base;

	/* Convert input to projective coordinates in Montgomery form */
	DBGC ( curve, "WEIERSTRASS %s base (", curve->name );
	for ( i = 0, offset = 0 ; i < WEIERSTRASS_AXES ; i++, offset += len ) {
		bigint_init ( &temp.multiple.axis[i], ( base + offset ), len );
		DBGC ( curve, "%s%s", ( i ? "," : "" ),
		       bigint_ntoa ( &temp.multiple.axis[i] ) );
		bigint_multiply ( &temp.multiple.axis[i], square,
				  &temp.product_in );
		bigint_montgomery_relaxed ( prime, &temp.product_in,
					    &temp.multiple.axis[i] );
	}
	bigint_copy ( one, &temp.multiple.z );
	DBGC ( curve, ")\n" );

	/* Verify point is on curve */
	if ( ( rc = weierstrass_verify ( curve, &temp.multiple ) ) != 0 )
		return rc;

	/* Construct identity element (the point at infinity) */
	memset ( &temp.result, 0, sizeof ( temp.result ) );
	bigint_copy ( one, &temp.result.y );

	/* Initialise scalar */
	bigint_init ( &temp.scalar, scalar, len );
	DBGC ( curve, "WEIERSTRASS %s scalar %s\n",
	       curve->name, bigint_ntoa ( &temp.scalar ) );

	/* Perform multiplication via Montgomery ladder */
	bigint_ladder ( &temp.result.all, &temp.multiple.all, &temp.scalar,
			weierstrass_add_ladder, curve, NULL );

	/* Invert result Z co-ordinate (via Fermat's little theorem) */
	bigint_copy ( one, &temp.multiple.z );
	bigint_ladder ( &temp.multiple.z, &temp.result.z, fermat,
			bigint_mod_exp_ladder, prime, &temp.product_out );

	/* Convert result back to affine co-ordinates */
	DBGC ( curve, "WEIERSTRASS %s result (", curve->name );
	for ( i = 0, offset = 0 ; i < WEIERSTRASS_AXES ; i++, offset += len ) {
		bigint_multiply ( &temp.result.axis[i], &temp.multiple.z,
				  &temp.product_out );
		bigint_montgomery_relaxed ( prime, &temp.product_out,
					    &temp.result.axis[i] );
		bigint_grow ( &temp.result.axis[i], &temp.product_out );
		bigint_montgomery ( prime, &temp.product_out,
				    &temp.result.axis[i] );
		DBGC ( curve, "%s%s", ( i ? "," : "" ),
		       bigint_ntoa ( &temp.result.axis[i] ) );
		bigint_done ( &temp.result.axis[i], ( result + offset ), len );
	}
	DBGC ( curve, ")\n" );

	return 0;
}