[crypto] Add X25519 key exchange algorithm

Add an implementation of the X25519 key exchange algorithm as defined
in RFC7748.

This implementation is inspired by and partially based upon the paper
"Implementing Curve25519/X25519: A Tutorial on Elliptic Curve
Cryptography" by Martin Kleppmann, available for download from
https://www.cl.cam.ac.uk/teaching/2122/Crypto/curve25519.pdf

The underlying modular addition, subtraction, and multiplication
operations are completely redesigned for substantially improved
efficiency compared to the TweetNaCl implementation studied in that
paper (approximately 5x-10x faster and with 70% less memory usage).

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

src/crypto/x25519.c

@@ -0,0 +1,808 @@
+/*
+ * 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
+ *
+ * X25519 key exchange
+ *
+ * This implementation is inspired by and partially based upon the
+ * paper "Implementing Curve25519/X25519: A Tutorial on Elliptic Curve
+ * Cryptography" by Martin Kleppmann, available for download from
+ * https://www.cl.cam.ac.uk/teaching/2122/Crypto/curve25519.pdf
+ *
+ * The underlying modular addition, subtraction, and multiplication
+ * operations are completely redesigned for substantially improved
+ * efficiency compared to the TweetNaCl implementation studied in that
+ * paper.
+ *
+ *				       TweetNaCl	    iPXE
+ *				       ---------	    ----
+ *
+ * Storage size of each big integer	     128	      40
+ * (in bytes)
+ *
+ * Stack usage for key exchange		    1144	     360
+ * (in bytes, large objects only)
+ *
+ * Cost of big integer addition		      16	       5
+ * (in number of 64-bit additions)
+ *
+ * Cost of big integer multiplication	     273	      31
+ * (in number of 64-bit multiplications)
+ *
+ * The implementation is constant-time (provided that the underlying
+ * big integer operations are also constant-time).
+ */
+
+#include <stdint.h>
+#include <string.h>
+#include <assert.h>
+#include <ipxe/init.h>
+#include <ipxe/x25519.h>
+
+/** X25519 reduction constant
+ *
+ * The X25519 field prime is p=2^255-19.  This gives us:
+ *
+ *               p = 2^255 - 19
+ *           2^255 = p + 19
+ *           2^255 = 19          (mod p)
+ *       k * 2^255 = k * 19      (mod p)
+ *
+ * We can therefore reduce a value modulo p by taking the high-order
+ * bits of the value from bit 255 and above, multiplying by 19, and
+ * adding this to the low-order 255 bits of the value.
+ *
+ * This would be cumbersome to do in practice since it would require
+ * partitioning the value at a 255-bit boundary (and hence would
+ * require some shifting and masking operations).  However, we can
+ * note that:
+ *
+ *       k * 2^255 = k * 19      (mod p)
+ *   k * 2 * 2^255 = k * 2 * 19  (mod p)
+ *       k * 2^256 = k * 38      (mod p)
+ *
+ * We can therefore simplify the reduction to taking the high order
+ * bits of the value from bit 256 and above, multiplying by 38, and
+ * adding this to the low-order 256 bits of the value.
+ *
+ * Since 256 will inevitably be a multiple of the big integer element
+ * size (typically 32 or 64 bits), this avoids the need to perform any
+ * shifting or masking operations.
+ */
+#define X25519_REDUCE_256 38
+
+/** X25519 multiplication step 1 result
+ *
+ * Step 1 of X25519 multiplication is to compute the product of two
+ * X25519 unsigned 258-bit integers.
+ *
+ * Both multiplication inputs are limited to 258 bits, and so the
+ * product will have at most 516 bits.
+ */
+union x25519_multiply_step1 {
+	/** Raw product
+	 *
+	 * Big integer multiplication produces a result with a number
+	 * of elements equal to the sum of the number of elements in
+	 * each input.
+	 */
+	bigint_t ( X25519_SIZE + X25519_SIZE ) product;
+	/** Partition into low-order and high-order bits
+	 *
+	 * Reduction modulo p requires separating the low-order 256
+	 * bits from the remaining high-order bits.
+	 *
+	 * Since the value will never exceed 516 bits (see above),
+	 * there will be at most 260 high-order bits.
+	 */
+	struct {
+		/** Low-order 256 bits */
+		bigint_t ( bigint_required_size ( ( 256 /* bits */ + 7 ) / 8 ) )
+			low_256bit;
+		/** High-order 260 bits */
+		bigint_t ( bigint_required_size ( ( 260 /* bits */ + 7 ) / 8 ) )
+			high_260bit;
+	} __attribute__ (( packed )) parts;
+};
+
+/** X25519 multiplication step 2 result
+ *
+ * Step 2 of X25519 multiplication is to multiply the high-order 260
+ * bits from step 1 with the 6-bit reduction constant 38, and to add
+ * this to the low-order 256 bits from step 1.
+ *
+ * The multiplication inputs are limited to 260 and 6 bits
+ * respectively, and so the product will have at most 266 bits.  After
+ * adding the low-order 256 bits from step 1, the result will have at
+ * most 267 bits.
+ */
+union x25519_multiply_step2 {
+	/** Raw product
+	 *
+	 * Big integer multiplication produces a result with a number
+	 * of elements equal to the sum of the number of elements in
+	 * each input.
+	 */
+	bigint_t ( bigint_required_size ( ( 260 /* bits */ + 7 ) / 8 ) +
+		   bigint_required_size ( ( 6 /* bits */ + 7 ) / 8 ) ) product;
+	/** Big integer value
+	 *
+	 * The value will never exceed 267 bits (see above), and so
+	 * may be consumed as a normal X25519 big integer.
+	 */
+	x25519_t value;
+	/** Partition into low-order and high-order bits
+	 *
+	 * Reduction modulo p requires separating the low-order 256
+	 * bits from the remaining high-order bits.
+	 *
+	 * Since the value will never exceed 267 bits (see above),
+	 * there will be at most 11 high-order bits.
+	 */
+	struct {
+		/** Low-order 256 bits */
+		bigint_t ( bigint_required_size ( ( 256 /* bits */ + 7 ) / 8 ) )
+			low_256bit;
+		/** High-order 11 bits */
+		bigint_t ( bigint_required_size ( ( 11 /* bits */ + 7 ) / 8 ) )
+			high_11bit;
+	} __attribute__ (( packed )) parts;
+};
+
+/** X25519 multiplication step 3 result
+ *
+ * Step 3 of X25519 multiplication is to multiply the high-order 11
+ * bits from step 2 with the 6-bit reduction constant 38, and to add
+ * this to the low-order 256 bits from step 2.
+ *
+ * The multiplication inputs are limited to 11 and 6 bits
+ * respectively, and so the product will have at most 17 bits.  After
+ * adding the low-order 256 bits from step 2, the result will have at
+ * most 257 bits.
+ */
+union x25519_multiply_step3 {
+	/** Raw product
+	 *
+	 * Big integer multiplication produces a result with a number
+	 * of elements equal to the sum of the number of elements in
+	 * each input.
+	 */
+	bigint_t ( bigint_required_size ( ( 11 /* bits */ + 7 ) / 8 ) +
+		   bigint_required_size ( ( 6 /* bits */ + 7 ) / 8 ) ) product;
+	/** Big integer value
+	 *
+	 * The value will never exceed 267 bits (see above), and so
+	 * may be consumed as a normal X25519 big integer.
+	 */
+	x25519_t value;
+};
+
+/** X25519 multiplication temporary working space
+ *
+ * We overlap the buffers used by each step of the multiplication
+ * calculation to reduce the total stack space required:
+ *
+ * |--------------------------------------------------------|
+ * | <- pad -> | <------------ step 1 result -------------> |
+ * |           | <- low 256 bits -> | <-- high 260 bits --> |
+ * | <------- step 2 result ------> | <-- step 3 result --> |
+ * |--------------------------------------------------------|
+ */
+union x25519_multiply_workspace {
+	/** Step 1 result */
+	struct {
+		/** Padding to avoid collision between steps 1 and 2
+		 *
+		 * The step 2 multiplication consumes the high 260
+		 * bits of step 1, and so the step 2 multiplication
+		 * result must not overlap this portion of the step 1
+		 * result.
+		 */
+		uint8_t pad[ sizeof ( union x25519_multiply_step2 ) -
+			     offsetof ( union x25519_multiply_step1,
+					parts.high_260bit ) ];
+		/** Step 1 result */
+		union x25519_multiply_step1 step1;
+	} __attribute__ (( packed ));
+	/** Steps 2 and 3 results */
+	struct {
+		/** Step 2 result */
+		union x25519_multiply_step2 step2;
+		/** Step 3 result */
+		union x25519_multiply_step3 step3;
+	} __attribute__ (( packed ));
+};
+
+/** An X25519 elliptic curve point in projective coordinates
+ *
+ * A point (x,y) on the Montgomery curve used in X25519 is represented
+ * using projective coordinates (X/Z,Y/Z) so that intermediate
+ * calculations may be performed on both numerator and denominator
+ * separately, with the division step performed only once at the end
+ * of the calculation.
+ *
+ * The group operation calculation is performed using a Montgomery
+ * ladder as:
+ *
+ *   X[2i]   = ( X[i]^2 - Z[i]^2 )^2
+ *   X[2i+1] = ( X[i] * X[i+1] - Z[i] * Z[i+1] )^2
+ *   Z[2i]   = 4 * X[i] * Z[i] * ( X[i]^2 + A * X[i] * Z[i] + Z[i]^2 )
+ *   Z[2i+1] = X[0] * ( X[i] * Z[i+1] - X[i+1] * Z[i] ) ^ 2
+ *
+ * It is therefore not necessary to store (or use) the value of Y.
+ */
+struct x25519_projective {
+	/** X coordinate */
+	union x25519_quad257 X;
+	/** Z coordinate */
+	union x25519_quad257 Z;
+};
+
+/** An X25519 Montgomery ladder step */
+struct x25519_step {
+	/** X[n]/Z[n] */
+	struct x25519_projective x_n;
+	/** X[n+1]/Z[n+1] */
+	struct x25519_projective x_n1;
+};
+
+/** Constant p=2^255-19 (the finite field prime) */
+static const uint8_t x25519_p_raw[] = {
+	0x7f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+	0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+	0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+	0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xed
+};
+
+/** Constant p=2^255-19 (the finite field prime) */
+static x25519_t x25519_p;
+
+/** Constant 2p=2^256-38 */
+static x25519_t x25519_2p;
+
+/** Constant 4p=2^257-76 */
+static x25519_t x25519_4p;
+
+/** Reduction constant (used during multiplication) */
+static const uint8_t x25519_reduce_256_raw[] = { X25519_REDUCE_256 };
+
+/** Reduction constant (used during multiplication) */
+static bigint_t ( bigint_required_size ( sizeof ( x25519_reduce_256_raw ) ) )
+	x25519_reduce_256;
+
+/** Constant 121665 (used in the Montgomery ladder) */
+static const uint8_t x25519_121665_raw[] = { 0x01, 0xdb, 0x41 };
+
+/** Constant 121665 (used in the Montgomery ladder) */
+static union x25519_oct258 x25519_121665;
+
+/**
+ * Initialise constants
+ *
+ */
+static void x25519_init_constants ( void ) {
+
+	/* Construct constant p */
+	bigint_init ( &x25519_p, x25519_p_raw, sizeof ( x25519_p_raw ) );
+
+	/* Construct constant 2p */
+	bigint_copy ( &x25519_p, &x25519_2p );
+	bigint_add ( &x25519_p, &x25519_2p );
+
+	/* Construct constant 4p */
+	bigint_copy ( &x25519_2p, &x25519_4p );
+	bigint_add ( &x25519_2p, &x25519_4p );
+
+	/* Construct reduction constant */
+	bigint_init ( &x25519_reduce_256, x25519_reduce_256_raw,
+		      sizeof ( x25519_reduce_256_raw ) );
+
+	/* Construct constant 121665 */
+	bigint_init ( &x25519_121665.value, x25519_121665_raw,
+		      sizeof ( x25519_121665_raw ) );
+}
+
+/** Initialisation function */
+struct init_fn x25519_init_fn __init_fn ( INIT_NORMAL ) = {
+	.initialise = x25519_init_constants,
+};
+
+/**
+ * Add big integers modulo field prime
+ *
+ * @v augend		Big integer to add
+ * @v addend		Big integer to add
+ * @v result		Big integer to hold result (may overlap augend)
+ */
+static inline __attribute__ (( always_inline )) void
+x25519_add ( const union x25519_quad257 *augend,
+	     const union x25519_quad257 *addend,
+	     union x25519_oct258 *result ) {
+	int copy;
+
+	/* Copy augend if necessary */
+	copy = ( result != &augend->oct258 );
+	build_assert ( __builtin_constant_p ( copy ) );
+	if ( copy ) {
+		build_assert ( result != &addend->oct258 );
+		bigint_copy ( &augend->oct258.value, &result->value );
+	}
+
+	/* Perform addition
+	 *
+	 * Both inputs are in the range [0,4p-1] and the resulting
+	 * sum is therefore in the range [0,8p-2].
+	 *
+	 * This range lies within the range [0,8p-1] and the result is
+	 * therefore a valid X25519 unsigned 258-bit integer, as
+	 * required.
+	 */
+	bigint_add ( &addend->value, &result->value );
+}
+
+/**
+ * Subtract big integers modulo field prime
+ *
+ * @v minuend		Big integer from which to subtract
+ * @v subtrahend	Big integer to subtract
+ * @v result		Big integer to hold result (may overlap minuend)
+ */
+static inline __attribute__ (( always_inline )) void
+x25519_subtract ( const union x25519_quad257 *minuend,
+		  const union x25519_quad257 *subtrahend,
+		  union x25519_oct258 *result ) {
+	int copy;
+
+	/* Copy minuend if necessary */
+	copy = ( result != &minuend->oct258 );
+	build_assert ( __builtin_constant_p ( copy ) );
+	if ( copy ) {
+		build_assert ( result != &subtrahend->oct258 );
+		bigint_copy ( &minuend->oct258.value, &result->value );
+	}
+
+	/* Perform subtraction
+	 *
+	 * Both inputs are in the range [0,4p-1] and the resulting
+	 * difference is therefore in the range [1-4p,4p-1].
+	 *
+	 * This range lies partially outside the range [0,8p-1] and
+	 * the result is therefore not yet a valid X25519 unsigned
+	 * 258-bit integer.
+	 */
+	bigint_subtract ( &subtrahend->value, &result->value );
+
+	/* Add constant multiple of field prime p
+	 *
+	 * Add the constant 4p to the result.  This brings the result
+	 * within the range [1,8p-1] (without changing the value
+	 * modulo p).
+	 *
+	 * This range lies within the range [0,8p-1] and the result is
+	 * therefore now a valid X25519 unsigned 258-bit integer, as
+	 * required.
+	 */
+	bigint_add ( &x25519_4p, &result->value );
+}
+
+/**
+ * Multiply big integers modulo field prime
+ *
+ * @v multiplicand	Big integer to be multiplied
+ * @v multiplier	Big integer to be multiplied
+ * @v result		Big integer to hold result (may overlap either input)
+ */
+void x25519_multiply ( const union x25519_oct258 *multiplicand,
+		       const union x25519_oct258 *multiplier,
+		       union x25519_quad257 *result ) {
+	union x25519_multiply_workspace tmp;
+	union x25519_multiply_step1 *step1 = &tmp.step1;
+	union x25519_multiply_step2 *step2 = &tmp.step2;
+	union x25519_multiply_step3 *step3 = &tmp.step3;
+
+	/* Step 1: perform raw multiplication
+	 *
+	 *   step1 = multiplicand * multiplier
+	 *
+	 * Both inputs are 258-bit numbers and the step 1 result is
+	 * therefore 258+258=516 bits.
+	 */
+	static_assert ( sizeof ( step1->product ) >= sizeof ( step1->parts ) );
+	bigint_multiply ( &multiplicand->value, &multiplier->value,
+			  &step1->product );
+
+	/* Step 2: reduce high-order 516-256=260 bits of step 1 result
+	 *
+	 * Use the identity 2^256=38 (mod p) to reduce the high-order
+	 * bits of the step 1 result.  We split the 516-bit result
+	 * from step 1 into its low-order 256 bits and high-order 260
+	 * bits:
+	 *
+	 *   step1 = step1(low 256 bits) + step1(high 260 bits) * 2^256
+	 *
+	 * and then perform the calculation:
+	 *
+	 *   step2 = step1                                              (mod p)
+	 *         = step1(low 256 bits) + step1(high 260 bits) * 2^256 (mod p)
+	 *         = step1(low 256 bits) + step1(high 260 bits) * 38    (mod p)
+	 *
+	 * There are 6 bits in the constant value 38.  The step 2
+	 * multiplication product will therefore have 260+6=266 bits,
+	 * and the step 2 result (after the addition) will therefore
+	 * have 267 bits.
+	 */
+	static_assert ( sizeof ( step2->product ) >= sizeof ( step2->value ) );
+	static_assert ( sizeof ( step2->product ) >= sizeof ( step2->parts ) );
+	bigint_grow ( &step1->parts.low_256bit, &result->value );
+	bigint_multiply ( &step1->parts.high_260bit, &x25519_reduce_256,
+			  &step2->product );
+	bigint_add ( &result->value, &step2->value );
+
+	/* Step 3: reduce high-order 267-256=11 bits of step 2 result
+	 *
+	 * Use the identity 2^256=38 (mod p) again to reduce the
+	 * high-order bits of the step 2 result.  As before, we split
+	 * the 267-bit result from step 2 into its low-order 256 bits
+	 * and high-order 11 bits:
+	 *
+	 *   step2 = step2(low 256 bits) + step2(high 11 bits) * 2^256
+	 *
+	 * and then perform the calculation:
+	 *
+	 *   step3 = step2                                             (mod p)
+	 *         = step2(low 256 bits) + step2(high 11 bits) * 2^256 (mod p)
+	 *         = step2(low 256 bits) + step2(high 11 bits) * 38    (mod p)
+	 *
+	 * There are 6 bits in the constant value 38.  The step 3
+	 * multiplication product will therefore have 11+6=19 bits,
+	 * and the step 3 result (after the addition) will therefore
+	 * have 257 bits.
+	 *
+	 * A loose upper bound for the step 3 result (after the
+	 * addition) is given by:
+	 *
+	 *   step3 < ( 2^256 - 1 ) + ( 2^19 - 1 )
+	 *         < ( 2^257 - 2^256 - 1 ) + ( 2^19 - 1 )
+	 *         < ( 2^257 - 76 ) - 2^256 + 2^19 + 74
+	 *         < 4 * ( 2^255 - 19 ) - 2^256 + 2^19 + 74
+	 *         < 4p - 2^256 + 2^19 + 74
+	 *
+	 * and so the step 3 result is strictly less than 4p, and
+	 * therefore lies within the range [0,4p-1].
+	 */
+	memset ( &step3->value, 0, sizeof ( step3->value ) );
+	bigint_grow ( &step2->parts.low_256bit, &result->value );
+	bigint_multiply ( &step2->parts.high_11bit, &x25519_reduce_256,
+			  &step3->product );
+	bigint_add ( &step3->value, &result->value );
+
+	/* Step 1 calculates the product of the input operands, and
+	 * each subsequent step reduces the number of bits in the
+	 * result while preserving this value (modulo p).  The final
+	 * result is therefore equal to the product of the input
+	 * operands (modulo p), as required.
+	 *
+	 * The step 3 result lies within the range [0,4p-1] and the
+	 * final result is therefore a valid X25519 unsigned 257-bit
+	 * integer, as required.
+	 */
+}
+
+/**
+ * Compute multiplicative inverse
+ *
+ * @v invertend		Big integer to be inverted
+ * @v result		Big integer to hold result (may not overlap input)
+ */
+void x25519_invert ( const union x25519_oct258 *invertend,
+		     union x25519_quad257 *result ) {
+	int i;
+
+	/* Sanity check */
+	assert ( invertend != &result->oct258 );
+
+	/* Calculate inverse as x^(-1)=x^(p-2) where p is the field prime
+	 *
+	 * The field prime is p=2^255-19 and so:
+	 *
+	 *   p - 2 = 2^255 - 21
+	 *         = (2^255 - 1) - 2^4 - 2^2
+	 *
+	 * i.e. p-2 is a 254-bit number in which all bits are set
+	 * apart from bit 2 and bit 4.
+	 *
+	 * We use the square-and-multiply method to compute x^(p-2).
+	 */
+	bigint_copy ( &invertend->value, &result->value );
+	for ( i = 253 ; i >= 0 ; i-- ) {
+
+		/* Square running total */
+		x25519_multiply ( &result->oct258, &result->oct258, result );
+
+		/* For each set bit in the exponent, multiply by invertend */
+		if ( ( i != 2 ) && ( i != 4 ) ) {
+			x25519_multiply ( invertend, &result->oct258, result );
+		}
+	}
+}
+
+/**
+ * Reduce big integer via conditional subtraction
+ *
+ * @v subtrahend	Big integer to subtract
+ * @v value		Big integer to be subtracted from, if possible
+ */
+static void x25519_reduce_by ( const x25519_t *subtrahend, x25519_t *value ) {
+	unsigned int max_bit = ( ( 8 * sizeof ( *value ) ) - 1 );
+	x25519_t tmp;
+
+	/* Conditionally subtract subtrahend
+	 *
+	 * Subtract the subtrahend, discarding the result (in constant
+	 * time) if the subtraction underflows.
+	 */
+	bigint_copy ( value, &tmp );
+	bigint_subtract ( subtrahend, value );
+	bigint_swap ( value, &tmp, bigint_bit_is_set ( value, max_bit ) );
+}
+
+/**
+ * Reduce big integer to canonical range
+ *
+ * @v value		Big integer to be reduced
+ */
+void x25519_reduce ( union x25519_quad257 *value ) {
+
+	/* Conditionally subtract 2p
+	 *
+	 * Subtract twice the field prime, discarding the result (in
+	 * constant time) if the subtraction underflows.
+	 *
+	 * The input value is in the range [0,4p-1].  After this
+	 * conditional subtraction, the value is in the range
+	 * [0,2p-1].
+	 */
+	x25519_reduce_by ( &x25519_2p, &value->value );
+
+	/* Conditionally subtract p
+	 *
+	 * Subtract the field prime, discarding the result (in
+	 * constant time) if the subtraction underflows.
+	 *
+	 * The value is already in the range [0,2p-1].  After this
+	 * conditional subtraction, the value is in the range [0,p-1]
+	 * and is therefore the canonical representation.
+	 */
+	x25519_reduce_by ( &x25519_p, &value->value );
+}
+
+/**
+ * Compute next step of the Montgomery ladder
+ *
+ * @v base		Base point
+ * @v bit		Bit value
+ * @v step		Ladder step
+ */
+static void x25519_step ( const union x25519_quad257 *base, int bit,
+			  struct x25519_step *step ) {
+	union x25519_quad257 *a = &step->x_n.X;
+	union x25519_quad257 *b = &step->x_n1.X;
+	union x25519_quad257 *c = &step->x_n.Z;
+	union x25519_quad257 *d = &step->x_n1.Z;
+	union x25519_oct258 e;
+	union x25519_quad257 f;
+	union x25519_oct258 *v1_e;
+	union x25519_oct258 *v2_a;
+	union x25519_oct258 *v3_c;
+	union x25519_oct258 *v4_b;
+	union x25519_quad257 *v5_d;
+	union x25519_quad257 *v6_f;
+	union x25519_quad257 *v7_a;
+	union x25519_quad257 *v8_c;
+	union x25519_oct258 *v9_e;
+	union x25519_oct258 *v10_a;
+	union x25519_quad257 *v11_b;
+	union x25519_oct258 *v12_c;
+	union x25519_quad257 *v13_a;
+	union x25519_oct258 *v14_a;
+	union x25519_quad257 *v15_c;
+	union x25519_quad257 *v16_a;
+	union x25519_quad257 *v17_d;
+	union x25519_quad257 *v18_b;
+
+	/* See the referenced paper "Implementing Curve25519/X25519: A
+	 * Tutorial on Elliptic Curve Cryptography" for the reasoning
+	 * behind this calculation.
+	 */
+
+	/* Reuse storage locations for intermediate results where possible */
+	v1_e = &e;
+	v2_a = container_of ( &a->value, union x25519_oct258, value );
+	v3_c = container_of ( &c->value, union x25519_oct258, value );
+	v4_b = container_of ( &b->value, union x25519_oct258, value );
+	v5_d = d;
+	v6_f = &f;
+	v7_a = a;
+	v8_c = c;
+	v9_e = &e;
+	v10_a = container_of ( &a->value, union x25519_oct258, value );
+	v11_b = b;
+	v12_c = container_of ( &c->value, union x25519_oct258, value );
+	v13_a = a;
+	v14_a = container_of ( &a->value, union x25519_oct258, value );
+	v15_c = c;
+	v16_a = a;
+	v17_d = d;
+	v18_b = b;
+
+	/* Select inputs */
+	bigint_swap ( &a->value, &b->value, bit );
+	bigint_swap ( &c->value, &d->value, bit );
+
+	/* v1 = a + c */
+	x25519_add ( a, c, v1_e );
+
+	/* v2 = a - c */
+	x25519_subtract ( a, c, v2_a );
+
+	/* v3 = b + d */
+	x25519_add ( b, d, v3_c );
+
+	/* v4 = b - d */
+	x25519_subtract ( b, d, v4_b );
+
+	/* v5 = v1^2 = (a + c)^2 = a^2 + 2ac + c^2 */
+	x25519_multiply ( v1_e, v1_e, v5_d );
+
+	/* v6 = v2^2 = (a - c)^2 = a^2 - 2ac + c^2 */
+	x25519_multiply ( v2_a, v2_a, v6_f );
+
+	/* v7 = v3 * v2 = (b + d) * (a - c) = ab - bc + ad - cd */
+	x25519_multiply ( v3_c, v2_a, v7_a );
+
+	/* v8 = v4 * v1 = (b - d) * (a + c) = ab + bc - ad - cd */
+	x25519_multiply ( v4_b, v1_e, v8_c );
+
+	/* v9 = v7 + v8 = 2 * (ab - cd) */
+	x25519_add ( v7_a, v8_c, v9_e );
+
+	/* v10 = v7 - v8 = 2 * (ad - bc) */
+	x25519_subtract ( v7_a, v8_c, v10_a );
+
+	/* v11 = v10^2 = 4 * (ad - bc)^2 */
+	x25519_multiply ( v10_a, v10_a, v11_b );
+
+	/* v12 = v5 - v6 = (a + c)^2 - (a - c)^2 = 4ac */
+	x25519_subtract ( v5_d, v6_f, v12_c );
+
+	/* v13 = v12 * 121665 = 486660ac = (A-2) * ac */
+	x25519_multiply ( v12_c, &x25519_121665, v13_a );
+
+	/* v14 = v13 + v5 = (A-2) * ac + a^2 + 2ac + c^2 = a^2 + A * ac + c^2 */
+	x25519_add ( v13_a, v5_d, v14_a );
+
+	/* v15 = v12 * v14 = 4ac * (a^2 + A * ac + c^2) */
+	x25519_multiply ( v12_c, v14_a, v15_c );
+
+	/* v16 = v5 * v6 = (a + c)^2 * (a - c)^2 = (a^2 - c^2)^2 */
+	x25519_multiply ( &v5_d->oct258, &v6_f->oct258, v16_a );
+
+	/* v17 = v11 * base = 4 * base * (ad - bc)^2 */
+	x25519_multiply ( &v11_b->oct258, &base->oct258, v17_d );
+
+	/* v18 = v9^2 = 4 * (ab - cd)^2 */
+	x25519_multiply ( v9_e, v9_e, v18_b );
+
+	/* Select outputs */
+	bigint_swap ( &a->value, &b->value, bit );
+	bigint_swap ( &c->value, &d->value, bit );
+}
+
+/**
+ * Multiply X25519 elliptic curve point
+ *
+ * @v base		Base point
+ * @v scalar		Scalar multiple
+ * @v result		Point to hold result (may overlap base point)
+ */
+static void x25519_ladder ( const union x25519_quad257 *base,
+			    struct x25519_value *scalar,
+			    union x25519_quad257 *result ) {
+	static const uint8_t zero[] = { 0 };
+	static const uint8_t one[] = { 1 };
+	struct x25519_step step;
+	union x25519_quad257 *tmp;
+	int bit;
+	int i;
+
+	/* Initialise ladder */
+	bigint_init ( &step.x_n.X.value, one, sizeof ( one ) );
+	bigint_init ( &step.x_n.Z.value, zero, sizeof ( zero ) );
+	bigint_copy ( &base->value, &step.x_n1.X.value );
+	bigint_init ( &step.x_n1.Z.value, one, sizeof ( one ) );
+
+	/* Use ladder */
+	for ( i = 254 ; i >= 0 ; i-- ) {
+		bit = ( ( scalar->raw[ i / 8 ] >> ( i % 8 ) ) & 1 );
+		x25519_step ( base, bit, &step );
+	}
+
+	/* Convert back to affine coordinate */
+	tmp = &step.x_n1.X;
+	x25519_invert ( &step.x_n.Z.oct258, tmp );
+	x25519_multiply ( &step.x_n.X.oct258, &tmp->oct258, result );
+	x25519_reduce ( result );
+}
+
+/**
+ * Reverse X25519 value endianness
+ *
+ * @v value		Value to reverse
+ */
+static void x25519_reverse ( struct x25519_value *value ) {
+	uint8_t *low = value->raw;
+	uint8_t *high = &value->raw[ sizeof ( value->raw ) - 1 ];
+	uint8_t tmp;
+
+	/* Reverse bytes */
+	do {
+		tmp = *low;
+		*low = *high;
+		*high = tmp;
+	} while ( ++low < --high );
+}
+
+/**
+ * Calculate X25519 key
+ *
+ * @v base		Base point
+ * @v scalar		Scalar multiple
+ * @v result		Point to hold result (may overlap base point)
+ */
+void x25519_key ( const struct x25519_value *base,
+		  const struct x25519_value *scalar,
+		  struct x25519_value *result ) {
+	struct x25519_value *tmp = result;
+	union x25519_quad257 point;
+
+	/* Reverse base point and clear high bit as required by RFC7748 */
+	memcpy ( tmp, base, sizeof ( *tmp ) );
+	x25519_reverse ( tmp );
+	tmp->raw[0] &= 0x7f;
+	bigint_init ( &point.value, tmp->raw, sizeof ( tmp->raw ) );
+
+	/* Clamp scalar as required by RFC7748 */
+	memcpy ( tmp, scalar, sizeof ( *tmp ) );
+	tmp->raw[0] &= 0xf8;
+	tmp->raw[31] |= 0x40;
+
+	/* Multiply elliptic curve point */
+	x25519_ladder ( &point, tmp, &point );
+
+	/* Reverse result */
+	bigint_done ( &point.value, result->raw, sizeof ( result->raw ) );
+	x25519_reverse ( result );
+}

