betrusted-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_][hwcd08], 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.

Currently, this implementation uses only the first two algorithms.

# Parallel formulas

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.

# Implementation strategy

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 share a single
instruction.

Our strategy is to implement 4-wide multiplication and squaring
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 amounts to regrouping digits.

The addition and subtraction steps are done largely serially, using
masking to handle the instruction divergence.  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).

The 4-wide formulas of the HWCD paper do not seem to have been
implemented using SIMD before.  The 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][sandy2x] also
implements Edwards arithmetic, and cites the HWCD paper, it doesn't
mention or discuss the parallel formulas from HWCD, or that the
2-wide Montgomery formulas it uses were previously published there.
There is also a 2015 paper by Hernández and López on 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][avx2trac], suggesting that it also overlooked the
HWCD 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.

# 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 (SSA) form.

## Addition

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 \\).

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} \\) or \\(
2\^{26} \\).

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 \\).

# Field element representation

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.

```text,no_run
(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

```text,no_run
(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

```text,no_run
(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.

We don't attempt to use AVX2 for serial field element computations
such as inversion, since wherever we have AVX2 we also have `mulx`.
However, it might be useful for batched inverse square-root
computations, which can't be batched in the same way inversions can.

# Implementation details

The implementation uses the unstable `stdsimd` crate to provide AVX2
intrinsics, and the code is not yet cleanly factored between the
field element parts and the point parts.

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.

The addition and subtraction steps involve masking, to apply
operations to a single lane of the vector.  AVX512VL extends the
predication features of AVX512 to AVX2 code and would probably be
beneficial.  Unfortunately, LLVM is currently unable to lower `op +
blend` into an AVX512VL masked operation.  However, the explicitly
masked versions of the intrinsics seem to produce the same LLVM IR
as an `op + blend`, so hopefully this will improve as the AVX512
support in LLVM improves.

When used for constant-time variable-base scalar multiplication,
this strategy (using AVX2) gives a significant speedup over the
serial implementation (using the \\(64 \times 64\\) multiplier) of
approximately 1.6x for Skylake-X with `target_cpu=skylake` (using AVX2), of
approximately 1.8x for Skylake-X with `target_cpu=skylake-avx512` (using the extra
`ymm16..ymm31` registers from AVX512VL), and of approximately 1.0x
for Ryzen (which implements AVX2 at half rate).

When used for variable-time double-base scalar multiplication 
\\( aA + bB \\) for fixed \\(B\\) (as in, e.g., signature verification),
this strategy provides a 1.4x speedup on Skylake-X over the same
operation as implemented in `ed25519-donna`, the fastest
production-quality Ed25519 implementation.

[sandy2x]: https://eprint.iacr.org/2015/943.pdf
[avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28
[hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf