diff --git a/src/backend/avx2/mod.rs b/src/backend/avx2/mod.rs index 456c451..7a06512 100644 --- a/src/backend/avx2/mod.rs +++ b/src/backend/avx2/mod.rs @@ -8,6 +8,305 @@ // - Isis Agora Lovecruft // - Henry de Valence +//! 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. +//! +//! 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. 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 this can be done in parallel as a scaling by \\( (121666, 121666, +//! 2\cdot 121665, 2\cdot 121666) \\). To handle the sign, we use +//! masking to negate one of the field elements. +//! +//! 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. +//! +//! 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; this strategy was +//! used by Tung Chou's `sandy2x` implementation, which used a 2-wide +//! field implementation in 128-bit registers. Curiously, however, +//! although the `sandy2x` paper cites the HWCD paper for extended +//! twisted Edwards coordinates, it does not mention the 4-wide HWCD +//! Edwards formulas or that the 2-wide Montgomery formulas it uses were +//! previously published there. +//! +//! 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 used for the parallel formulas, which are +//! therefore slightly less efficient when counting the total number of +//! multiplications and squarings. In addition, the parallel formulas +//! can only use a \\( 32 \times 32 \rightarrow 64 \\)-bit multiplier +//! instead of a \\( 64 \times 64 \rightarrow 128\\)-bit multiplier. +//! +//! 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). +//! +//! # 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 \\). +//! +//! ## 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 \\). +//! +//! In practice, we compute \\( (S\_5, S\_6, S\_7, 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} +//! $$ +//! +//! adding multiples of \\(p\\) to prevent underflow. This results in +//! 32-bit limbs which are just too large for multiplication, so we +//! perform a reduction. However, since we just need to reduce the +//! excess in each limb, not a full reduction, it's enough to perform +//! each carry in parallel. +//! +//! With some finesse, it may be possible to rearrange this +//! computation to avoid the extra carry pass, but this is not yet +//! implemented. +//! +//! # 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. +//! +//! 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. + + pub(crate) mod field; pub(crate) mod edwards;