curve25519-dalek-source/docs/avx2-notes.md

An implementation of group operations on the twisted Edwards form of Curve25519, using AVX2 to implement the 4-way parallel formulas of Hisil, Wong, Carter, and Dawson (HWCD). Their 2008 paper Twisted Edwards Curves Revisited, which introduced the extended coordinates used in other parts of -dalek, also describes 4-way parallel formulas for point addition and doubling:

• a unified addition algorithm taking an effective \(2\mathbf M + 1\mathbf D\);

• a doubling algorithm taking an effective \(1\mathbf M + 1\mathbf S\);

• a dedicated (i.e., for distinct points) addition algorithm taking an effective \(2 \mathbf M \).

Here \(\mathbf M\) and \(\mathbf S\) represent the cost of multiplication and squaring of generic field elements and \(\mathbf D\) represents the cost of multiplication by a curve constant.

These formulas do not seem to have been implemented using SIMD before. A 2015 paper by Hernández and López mentions using AVX2 for the X25519 Montgomery ladder, but neither the paper nor the code are publicly available, and it apparently gives only a slight speedup. The 2008 HWCD paper also describes and analyzes a 2-wide variant of the Montgomery ladder (for comparison with parallel Edwards formulas); this strategy was used in 2015 by Tung Chou's sandy2x implementation, which used a 2-wide field implementation in 128-bit vector registers. Curiously, however, although the sandy2x paper also implements Edwards arithmetic, and cites the HWCD paper, it doesn't mention the parallel formulas from HWCD, suggesting that they have been overlooked for software implementations.

The notes below describe a tweak to the \( 2\mathbf M + 1\mathbf D \) unified addition formulas to give \( 2\mathbf M \) readdition with \(1\mathbf D\) precomputation, and a tweak to the doubling formulas to avoid an extra reduction. These tweaked formulas are the ones used by the avx2 backend of curve25519-dalek.

Parallel formulas in HWCD'08

The doubling formula is presented in the HWCD paper as follows:

Cost	Processor 1	Processor 2	Processor 3	Processor 4
	idle	idle	idle	\( R_1 \gets X_1 + Y_1 \)
\(1\mathbf S\)	\( R_2 \gets X_1^2 \)	\( R_3 \gets Y_1^2 \)	\( R_4 \gets Z_1^2 \)	\( R_5 \gets R_1^2 \)
	\( R_6 \gets R_2 + R_3 \)	\( R_7 \gets R_2 - R_3 \)	\( R_4 \gets 2 R_4 \)	idle
	idle	\( R_1 \gets R_4 + R_7 \)	idle	\( R_2 \gets R_6 - R_5 \)
\(1\mathbf M\)	\( X_3 \gets R_1 R_2 \)	\( Y_3 \gets R_6 R_7 \)	\( T_3 \gets R_2 R_6 \)	\( Z_3 \gets R_1 R_7 \)

and the unified addition algorithm is presented as follows:

Cost	Processor 1	Processor 2	Processor 3	Processor 4
	\( R_1 \gets Y_1 - X_1 \)	\( R_2 \gets Y_2 - X_2 \)	\( R_3 \gets Y_1 + X_1 \)	\( R_4 \gets Y_2 + X_2 \)
\(1\mathbf M\)	\( R_5 \gets R_1 R_2 \)	\( R_6 \gets R_3 R_4 \)	\( R_7 \gets T_1 T_2 \)	\( R_8 \gets Z_1 Z_2 \)
\(1\mathbf D\)	idle	idle	\( R_7 \gets k R_7 \)	\( R_8 \gets 2 R_8 \)
	\( R_1 \gets R_6 - R_5 \)	\( R_2 \gets R_8 - R_7 \)	\( R_3 \gets R_8 + R_7 \)	\( R_4 \gets R_6 + R_5 \)
\(1\mathbf M\)	\( X_3 \gets R_1 R_2 \)	\( Y_3 \gets R_3 R_4 \)	\( T_3 \gets R_1 R_4 \)	\( Z_3 \gets R_2 R_3 \)

Here \( k = 2d \) is a curve constant.

For a software implementation, each processor's operations are too low-latency to parallelize across threads. However, the main cost is in the multiplication and squaring steps, which are uniform, while the divergent steps involve inexpensive additions and subtractions.

This means we can use SIMD to implement the expensive portions in parallel, and handle the instruction divergence on the inexpensive parts using masking.

The remaining obstacle to parallelism is the multiplication by the curve constant \(k = 2d\). In the Curve25519 case, this is

 k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. 

