[crypto] Add support for Weierstrass elliptic curve point multiplication

The NIST elliptic curves are Weierstrass curves and have the form

  y^2 = x^3 + ax + b

with each curve defined by its field prime, the constants "a" and "b",
and a generator base point.

Implement a constant-time algorithm for point addition, based upon
Algorithm 1 from "Complete addition formulas for prime order elliptic
curves" (Joost Renes, Craig Costello, and Lejla Batina), and use this
as a Montgomery ladder commutative operation to perform constant-time
point multiplication.

The code for point addition is implemented using a custom bytecode
interpreter with 16-bit instructions, since this results in
substantially smaller code than compiling the somewhat lengthy
sequence of arithmetic operations directly.  Values are calculated
modulo small multiples of the field prime in order to allow for the
use of relaxed Montgomery reduction.

Signed-off-by: Michael Brown <mcb30@ipxe.org>

3 files changed, 1044 insertions, 0 deletions

diff --git a/src/crypto/weierstrass.c b/src/crypto/weierstrass.c
new file mode 100644
index 000000000..be3909542
--- /dev/null
+++ b/src/crypto/weierstrass.c
@@ -0,0 +1,877 @@
+/*
+ * 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 )) *prime_double =
+		( ( void * ) prime );
+	bigint_t ( size * 2 ) __attribute__ (( may_alias )) *square_double =
+		( ( void * ) square );
+	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
+	 *
+	 * We rely on the fact that the subsequent big integers in the
+	 * cache (i.e. the first prime multiple, and the constant "1")
+	 * have not yet been written to, and so can be treated as
+	 * being the (zero) upper halves required to hold the
+	 * double-width value R^2.
+	 */
+	bigint_reduce_supremum ( prime_double, square_double );
+	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;
+}
diff --git a/src/include/ipxe/errfile.h b/src/include/ipxe/errfile.h
index 48b35c537..b826a4a6f 100644
--- a/src/include/ipxe/errfile.h
+++ b/src/include/ipxe/errfile.h
@@ -426,6 +426,7 @@ FILE_LICENCE ( GPL2_OR_LATER_OR_UBDL );
 #define ERRFILE_form_ui		      ( ERRFILE_OTHER | 0x00630000 )
 #define ERRFILE_usb_cmd		      ( ERRFILE_OTHER | 0x00640000 )
 #define ERRFILE_usb_settings	      ( ERRFILE_OTHER | 0x00650000 )
+#define ERRFILE_weierstrass	      ( ERRFILE_OTHER | 0x00660000 )
 
 /** @} */
 