src/include/ipxe/x25519.h

@@ -0,0 +1,91 @@
+#ifndef _IPXE_X25519_H
+#define _IPXE_X25519_H
+
+/** @file
+ *
+ * X25519 key exchange
+ *
+ */
+
+FILE_LICENCE ( GPL2_OR_LATER_OR_UBDL );
+
+#include <stdint.h>
+#include <ipxe/bigint.h>
+
+/** X25519 unsigned big integer size
+ *
+ * X25519 uses the finite field of integers modulo the prime
+ * p=2^255-19.  The canonical representations of integers in this
+ * field therefore require only 255 bits.
+ *
+ * For internal calculations we use big integers containing up to 267
+ * bits, since this ends up allowing us to avoid some unnecessary (and
+ * expensive) intermediate reductions modulo p.
+ */
+#define X25519_SIZE bigint_required_size ( ( 267 /* bits */ + 7 ) / 8 )
+
+/** An X25519 unsigned big integer used in internal calculations */
+typedef bigint_t ( X25519_SIZE ) x25519_t;
+
+/** An X25519 unsigned 258-bit integer
+ *
+ * This is an unsigned integer N in the finite field of integers
+ * modulo the prime p=2^255-19.
+ *
+ * In this representation, N is encoded as any big integer that is in
+ * the same congruence class as N (i.e that has the same value as N
+ * modulo p) and that lies within the 258-bit range [0,8p-1].
+ *
+ * This type can be used as an input for multiplication (but not for
+ * addition or subtraction).
+ *
+ * Addition or subtraction will produce an output of this type.
+ */
+union x25519_oct258 {
+	/** Big integer value */
+	x25519_t value;
+};
+
+/** An X25519 unsigned 257-bit integer
+ *
+ * This is an unsigned integer N in the finite field of integers
+ * modulo the prime p=2^255-19.
+ *
+ * In this representation, N is encoded as any big integer that is in
+ * the same congruence class as N (i.e that has the same value as N
+ * modulo p) and that lies within the 257-bit range [0,4p-1].
+ *
+ * This type can be used as an input for addition, subtraction, or
+ * multiplication.
+ *
+ * Multiplication will produce an output of this type.
+ */
+union x25519_quad257 {
+	/** Big integer value */
+	x25519_t value;
+	/** X25519 unsigned 258-bit integer
+	 *
+	 * Any value in the range [0,4p-1] is automatically also
+	 * within the range [0,8p-1] and so may be consumed as an
+	 * unsigned 258-bit integer.
+	 */
+	const union x25519_oct258 oct258;
+};
+
+/** An X25519 32-byte value */
+struct x25519_value {
+	/** Raw value */
+	uint8_t raw[32];
+};
+
+extern void x25519_multiply ( const union x25519_oct258 *multiplicand,
+			      const union x25519_oct258 *multiplier,
+			      union x25519_quad257 *result );
+extern void x25519_invert ( const union x25519_oct258 *invertend,
+			    union x25519_quad257 *result );
+extern void x25519_reduce ( union x25519_quad257 *value );
+extern void x25519_key ( const struct x25519_value *base,
+			 const struct x25519_value *scalar,
+			 struct x25519_value *result );
+
+#endif /* _IPXE_X25519_H */