HWCD suggest parallelising this step by breaking \(k\) into four parts as \(k = k_0 + 2^n k_1 + 2^{2n} k_2 + 2^{3n} k_3 \) and computing \(k_i R_7 \) in parallel. However, this would be somewhat awkward in our case, since we would normally represent \(k\) as \( 10 \) 32-bit limbs, and \(10 \) is not divisible by \(4\), so we would need a specialized routine to perform a vectorized multiplication by 64-bit constants.

Instead, since we are working projectively, we can multiply \(R_7\) by \( -2\cdot 121665 \) and multiply the other three variables by \(121666\). This trick was suggested by Mike Hamburg. Ignoring the sign for the moment, since \(2 \cdot 121666 < 2^{18}\), all these constants fit in 32 bits, so (up to sign) this can be done in parallel as four multiplications by small constants \( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \).

How do we handle the sign? Since we're primarily interested in Ristretto performance, not Curve25519 performance, we could alternately work on the \(4\)-isogenous "IsoEd25519" curve, which has \(d = 121665\). However, this would only save the negation step, since multiplying one field element by a 32-bit constant is not much easier than multiplying four field elements by 32-bit constants, and it would prevent accelerating Curve25519, so we don't make this choice. Instead, we just negate one lane, and move the \(1 \mathbf D\) into precomputation (see below).

Tweaked formulas

After tweaking the formulas as described above, we obtain the following. To avoid confusion with the original HWCD formulas, temporary variables are named \(S\) instead of \(R\) and are in static single-assignment form.

Addition

This implementation only implements readdition, but the tweaked addition formulas are described first. To add points \(P_1 = (X_1 : Y_1 : Z_1 : T_1) \) and \(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \), we compute

\begin{aligned} S_0 &\gets Y_1 - X_1 \\ S_1 &\gets Y_1 + X_1 \\ S_2 &\gets Y_2 - X_2 \\ S_3 &\gets Y_2 + X_2 \end{aligned}

\begin{aligned} S_4 &\gets S_0 S_2 \\ S_5 &\gets S_1 S_3 \\ S_6 &\gets Z_1 Z_2 \\ S_7 &\gets T_1 T_2 \end{aligned}

\begin{aligned} S_8 &\gets S_4 \cdot 121666 \\ S_9 &\gets S_5 \cdot 121666 \\ S_{10} &\gets S_6 \cdot 2 \cdot 121666 \\ S_{11} &\gets S_7 \cdot -2 \cdot 121665 \end{aligned}

\begin{aligned} S_{12} &\gets S_9 - S_8 \\ S_{13} &\gets S_9 + S_8 \\ S_{14} &\gets S_{10} - S_{11} \\ S_{15} &\gets S_{10} + S_{11} \end{aligned}

\begin{aligned} X_3 &\gets S_{12} S_{14} \\ Y_3 &\gets S_{15} S_{13} \\ Z_3 &\gets S_{15} S_{14} \\ T_3 &\gets S_{12} S_{13} \end{aligned}

to obtain \( P_3 = (X_3 : Y_3 : Z_3 : T_3) = P_1 + P_2 \).

Readdition

If the point \( P_2 = (X_2 : Y_2 : Z_2 : T_2) \) is fixed, we can precompute

\begin{aligned} S_2 &\gets Y_2 - X_2 \\ S_3 &\gets Y_2 + X_2 \end{aligned}

\begin{aligned} S_2' &\gets S_2 \cdot 121666 \\ S_3' &\gets S_3 \cdot 121666 \\ Z_2' &\gets Z_2 \cdot 2 \cdot 121666 \\ T_2' &\gets T_2 \cdot -2 \cdot 121665 \\ \end{aligned}

to obtain the CachedPoint \( (S_2', S_3', Z_2', T_2') \). This precomputation is essentially the same as that suggested in §3.1 of HWCD, with the difference that the multiplication by the curve constant \( -121665 / 121666 \) is spread over all four coordinates, to allow a vectorized computation of four multiplications of small constants instead of a serial computation of multiplication by a large constant.

To perform readdition of \(P_1 = (X_1 : Y_1 : Z_1 : T_1) \) and \(P_2 = (S_2', S_3', Z_2', T_2') \), we compute

\begin{aligned} S_0 &\gets Y_1 - X_1 \\ S_1 &\gets Y_1 + X_1 \end{aligned}

\begin{aligned} S_8 &\gets S_0 S_2' \\ S_9 &\gets S_1 S_3' \\ S_{10} &\gets Z_1 Z_2' \\ S_{11} &\gets T_1 T_2' \end{aligned}

\begin{aligned} S_{12} &\gets S_9 - S_8 \\ S_{13} &\gets S_9 + S_8 \\ S_{14} &\gets S_{10} - S_{11} \\ S_{15} &\gets S_{10} + S_{11} \end{aligned}

\begin{aligned} X_3 &\gets S_{12} S_{14} \\ Y_3 &\gets S_{15} S_{13} \\ Z_3 &\gets S_{15} S_{14} \\ T_3 &\gets S_{12} S_{13} \end{aligned}

to obtain \( P_3 = (X_3 : Y_3 : Z_3 : T_3) = P_1 + P_2 \).

Compared to the addition formulas above, this saves \( 1\mathbf D \).

Doubling

To double a point \( P = (X_1 : Y_1 : Z_1 : T_1) \), we compute

 S\_0 \gets X\_1 + Y\_1 

\begin{aligned} S_1 &\gets X_1^2 \\ S_2 &\gets Y_1^2 \\ S_3 &\gets Z_1^2 \\ S_4 &\gets S_0^2 \end{aligned}

\begin{aligned} S_5 &\gets S_1 + S_2 \\ S_6 &\gets S_1 - S_2 \\ S_7 &\gets 2S_3 \\ S_8 &\gets S_7 + S_6 = S_1 + 2S_3 - S_2 \\ S_9 &\gets S_5 - S_4 = S_1 + S_2 - S_4 \end{aligned}

\begin{aligned} X_3 &\gets S_8 S_9 \\ Y_3 &\gets S_5 S_6 \\ Z_3 &\gets S_8 S_6 \\ T_3 &\gets S_5 S_9 \end{aligned}

to obtain \( P_3 = (X_3 : Y_3 : Z_3 : T_3) = [2]P_1 \).

Unlike the (re)addition formulas, the divergent parts of these formulas are less nice. However, with some careful bounds-juggling, it is possible to implement them without inserting extra carry chains, as described below.

Field element representation

Our strategy is to implement 4-wide multiplication and squaring by wordslicing, using one 64-bit AVX2 lane for each field element. Field elements are represented in the usual way as 10 u32 limbs in radix \(25.5\) (i.e., alternating between \(2^{26}\) for even limbs and \(2^{25}\) for odd limbs). This has the effect that passing between the parallel 32-bit AVX2 representation and the serial 64-bit representation (which uses radix \(2^{51}\)) amounts to regrouping digits.

The field element representation is oriented around the AVX2 vpmuluqdq instruction, which multiplies the low 32 bits of each 64-bit lane of each operand to produce a 64-bit result.

(a1 ?? b1 ?? c1 ?? d1 ??)
(a2 ?? b2 ?? c2 ?? d2 ??)

(a1*a2 b1*b2 c1*c2 d1*d2)

To unpack 32-bit values into 64-bit lanes for use in multiplication it would be convenient to use the vpunpck[lh]dq instructions, which unpack and interleave the low and high 32-bit lanes of two source vectors. However, the AVX2 versions of these instructions are designed to operate only within 128-bit lanes of the 256-bit vectors, so that interleaving the low lanes of (a0 b0 c0 d0 a1 b1 c1 d1) with zero gives (a0 00 b0 00 a1 00 b1 00). Instead, we pre-shuffle the data layout as (a0 b0 a1 b1 c0 d0 c1 d1) so that we can unpack the "low" and "high" parts as

(a0 00 b0 00 c0 00 d0 00)
(a1 00 b1 00 c1 00 d1 00)

The data layout for a vector of four field elements \( (a,b,c,d) \) with limbs \( a_0, a_1, \ldots, a_9 \) is as [u32x8; 5] in the form

(a0 b0 a1 b1 c0 d0 c1 d1)
(a2 b2 a3 b3 c2 d2 c3 d3)
(a4 b4 a5 b5 c4 d4 c5 d5)
(a6 b6 a7 b7 c6 d6 c7 d7)
(a8 b8 a9 b9 c8 d8 c9 d9)

Since this breaks cleanly into two 128-bit lanes, it may be possible to adapt it to 128-bit vector instructions such as NEON without too much difficulty. Going the other direction, to extend this to AVX512, we could either run two point operations in parallel in lower and upper halves of the registers, or use 2-way parallelism within a field operation.

Handling the Doubling Formulas

The non-parallel portion of the doubling formulas is

\begin{aligned} S_5 &\gets S_1 + S_2 \\ S_6 &\gets S_1 - S_2 \\ S_7 &\gets 2S_3 \\ S_8 &\gets S_7 + S_6 = S_1 + 2S_3 - S_2 \\ S_9 &\gets S_5 - S_4 = S_1 + S_2 - S_4 \end{aligned}