diff --git a/src/include/ipxe/weierstrass.h b/src/include/ipxe/weierstrass.h
new file mode 100644
index 000000000..fb09fa6fa
--- /dev/null
+++ b/src/include/ipxe/weierstrass.h
@@ -0,0 +1,166 @@
+#ifndef _IPXE_WEIERSTRASS_H
+#define _IPXE_WEIERSTRASS_H
+
+/** @file
+ *
+ * Weierstrass elliptic curves
+ *
+ */
+
+FILE_LICENCE ( GPL2_OR_LATER_OR_UBDL );
+
+#include <ipxe/bigint.h>
+#include <ipxe/crypto.h>
+
+/** Number of axes in Weierstrass curve point representation */
+#define WEIERSTRASS_AXES 2
+
+/**
+ * Maximum multiple of field prime encountered during calculations
+ *
+ * Calculations are performed using values modulo a small multiple of
+ * the field prime, rather than modulo the field prime itself.  This
+ * allows explicit reductions after additions, subtractions, and
+ * relaxed Montgomery multiplications to be omitted entirely, provided
+ * that we keep careful track of the field prime multiple for each
+ * intermediate value.
+ *
+ * Relaxed Montgomery multiplication will produce a result in the
+ * range t < (1+m/k)N, where m is this maximum multiple of the field
+ * prime, and k is the constant in R > kN representing the leading
+ * zero padding in the big integer representation of the field prime.
+ * We choose to set k=m so that multiplications will always produce a
+ * result in the range t < 2N.
+ *
+ * This is expressed as the base-two logarithm of the multiple
+ * (rounded up), to simplify compile-time calculations.
+ */
+#define WEIERSTRASS_MAX_MULTIPLE_LOG2 5 /* maximum reached is mod 20N */
+
+/**
+ * Determine number of elements in scalar values for a Weierstrass curve
+ *
+ * @v len		Length of field prime, in bytes
+ * @ret size		Number of elements
+ */
+#define weierstrass_size( len )						\
+	bigint_required_size ( (len) +					\
+			       ( ( WEIERSTRASS_MAX_MULTIPLE_LOG2 + 7 )	\
+				 / 8 ) )
+
+/**
+ * Define a Weierstrass projective co-ordinate type
+ *
+ * @v size		Number of elements in scalar values
+ * @ret weierstrass_t	Projective co-ordinate type
+ */
+#define weierstrass_t( size )						\
+	union {								\
+		bigint_t ( size ) axis[3];				\
+		struct {						\
+			bigint_t ( size ) x;				\
+			bigint_t ( size ) y;				\
+			bigint_t ( size ) z;				\
+		};							\
+		bigint_t ( size * 3 ) all;				\
+	}
+
+/** Indexes for stored multiples of the field prime */
+enum weierstrass_multiple {
+	WEIERSTRASS_N = 0,
+	WEIERSTRASS_2N,
+	WEIERSTRASS_4N,
+	WEIERSTRASS_NUM_MULTIPLES
+};
+
+/** Number of cached in Montgomery form for each Weierstrass curve */
+#define WEIERSTRASS_NUM_MONT 3
+
+/** Number of cached big integers for each Weierstrass curve */
+#define WEIERSTRASS_NUM_CACHED						\
+	( WEIERSTRASS_NUM_MULTIPLES +					\
+	  1 /* fermat */ + 1 /* mont */ +				\
+	  WEIERSTRASS_NUM_MONT )
+
+/**
+ * A Weierstrass elliptic curve
+ *
+ * This is an elliptic curve y^2 = x^3 + ax + b
+ */
+struct weierstrass_curve {
+	/** Number of elements in scalar values */
+	const unsigned int size;
+	/** Curve name */
+	const char *name;
+	/** Length of raw scalar values */
+	size_t len;
+	/** Field prime */
+	const uint8_t *prime_raw;
+	/** Constant "a" */
+	const uint8_t *a_raw;
+	/** Constant "b" */
+	const uint8_t *b_raw;
+	/** Base point */
+	const uint8_t *base;
+
+	/** Cached field prime "N" (and multiples thereof) */
+	bigint_element_t *prime[WEIERSTRASS_NUM_CACHED];
+	/** Cached constant "N-2" (for Fermat's little theorem) */
+	bigint_element_t *fermat;
+	/** Cached Montgomery constant (R^2 mod N) */
+	bigint_element_t *square;
+	/** Cached constants in Montgomery form */
+	union {
+		struct {
+			/** Cached constant "1", in Montgomery form */
+			bigint_element_t *one;
+			/** Cached constant "a", in Montgomery form */
+			bigint_element_t *a;
+			/** Cached constant "3b", in Montgomery form */
+			bigint_element_t *b3;
+		};
+		bigint_element_t *mont[WEIERSTRASS_NUM_MONT];
+	};
+};
+
+extern int weierstrass_multiply ( struct weierstrass_curve *curve,
+				  const void *base, const void *scalar,
+				  void *result );
+
+/** Define a Weierstrass curve */
+#define WEIERSTRASS_CURVE( _name, _curve, _len, _prime, _a, _b, _base )	\
+	static bigint_t ( weierstrass_size(_len) )			\
+		_name ## _cache[WEIERSTRASS_NUM_CACHED];		\
+	static struct weierstrass_curve _name ## _weierstrass = {	\
+		.size = weierstrass_size(_len),				\
+		.name = #_name,						\
+		.len = (_len),						\
+		.prime_raw = (_prime),					\
+		.a_raw = (_a),						\
+		.b_raw = (_b),						\
+		.base = (_base),					\
+		.prime = {						\
+			(_name ## _cache)[0].element,			\
+			(_name ## _cache)[1].element,			\
+			(_name ## _cache)[2].element,			\
+		},							\
+		.fermat = (_name ## _cache)[3].element,			\
+		.square = (_name ## _cache)[4].element,			\
+		.one = (_name ## _cache)[5].element,			\
+		.a = (_name ## _cache)[6].element,			\
+		.b3 = (_name ## _cache)[7].element,			\
+	};								\
+	static int _name ## _multiply ( const void *base,		\
+					const void *scalar,		\
+					void *result ) {		\
+		return weierstrass_multiply ( &_name ## _weierstrass,	\
+					      base, scalar, result );	\
+	}								\
+	struct elliptic_curve _curve = {				\
+		.name = #_name,						\
+		.pointsize = ( WEIERSTRASS_AXES * (_len) ),		\
+		.keysize = (_len),					\
+		.multiply = _name ## _multiply,				\
+	}
+
+#endif /* _IPXE_WEIERSTRASS_H */