src/tests/tests.c

@@ -81,3 +81,4 @@ REQUIRE_OBJECT ( hmac_test );
 REQUIRE_OBJECT ( dhe_test );
 REQUIRE_OBJECT ( gcm_test );
 REQUIRE_OBJECT ( nap_test );
+REQUIRE_OBJECT ( x25519_test );

src/tests/x25519_test.c

@@ -0,0 +1,570 @@
+/*
+ * 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
+ *
+ * X25519 key exchange test
+ *
+ * Full key exchange test vectors are taken from RFC 7748.
+ */
+
+/* Forcibly enable assertions */
+#undef NDEBUG
+
+#include <stdint.h>
+#include <string.h>
+#include <ipxe/x25519.h>
+#include <ipxe/test.h>
+
+/** Define inline multiplicand */
+#define MULTIPLICAND(...) { __VA_ARGS__ }
+
+/** Define inline multiplier */
+#define MULTIPLIER(...) { __VA_ARGS__ }
+
+/** Define inline invertend */
+#define INVERTEND(...) { __VA_ARGS__ }
+
+/** Define inline base point */
+#define BASE(...) { __VA_ARGS__ }
+
+/** Define inline scalar multiple */
+#define SCALAR(...) { __VA_ARGS__ }
+
+/** Define inline expected result */
+#define EXPECTED(...) { __VA_ARGS__ }
+
+/** An X25519 multiplication self-test */
+struct x25519_multiply_test {
+	/** Multiplicand */
+	const void *multiplicand;
+	/** Length of multiplicand */
+	size_t multiplicand_len;
+	/** Multiplier */
+	const void *multiplier;
+	/** Length of multiplier */
+	size_t multiplier_len;
+	/** Expected result */
+	const void *expected;
+	/** Length of expected result */
+	size_t expected_len;
+};
+
+/**
+ * Define an X25519 multiplication test
+ *
+ * @v name		Test name
+ * @v MULTIPLICAND	258-bit multiplicand
+ * @v MULTIPLIER	258-bit multiplier
+ * @v EXPECTED		255-bit expected result
+ * @ret test		X25519 multiplication test
+ */
+#define X25519_MULTIPLY_TEST( name, MULTIPLICAND, MULTIPLIER,		\
+			      EXPECTED )				\
+	static const uint8_t name ## _multiplicand[] = MULTIPLICAND;	\
+	static const uint8_t name ## _multiplier[] = MULTIPLIER;	\
+	static const uint8_t name ## _expected[] = EXPECTED;		\
+	static struct x25519_multiply_test name = {			\
+		.multiplicand = name ## _multiplicand,			\
+		.multiplicand_len = sizeof ( name ## _multiplicand ),	\
+		.multiplier = name ## _multiplier,			\
+		.multiplier_len = sizeof ( name ## _multiplier ),	\
+		.expected = name ## _expected,				\
+		.expected_len = sizeof ( name ## _expected ),		\
+	}
+
+/** An X25519 multiplicative inversion self-test */
+struct x25519_invert_test {
+	/** Invertend */
+	const void *invertend;
+	/** Length of invertend */
+	size_t invertend_len;
+	/** Expected result */
+	const void *expected;
+	/** Length of expected result */
+	size_t expected_len;
+};
+
+/**
+ * Define an X25519 multiplicative inversion test
+ *
+ * @v name		Test name
+ * @v INVERTEND		258-bit invertend
+ * @v EXPECTED		255-bit expected result
+ * @ret test		X25519 multiplicative inversion test
+ */
+#define X25519_INVERT_TEST( name, INVERTEND, EXPECTED )			\
+	static const uint8_t name ## _invertend[] = INVERTEND;		\
+	static const uint8_t name ## _expected[] = EXPECTED;		\
+	static struct x25519_invert_test name = {			\
+		.invertend = name ## _invertend,			\
+		.invertend_len = sizeof ( name ## _invertend ),		\
+		.expected = name ## _expected,				\
+		.expected_len = sizeof ( name ## _expected ),		\
+	}
+
+/** An X25519 key exchange self-test */
+struct x25519_key_test {
+	/** Base */
+	struct x25519_value base;
+	/** Scalar */
+	struct x25519_value scalar;
+	/** Expected result */
+	struct x25519_value expected;
+	/** Number of iterations */
+	unsigned int count;
+};
+
+/**
+ * Define an X25519 key exchange test
+ *
+ * @v name		Test name
+ * @v COUNT		Number of iterations
+ * @v BASE		Base point
+ * @v SCALAR		Scalar multiple
+ * @v EXPECTED		Expected result
+ * @ret test		X25519 key exchange test
+ */
+#define X25519_KEY_TEST( name, COUNT, BASE, SCALAR, EXPECTED )		\
+	static struct x25519_key_test name = {				\
+		.count = COUNT,						\
+		.base = { .raw = BASE },				\
+		.scalar = { .raw = SCALAR },				\
+		.expected = { .raw = EXPECTED },			\
+	}
+
+/**
+ * Report an X25519 multiplication test result
+ *
+ * @v test		X25519 multiplication test
+ * @v file		Test code file
+ * @v line		Test code line
+ */
+static void x25519_multiply_okx ( struct x25519_multiply_test *test,
+				  const char *file, unsigned int line ) {
+	union x25519_oct258 multiplicand;
+	union x25519_oct258 multiplier;
+	union x25519_quad257 expected;
+	union x25519_quad257 actual;
+
+	/* Construct big integers */
+	bigint_init ( &multiplicand.value, test->multiplicand,
+		      test->multiplicand_len );
+	DBGC ( test, "X25519 multiplicand:\n" );
+	DBGC_HDA ( test, 0, &multiplicand, sizeof ( multiplicand ) );
+	bigint_init ( &multiplier.value, test->multiplier,
+		      test->multiplier_len );
+	DBGC ( test, "X25519 multiplier:\n" );
+	DBGC_HDA ( test, 0, &multiplier, sizeof ( multiplier ) );
+	bigint_init ( &expected.value, test->expected, test->expected_len );
+	DBGC ( test, "X25519 expected product:\n" );
+	DBGC_HDA ( test, 0, &expected, sizeof ( expected ) );
+
+	/* Perform multiplication */
+	x25519_multiply ( &multiplicand, &multiplier, &actual );
+
+	/* Reduce result to allow for comparison */
+	x25519_reduce ( &actual );
+	DBGC ( test, "X25519 actual product:\n" );
+	DBGC_HDA ( test, 0, &actual, sizeof ( actual ) );
+
+	/* Compare against expected result */
+	okx ( memcmp ( &actual, &expected, sizeof ( expected ) ) == 0,
+	      file, line );
+}
+#define x25519_multiply_ok( test ) \
+	x25519_multiply_okx ( test, __FILE__, __LINE__ )
+
+/**
+ * Report an X25519 multiplicative inversion test result
+ *
+ * @v test		X25519 multiplicative inversion test
+ * @v file		Test code file
+ * @v line		Test code line
+ */
+static void x25519_invert_okx ( struct x25519_invert_test *test,
+				const char *file, unsigned int line ) {
+	static const uint8_t one[] = { 1 };
+	union x25519_oct258 invertend;
+	union x25519_quad257 expected;
+	union x25519_quad257 actual;
+	union x25519_quad257 product;
+	union x25519_quad257 identity;
+
+	/* Construct big integers */
+	bigint_init ( &invertend.value, test->invertend, test->invertend_len );
+	DBGC ( test, "X25519 invertend:\n" );
+	DBGC_HDA ( test, 0, &invertend, sizeof ( invertend ) );
+	bigint_init ( &expected.value, test->expected, test->expected_len );
+	DBGC ( test, "X25519 expected inverse:\n" );
+	DBGC_HDA ( test, 0, &expected, sizeof ( expected ) );
+	bigint_init ( &identity.value, one, sizeof ( one ) );
+
+	/* Perform inversion */
+	x25519_invert ( &invertend, &actual );
+
+	/* Multiply invertend by inverse */
+	x25519_multiply ( &invertend, &actual.oct258, &product );
+
+	/* Reduce results to allow for comparison */
+	x25519_reduce ( &actual );
+	DBGC ( test, "X25519 actual inverse:\n" );
+	DBGC_HDA ( test, 0, &actual, sizeof ( actual ) );
+	x25519_reduce ( &product );
+	DBGC ( test, "X25519 actual product:\n" );
+	DBGC_HDA ( test, 0, &product, sizeof ( product ) );
+
+	/* Compare against expected results */
+	okx ( memcmp ( &actual, &expected, sizeof ( expected ) ) == 0,
+	      file, line );
+	okx ( memcmp ( &product, &identity, sizeof ( identity ) ) == 0,
+	      file, line );
+}
+#define x25519_invert_ok( test ) \
+	x25519_invert_okx ( test, __FILE__, __LINE__ )
+
+/**
+ * Report an X25519 key exchange test result
+ *
+ * @v test		X25519 key exchange test
+ * @v file		Test code file
+ * @v line		Test code line
+ */
+static void x25519_key_okx ( struct x25519_key_test *test,
+			     const char *file, unsigned int line ) {
+	struct x25519_value base;
+	struct x25519_value scalar;
+	struct x25519_value actual;
+	unsigned int i;
+
+	/* Construct input values */
+	memcpy ( &base, &test->base, sizeof ( test->base ) );
+	memcpy ( &scalar, &test->scalar, sizeof ( test->scalar ) );
+	DBGC ( test, "X25519 base:\n" );
+	DBGC_HDA ( test, 0, &base, sizeof ( base ) );
+	DBGC ( test, "X25519 scalar:\n" );
+	DBGC_HDA ( test, 0, &scalar, sizeof ( scalar ) );
+	DBGC ( test, "X25519 expected result (x%d):\n", test->count );
+	DBGC_HDA ( test, 0, &test->expected, sizeof ( test->expected ) );
+
+	/* Calculate key */
+	for ( i = 0 ; i < test->count ; i++ ) {
+		x25519_key ( &base, &scalar, &actual );
+		memcpy ( &base, &scalar, sizeof ( base ) );
+		memcpy ( &scalar, &actual, sizeof ( scalar ) );
+	}
+	DBGC ( test, "X25519 actual result (x%d):\n", test->count );
+	DBGC_HDA ( test, 0, &actual, sizeof ( actual ) );
+
+	/* Compare against expected result */
+	okx ( memcmp ( &actual, &test->expected,
+		       sizeof ( test->expected ) ) == 0, file, line );
+}
+#define x25519_key_ok( test ) \
+	x25519_key_okx ( test, __FILE__, __LINE__ )
+
+/* Test multiplying small numbers */
+X25519_MULTIPLY_TEST ( multiply_small, MULTIPLICAND ( 6 ),
+		       MULTIPLIER ( 9 ), EXPECTED ( 6 * 9 ) );
+
+/* Test exact multiple of field prime */
+X25519_MULTIPLY_TEST ( multiply_k_p,
+	MULTIPLICAND ( 0x7f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+		       0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+		       0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+		       0xff, 0xff, 0xff, 0xff, 0xed ),
+	MULTIPLIER ( 0x00, 0xe8, 0x0d, 0x83, 0xd4, 0xe9, 0x1e, 0xdd, 0x7a,
+		     0x45, 0x14, 0x87, 0xb7, 0xfc, 0x62, 0x54, 0x1f, 0xb2,
+		     0x97, 0x24, 0xde, 0xfa, 0xd3, 0xe7, 0x3e, 0x83, 0x93,
+		     0x60, 0xbc, 0x20, 0x97, 0x9b, 0x22 ),
+	EXPECTED ( 0x00 ) );
+
+/* 0x0223b8c1e9392456de3eb13b9046685257bdd640fb06671ad11c80317fa3b1799d *
+ * 0x006c031199972a846916419f828b9d2434e465e150bd9c66b3ad3c2d6d1a3d1fa7 =
+ * 0x1ba87e982f7c477616b4d5136ba54733e40081c1c2e27d864aa178ce893d1297 (mod p)
+ */
+X25519_MULTIPLY_TEST ( multiply_1,
+	MULTIPLICAND ( 0x02, 0x23, 0xb8, 0xc1, 0xe9, 0x39, 0x24, 0x56, 0xde,
+		       0x3e, 0xb1, 0x3b, 0x90, 0x46, 0x68, 0x52, 0x57, 0xbd,
+		       0xd6, 0x40, 0xfb, 0x06, 0x67, 0x1a, 0xd1, 0x1c, 0x80,
+		       0x31, 0x7f, 0xa3, 0xb1, 0x79, 0x9d ),
+	MULTIPLIER ( 0x00, 0x6c, 0x03, 0x11, 0x99, 0x97, 0x2a, 0x84, 0x69,
+		     0x16, 0x41, 0x9f, 0x82, 0x8b, 0x9d, 0x24, 0x34, 0xe4,
+		     0x65, 0xe1, 0x50, 0xbd, 0x9c, 0x66, 0xb3, 0xad, 0x3c,
+		     0x2d, 0x6d, 0x1a, 0x3d, 0x1f, 0xa7 ),
+	EXPECTED ( 0x1b, 0xa8, 0x7e, 0x98, 0x2f, 0x7c, 0x47, 0x76, 0x16, 0xb4,
+		   0xd5, 0x13, 0x6b, 0xa5, 0x47, 0x33, 0xe4, 0x00, 0x81, 0xc1,
+		   0xc2, 0xe2, 0x7d, 0x86, 0x4a, 0xa1, 0x78, 0xce, 0x89, 0x3d,
+		   0x12, 0x97 ) );
+
+/* 0x008fadc1a606cb0fb39a1de644815ef6d13b8faa1837f8a88b17fc695a07a0ca6e *
+ * 0x0196da1dac72ff5d2a386ecbe06b65a6a48b8148f6b38a088ca65ed389b74d0fb1 =
+ * 0x351f7bf75ef580249ed6f9ff3996463b0730a1d49b5d36b863e192591157e950 (mod p)
+ */
+X25519_MULTIPLY_TEST ( multiply_2,
+	MULTIPLICAND ( 0x00, 0x8f, 0xad, 0xc1, 0xa6, 0x06, 0xcb, 0x0f, 0xb3,
+		       0x9a, 0x1d, 0xe6, 0x44, 0x81, 0x5e, 0xf6, 0xd1, 0x3b,
+		       0x8f, 0xaa, 0x18, 0x37, 0xf8, 0xa8, 0x8b, 0x17, 0xfc,
+		       0x69, 0x5a, 0x07, 0xa0, 0xca, 0x6e ),
+	MULTIPLIER ( 0x01, 0x96, 0xda, 0x1d, 0xac, 0x72, 0xff, 0x5d, 0x2a,
+		     0x38, 0x6e, 0xcb, 0xe0, 0x6b, 0x65, 0xa6, 0xa4, 0x8b,
+		     0x81, 0x48, 0xf6, 0xb3, 0x8a, 0x08, 0x8c, 0xa6, 0x5e,
+		     0xd3, 0x89, 0xb7, 0x4d, 0x0f, 0xb1 ),
+	EXPECTED ( 0x35, 0x1f, 0x7b, 0xf7, 0x5e, 0xf5, 0x80, 0x24, 0x9e, 0xd6,
+		   0xf9, 0xff, 0x39, 0x96, 0x46, 0x3b, 0x07, 0x30, 0xa1, 0xd4,
+		   0x9b, 0x5d, 0x36, 0xb8, 0x63, 0xe1, 0x92, 0x59, 0x11, 0x57,
+		   0xe9, 0x50 ) );
+
+/* 0x016c307511b2b9437a28df6ec4ce4a2bbdc241330b01a9e71fde8a774bcf36d58b *
+ * 0x0117be31111a2a73ed562b0f79c37459eef50bea63371ecd7b27cd813047229389 =
+ * 0x6b43b5185965f8f0920f31ae1b2cefadd7b078fecf68dbeaa17b9c385b558329 (mod p)
+ */
+X25519_MULTIPLY_TEST ( multiply_3,
+	MULTIPLICAND ( 0x01, 0x6c, 0x30, 0x75, 0x11, 0xb2, 0xb9, 0x43, 0x7a,
+		       0x28, 0xdf, 0x6e, 0xc4, 0xce, 0x4a, 0x2b, 0xbd, 0xc2,
+		       0x41, 0x33, 0x0b, 0x01, 0xa9, 0xe7, 0x1f, 0xde, 0x8a,
+		       0x77, 0x4b, 0xcf, 0x36, 0xd5, 0x8b ),
+	MULTIPLIER ( 0x01, 0x17, 0xbe, 0x31, 0x11, 0x1a, 0x2a, 0x73, 0xed,
+		     0x56, 0x2b, 0x0f, 0x79, 0xc3, 0x74, 0x59, 0xee, 0xf5,
+		     0x0b, 0xea, 0x63, 0x37, 0x1e, 0xcd, 0x7b, 0x27, 0xcd,
+		     0x81, 0x30, 0x47, 0x22, 0x93, 0x89 ),
+	EXPECTED ( 0x6b, 0x43, 0xb5, 0x18, 0x59, 0x65, 0xf8, 0xf0, 0x92, 0x0f,
+		   0x31, 0xae, 0x1b, 0x2c, 0xef, 0xad, 0xd7, 0xb0, 0x78, 0xfe,
+		   0xcf, 0x68, 0xdb, 0xea, 0xa1, 0x7b, 0x9c, 0x38, 0x5b, 0x55,
+		   0x83, 0x29 ) );
+
+/* 0x020b1f9163ce9ff57f43b7a3a69a8dca03580d7b71d8f564135be6128e18c26797 *
+ * 0x018d5288f1142c3fe860e7a113ec1b8ca1f91e1d4c1ff49b7889463e85759cde66 =
+ * 0x28a77d3c8a14323d63b288dbd40315b3f192b8485d86a02cb87d3dfb7a0b5447 (mod p)
+ */
+X25519_MULTIPLY_TEST ( multiply_4,
+	MULTIPLICAND ( 0x02, 0x0b, 0x1f, 0x91, 0x63, 0xce, 0x9f, 0xf5, 0x7f,
+		       0x43, 0xb7, 0xa3, 0xa6, 0x9a, 0x8d, 0xca, 0x03, 0x58,
+		       0x0d, 0x7b, 0x71, 0xd8, 0xf5, 0x64, 0x13, 0x5b, 0xe6,
+		       0x12, 0x8e, 0x18, 0xc2, 0x67, 0x97 ),
+	MULTIPLIER ( 0x01, 0x8d, 0x52, 0x88, 0xf1, 0x14, 0x2c, 0x3f, 0xe8,
+		     0x60, 0xe7, 0xa1, 0x13, 0xec, 0x1b, 0x8c, 0xa1, 0xf9,
+		     0x1e, 0x1d, 0x4c, 0x1f, 0xf4, 0x9b, 0x78, 0x89, 0x46,
+		     0x3e, 0x85, 0x75, 0x9c, 0xde, 0x66 ),
+	EXPECTED ( 0x28, 0xa7, 0x7d, 0x3c, 0x8a, 0x14, 0x32, 0x3d, 0x63, 0xb2,
+		   0x88, 0xdb, 0xd4, 0x03, 0x15, 0xb3, 0xf1, 0x92, 0xb8, 0x48,
+		   0x5d, 0x86, 0xa0, 0x2c, 0xb8, 0x7d, 0x3d, 0xfb, 0x7a, 0x0b,
+		   0x54, 0x47 ) );
+
+/* 0x023139d32c93cd59bf5c941cf0dc98d2c1e2acf72f9e574f7aa0ee89aed453dd32 *
+ * 0x03146d3f31fc377a4c4a15544dc5e7ce8a3a578a8ea9488d990bbb259911ce5dd2 =
+ * 0x4bdb7a35c0a5182000aa67554741e88cfdf460a78c6fae07adf83d2f005d2767 (mod p)
+ */
+X25519_MULTIPLY_TEST ( multiply_5,
+	MULTIPLICAND ( 0x02, 0x31, 0x39, 0xd3, 0x2c, 0x93, 0xcd, 0x59, 0xbf,
+		       0x5c, 0x94, 0x1c, 0xf0, 0xdc, 0x98, 0xd2, 0xc1, 0xe2,
+		       0xac, 0xf7, 0x2f, 0x9e, 0x57, 0x4f, 0x7a, 0xa0, 0xee,
+		       0x89, 0xae, 0xd4, 0x53, 0xdd, 0x32 ),
+	MULTIPLIER ( 0x03, 0x14, 0x6d, 0x3f, 0x31, 0xfc, 0x37, 0x7a, 0x4c,
+		     0x4a, 0x15, 0x54, 0x4d, 0xc5, 0xe7, 0xce, 0x8a, 0x3a,
+		     0x57, 0x8a, 0x8e, 0xa9, 0x48, 0x8d, 0x99, 0x0b, 0xbb,
+		     0x25, 0x99, 0x11, 0xce, 0x5d, 0xd2 ),
+	EXPECTED ( 0x4b, 0xdb, 0x7a, 0x35, 0xc0, 0xa5, 0x18, 0x20, 0x00, 0xaa,
+		   0x67, 0x55, 0x47, 0x41, 0xe8, 0x8c, 0xfd, 0xf4, 0x60, 0xa7,
+		   0x8c, 0x6f, 0xae, 0x07, 0xad, 0xf8, 0x3d, 0x2f, 0x00, 0x5d,
+		   0x27, 0x67 ) );
+
+/* 0x01d58842dea2bc372f7412b29347294739614ff3d719db3ad0ddd1dfb23b982ef8 ^ -1 =
+ * 0x093ff51750809d181a9a5481c564e37cff618def8ec45f464b1a6e24f8b826bd (mod p)
+ */
+X25519_INVERT_TEST ( invert_1,
+	INVERTEND ( 0x01, 0xd5, 0x88, 0x42, 0xde, 0xa2, 0xbc, 0x37, 0x2f,
+		    0x74, 0x12, 0xb2, 0x93, 0x47, 0x29, 0x47, 0x39, 0x61,
+		    0x4f, 0xf3, 0xd7, 0x19, 0xdb, 0x3a, 0xd0, 0xdd, 0xd1,
+		    0xdf, 0xb2, 0x3b, 0x98, 0x2e, 0xf8 ),
+	EXPECTED ( 0x09, 0x3f, 0xf5, 0x17, 0x50, 0x80, 0x9d, 0x18, 0x1a, 0x9a,
+		   0x54, 0x81, 0xc5, 0x64, 0xe3, 0x7c, 0xff, 0x61, 0x8d, 0xef,
+		   0x8e, 0xc4, 0x5f, 0x46, 0x4b, 0x1a, 0x6e, 0x24, 0xf8, 0xb8,
+		   0x26, 0xbd ) );
+
+/* 0x02efc89849b3aa7efe4458a885ab9099a435a240ae5af305535ec42e0829a3b2e9 ^ -1 =
+ * 0x591607b163e89d0ac33a62c881e984a25d3826e3db5ce229af240dc58e5b579a (mod p)
+ */
+X25519_INVERT_TEST ( invert_2,
+	INVERTEND ( 0x02, 0xef, 0xc8, 0x98, 0x49, 0xb3, 0xaa, 0x7e, 0xfe,
+		    0x44, 0x58, 0xa8, 0x85, 0xab, 0x90, 0x99, 0xa4, 0x35,
+		    0xa2, 0x40, 0xae, 0x5a, 0xf3, 0x05, 0x53, 0x5e, 0xc4,
+		    0x2e, 0x08, 0x29, 0xa3, 0xb2, 0xe9 ),
+	EXPECTED ( 0x59, 0x16, 0x07, 0xb1, 0x63, 0xe8, 0x9d, 0x0a, 0xc3, 0x3a,
+		   0x62, 0xc8, 0x81, 0xe9, 0x84, 0xa2, 0x5d, 0x38, 0x26, 0xe3,
+		   0xdb, 0x5c, 0xe2, 0x29, 0xaf, 0x24, 0x0d, 0xc5, 0x8e, 0x5b,
+		   0x57, 0x9a ) );
+
+/* 0x003eabedcbbaa80dd488bd64072bcfbe01a28defe39bf0027312476f57a5e5a5ab ^ -1 =
+ * 0x7d87c2e565b27c5038181a0a7cae9ebe826c8afc1f77128a4d62cce96d2759a2 (mod p)
+ */
+X25519_INVERT_TEST ( invert_3,
+	INVERTEND ( 0x00, 0x3e, 0xab, 0xed, 0xcb, 0xba, 0xa8, 0x0d, 0xd4,
+		    0x88, 0xbd, 0x64, 0x07, 0x2b, 0xcf, 0xbe, 0x01, 0xa2,
+		    0x8d, 0xef, 0xe3, 0x9b, 0xf0, 0x02, 0x73, 0x12, 0x47,
+		    0x6f, 0x57, 0xa5, 0xe5, 0xa5, 0xab ),
+	EXPECTED ( 0x7d, 0x87, 0xc2, 0xe5, 0x65, 0xb2, 0x7c, 0x50, 0x38, 0x18,
+		   0x1a, 0x0a, 0x7c, 0xae, 0x9e, 0xbe, 0x82, 0x6c, 0x8a, 0xfc,
+		   0x1f, 0x77, 0x12, 0x8a, 0x4d, 0x62, 0xcc, 0xe9, 0x6d, 0x27,
+		   0x59, 0xa2 ) );
+
+/* 0x008e944239b02b61c4a3d70628ece66fa2fd5166e6451b4cf36123fdf77656af72 ^ -1 =
+ * 0x08e96161a0eee1b29af396f154950d5c715dc61aff66ee97377ab22adf3321d7 (mod p)
+ */
+X25519_INVERT_TEST ( invert_4,
+	INVERTEND ( 0x00, 0x8e, 0x94, 0x42, 0x39, 0xb0, 0x2b, 0x61, 0xc4,
+		    0xa3, 0xd7, 0x06, 0x28, 0xec, 0xe6, 0x6f, 0xa2, 0xfd,
+		    0x51, 0x66, 0xe6, 0x45, 0x1b, 0x4c, 0xf3, 0x61, 0x23,
+		    0xfd, 0xf7, 0x76, 0x56, 0xaf, 0x72 ),
+	EXPECTED ( 0x08, 0xe9, 0x61, 0x61, 0xa0, 0xee, 0xe1, 0xb2, 0x9a, 0xf3,
+		   0x96, 0xf1, 0x54, 0x95, 0x0d, 0x5c, 0x71, 0x5d, 0xc6, 0x1a,
+		   0xff, 0x66, 0xee, 0x97, 0x37, 0x7a, 0xb2, 0x2a, 0xdf, 0x33,
+		   0x21, 0xd7 ) );
+
+/* 0x00d261a7ab3aa2e4f90e51f30dc6a7ee39c4b032ccd7c524a55304317faf42e12f ^ -1 =
+ * 0x0738781c0aeabfbe6e840c85bd30996ef71bc54988ce16cedd5ab4f30c281597 (mod p)
+ */
+X25519_INVERT_TEST ( invert_5,
+	INVERTEND ( 0x00, 0xd2, 0x61, 0xa7, 0xab, 0x3a, 0xa2, 0xe4, 0xf9,
+		    0x0e, 0x51, 0xf3, 0x0d, 0xc6, 0xa7, 0xee, 0x39, 0xc4,
+		    0xb0, 0x32, 0xcc, 0xd7, 0xc5, 0x24, 0xa5, 0x53, 0x04,
+		    0x31, 0x7f, 0xaf, 0x42, 0xe1, 0x2f ),
+	EXPECTED ( 0x07, 0x38, 0x78, 0x1c, 0x0a, 0xea, 0xbf, 0xbe, 0x6e, 0x84,
+		   0x0c, 0x85, 0xbd, 0x30, 0x99, 0x6e, 0xf7, 0x1b, 0xc5, 0x49,
+		   0x88, 0xce, 0x16, 0xce, 0xdd, 0x5a, 0xb4, 0xf3, 0x0c, 0x28,
+		   0x15, 0x97 ) );
+
+/* Base:   0xe6db6867583030db3594c1a424b15f7c726624ec26b3353b10a903a6d0ab1c4c
+ * Scalar: 0xa546e36bf0527c9d3b16154b82465edd62144c0ac1fc5a18506a2244ba449ac4
+ * Result: 0xc3da55379de9c6908e94ea4df28d084f32eccf03491c71f754b4075577a28552
+ */
+X25519_KEY_TEST ( rfc7748_1, 1,
+	BASE ( 0xe6, 0xdb, 0x68, 0x67, 0x58, 0x30, 0x30, 0xdb, 0x35, 0x94,
+	       0xc1, 0xa4, 0x24, 0xb1, 0x5f, 0x7c, 0x72, 0x66, 0x24, 0xec,
+	       0x26, 0xb3, 0x35, 0x3b, 0x10, 0xa9, 0x03, 0xa6, 0xd0, 0xab,
+	       0x1c, 0x4c ),
+	SCALAR ( 0xa5, 0x46, 0xe3, 0x6b, 0xf0, 0x52, 0x7c, 0x9d, 0x3b, 0x16,
+		 0x15, 0x4b, 0x82, 0x46, 0x5e, 0xdd, 0x62, 0x14, 0x4c, 0x0a,
+		 0xc1, 0xfc, 0x5a, 0x18, 0x50, 0x6a, 0x22, 0x44, 0xba, 0x44,
+		 0x9a, 0xc4 ),
+	EXPECTED ( 0xc3, 0xda, 0x55, 0x37, 0x9d, 0xe9, 0xc6, 0x90, 0x8e, 0x94,
+		   0xea, 0x4d, 0xf2, 0x8d, 0x08, 0x4f, 0x32, 0xec, 0xcf, 0x03,
+		   0x49, 0x1c, 0x71, 0xf7, 0x54, 0xb4, 0x07, 0x55, 0x77, 0xa2,
+		   0x85, 0x52 ) );
+
+/* Base:   0xe5210f12786811d3f4b7959d0538ae2c31dbe7106fc03c3efc4cd549c715a493
+ * Scalar: 0x4b66e9d4d1b4673c5ad22691957d6af5c11b6421e0ea01d42ca4169e7918ba0d
+ * Result: 0x95cbde9476e8907d7aade45cb4b873f88b595a68799fa152e6f8f7647aac7957
+ */
+X25519_KEY_TEST ( rfc7748_2, 1,
+	BASE ( 0xe5, 0x21, 0x0f, 0x12, 0x78, 0x68, 0x11, 0xd3, 0xf4, 0xb7,
+	       0x95, 0x9d, 0x05, 0x38, 0xae, 0x2c, 0x31, 0xdb, 0xe7, 0x10,
+	       0x6f, 0xc0, 0x3c, 0x3e, 0xfc, 0x4c, 0xd5, 0x49, 0xc7, 0x15,
+	       0xa4, 0x93 ),
+	SCALAR ( 0x4b, 0x66, 0xe9, 0xd4, 0xd1, 0xb4, 0x67, 0x3c, 0x5a, 0xd2,
+		 0x26, 0x91, 0x95, 0x7d, 0x6a, 0xf5, 0xc1, 0x1b, 0x64, 0x21,
+		 0xe0, 0xea, 0x01, 0xd4, 0x2c, 0xa4, 0x16, 0x9e, 0x79, 0x18,
+		 0xba, 0x0d ),
+	EXPECTED ( 0x95, 0xcb, 0xde, 0x94, 0x76, 0xe8, 0x90, 0x7d, 0x7a, 0xad,
+		   0xe4, 0x5c, 0xb4, 0xb8, 0x73, 0xf8, 0x8b, 0x59, 0x5a, 0x68,
+		   0x79, 0x9f, 0xa1, 0x52, 0xe6, 0xf8, 0xf7, 0x64, 0x7a, 0xac,
+		   0x79, 0x57 ) );
+
+/* Base:   0x0900000000000000000000000000000000000000000000000000000000000000
+ * Scalar: 0x0900000000000000000000000000000000000000000000000000000000000000
+ * Result: 0x422c8e7a6227d7bca1350b3e2bb7279f7897b87bb6854b783c60e80311ae3079
+ */
+X25519_KEY_TEST ( rfc7748_3, 1,
+	BASE ( 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	       0x00, 0x00 ),
+	SCALAR ( 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		 0x00, 0x00 ),
+	EXPECTED ( 0x42, 0x2c, 0x8e, 0x7a, 0x62, 0x27, 0xd7, 0xbc, 0xa1, 0x35,
+		   0x0b, 0x3e, 0x2b, 0xb7, 0x27, 0x9f, 0x78, 0x97, 0xb8, 0x7b,
+		   0xb6, 0x85, 0x4b, 0x78, 0x3c, 0x60, 0xe8, 0x03, 0x11, 0xae,
+		   0x30, 0x79 ) );
+
+/* Base:   0x0900000000000000000000000000000000000000000000000000000000000000
+ * Scalar: 0x0900000000000000000000000000000000000000000000000000000000000000
+ * Result: 0xb1a5a73158904c020866c13939dd7e1aa26852ee1d2609c92e5a8f1debe2150a
+ * (after 100 iterations)
+ *
+ * RFC 7748 gives test vectors for 1000 and 1000000 iterations with
+ * these starting values.  This test case stops after 100 iterations
+ * to avoid a pointlessly slow test cycle in the common case of
+ * running tests under Valgrind.
+ */
+X25519_KEY_TEST ( rfc7748_4_100, 100,
+	BASE ( 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	       0x00, 0x00 ),
+	SCALAR ( 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		 0x00, 0x00 ),
+	EXPECTED ( 0xb1, 0xa5, 0xa7, 0x31, 0x58, 0x90, 0x4c, 0x02, 0x08, 0x66,
+		   0xc1, 0x39, 0x39, 0xdd, 0x7e, 0x1a, 0xa2, 0x68, 0x52, 0xee,
+		   0x1d, 0x26, 0x09, 0xc9, 0x2e, 0x5a, 0x8f, 0x1d, 0xeb, 0xe2,
+		   0x15, 0x0a ) );
+
+/**
+ * Perform X25519 self-tests
+ *
+ */
+static void x25519_test_exec ( void ) {
+
+	/* Perform multiplication tests */
+	x25519_multiply_ok ( &multiply_small );
+	x25519_multiply_ok ( &multiply_k_p );
+	x25519_multiply_ok ( &multiply_1 );
+	x25519_multiply_ok ( &multiply_2 );
+	x25519_multiply_ok ( &multiply_3 );
+	x25519_multiply_ok ( &multiply_4 );
+	x25519_multiply_ok ( &multiply_5 );
+
+	/* Perform multiplicative inversion tests */
+	x25519_invert_ok ( &invert_1 );
+	x25519_invert_ok ( &invert_2 );
+	x25519_invert_ok ( &invert_3 );
+	x25519_invert_ok ( &invert_4 );
+	x25519_invert_ok ( &invert_5 );
+
+	/* Perform key exchange tests */
+	x25519_key_ok ( &rfc7748_1 );
+	x25519_key_ok ( &rfc7748_2 );
+	x25519_key_ok ( &rfc7748_3 );
+	x25519_key_ok ( &rfc7748_4_100 );
+}
+
+/** X25519 self-test */
+struct self_test x25519_test __self_test = {
+	.name = "x25519",
+	.exec = x25519_test_exec,
+};