Performing too many intermediate additions and subtractions grows the bounds beyond what is allowed as input to multiplication, forcing an extra carry pass. However, it is just possible to avoid this by rearranging signs.

Assume that the bounds on the limbs of each field element are parameterized by \( b \in \mathbb R \) representing the excess bits, so that each limb is bounded by either \( 2^{25+b} \) or \( 2^{26+b} \).

The multiplication routine requires that its inputs are bounded by \( b < 1.75 \), in order to fit a multiplication by \( 19 \) into 32 bits. Since \( \lg 19 < 4.25 \), \( 19x < 2^{32} \) when \( x < 2^{27.75} = 2^{26 + 1.75} \). However, this is only required for one of the inputs; the other can grow up to \( b < 2.5 \).

Computing \( (S_5, S_6, S_8, S_9 ) \) as

\begin{matrix} & S_1 & S_1 & S_1 & S_1 \\ +& S_2 & & & S_2 \\ +& & & S_3 & \\ +& & & S_3 & \\ +& & 2p & 2p & 2p \\ -& & S_2 & S_2 & \\ -& & & & S_4 \\ =& S_5 & S_6 & S_8 & S_9 \end{matrix}

results in bit-excesses \( (1.00, 1.59, 2.33, 2.00)\) for \( (S_5, S_6, S_8, S_9 ) \). The products we want to compute are then

\begin{aligned} X_3 &\gets S_8 S_9 \leftrightarrow (2.33, 2.00) \\ Y_3 &\gets S_5 S_6 \leftrightarrow (1.00, 1.59) \\ Z_3 &\gets S_8 S_6 \leftrightarrow (2.33, 1.59) \\ T_3 &\gets S_5 S_9 \leftrightarrow (1.00, 2.00) \end{aligned}

which are too large. However, if we flip the sign of \( S_4 = S_0^2 \) during squaring, so that we output \(S_4' = -S_4 \pmod p\), then we can compute

\begin{matrix} & S_1 & S_1 & S_1 & S_1 \\ +& S_2 & & & S_2 \\ +& & & S_3 & \\ +& & & S_3 & \\ +& & & & S_4' \\ +& & 2p & 2p & \\ -& & S_2 & S_2 & \\ =& S_5 & S_6 & S_8 & S_9 \end{matrix}

resulting in bit-excesses \( (1.00, 1.59, 2.33, 1.59)\) for \( (S_5, S_6, S_8, S_9 ) \). The products we want to compute are then

\begin{aligned} X_3 &\gets S_8 S_9 \leftrightarrow (2.33, 1.59) \\ Y_3 &\gets S_5 S_6 \leftrightarrow (1.00, 1.59) \\ Z_3 &\gets S_8 S_6 \leftrightarrow (2.33, 1.59) \\ T_3 &\gets S_5 S_9 \leftrightarrow (1.00, 1.59) \end{aligned}

whose right-hand sides are all bounded with \( b < 1.75 \) and whose left-hand sides are all bounded with \( b < 2.5 \).

Comparison to non-vectorized formulas

HWCD also suggest using a mixed representation, passing between \( \mathbb P^3 \) "extended" coordinates and \( \mathbb P^2 \) "projective" coordinates, where doubling is slightly cheaper (saving about \(\mathbf 1M\). This approach is used for the non-vectorized u32 and u64 backends, and more details on the different coordinate systems can be found in the curve_models module documentation.

This optimization is not compatible with the parallel formulas, which are therefore slightly less efficient when counting the total number of field multiplications and squarings. In particular, vectorized doublings are less efficient than serial doublings.

In addition, the parallel formulas can only use a \( 32 \times 32 \rightarrow 64 \)-bit integer multiplier, so the speedup from vectorization must overcome the disadvantage of losing the \( 64 \times 64 \rightarrow 128\)-bit (serial) integer multiplier.

When compiling with AVX512VL, LLVM is able to use the extra ymm16..ymm31 registers to reduce register pressure, and avoid spills during field multiplication. This gives a small but noticeable speedup.

Another concern with AVX2 is that currently-available Intel processors (particularly Skylake and Skylake-X microarchitectures) perform thermal throttling when using wide vector instructions. For a mixed workload, where point operations are interspersed with other tasks, this can reduce overall performance. This probably means that this implementation is not suitable for basic applications, like signatures, but could still be worthwhile for complex applications, like zero-knowledge proofs, which do enough work to make it worthwhile.

On AMD's Zen microarchitecture, thermal throttling is not a concern, since AVX2 is implemented at half rate, so there is no penalty for mixed workloads (but also no speedup).