diff --git a/.travis.yml b/.travis.yml index cf6db51..50e9a3c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -10,7 +10,7 @@ env: # Tests the u64 backend - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std u64_backend' # Tests the avx2 backend - - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend yolocrypto' + - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend' # Tests serde support and default feature selection - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='serde' # Tests building without std. We have to select a backend, so we select the one @@ -21,7 +21,7 @@ matrix: exclude: # Test the avx2 backend only on nightly - rust: stable - env: TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend yolocrypto' + env: TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend' # Test no_std only on nightly. - rust: stable env: TEST_COMMAND=build EXTRA_FLAGS=--no-default-features FEATURES='u32_backend' diff --git a/Makefile b/Makefile index 8b376cb..e447175 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -FEATURES := nightly yolocrypto +FEATURES := nightly yolocrypto avx2_backend doc: cargo rustdoc --features "$(FEATURES)" -- --html-in-header docs/assets/rustdoc-include-katex-header.html diff --git a/README.md b/README.md index 7dba401..0943fd7 100644 --- a/README.md +++ b/README.md @@ -71,15 +71,17 @@ Curve arithmetic is implemented using one of the following backends: * a `u32` backend using `u64` products; * a `u64` backend using `u128` products; -* an experimental AVX2 backend, available using the `yolocrypto` feature when - compiling for a target with `target_feature=+avx2`. +* an `avx2` backend using parallel formulas, available when compiling for a + target with `target_feature=+avx2`. By default the `u64` backend is selected. To select a specific backend, use: ```sh cargo build --no-default-features --features "std u32_backend" cargo build --no-default-features --features "std u64_backend" -cargo build --no-default-features --features "std avx2_backend yolocrypto" +cargo build --no-default-features --features "std avx2_backend" ``` +Crates using `curve25519-dalek` can either select a backend on behalf of their +users, or expose feature flags that control the `curve25519-dalek` backend. Benchmarks are run using [`criterion.rs`][criterion]: @@ -88,7 +90,7 @@ Benchmarks are run using [`criterion.rs`][criterion]: export RUSTFLAGS="-C target_cpu=native" cargo bench --no-default-features --features "std u32_backend" cargo bench --no-default-features --features "std u64_backend" -cargo bench --no-default-features --features "std avx2_backend yolocrypto" +cargo bench --no-default-features --features "std avx2_backend" ``` # Contributing @@ -117,7 +119,8 @@ to the Dalek race.* Portions of this library were originally a port of [Adam Langley's Golang ed25519 library](https://github.com/agl/ed25519), which was in -turn a port of the reference `ref10` implementation. +turn a port of the reference `ref10` implementation. Most of this code, +including the 32-bit field arithmetic, has since been rewritten. The fast `u32` and `u64` scalar arithmetic was implemented by Andrew Moon, and the addition chain for scalar inversion was provided by Brian Smith. diff --git a/build.rs b/build.rs index 64aec5d..e0a2da7 100644 --- a/build.rs +++ b/build.rs @@ -1,5 +1,5 @@ #![cfg_attr(feature = "nightly", feature(cfg_target_feature))] -#![cfg_attr(all(feature = "nightly", feature = "yolocrypto"), feature(stdsimd))] +#![cfg_attr(all(feature = "nightly", feature = "avx2_backend"), feature(stdsimd))] #![allow(unused_variables)] #![allow(non_snake_case)] #![allow(dead_code)] diff --git a/docs/avx2-notes.md b/docs/avx2-notes.md index fb8085d..d902650 100644 --- a/docs/avx2-notes.md +++ b/docs/avx2-notes.md @@ -1,42 +1,48 @@ -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 vectorized implementation of group operations on the twisted Edwards +form of Curve25519, using a modification of the 4-way parallel +formulas of Hisil, Wong, Carter, and Dawson. -* a unified addition algorithm taking an effective \\(2\mathbf M + -1\mathbf D\\); +# Overview -* a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf -S\\); +The 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08] by Hisil, +Wong, Carter, and Dawson (HWCD) introduced the “extended coordinates” +and mixed-model representations which are used by most Edwards curve +implementations. -* a dedicated (i.e., for distinct points) addition algorithm taking -an effective \\(2 \mathbf M \\). +However, they also describe 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\\), and a dedicated (i.e., for +distinct points) addition algorithm taking an effective \\(2 \mathbf M +\\). They compare these formulas with a 2-way parallel variant of the +Montgomery ladder. -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. +Unlike their serial formulas, which are used widely, their parallel +formulas do not seem to have been implemented in software before. The +2-way parallel Montgomery ladder was used in 2015 by Tung Chou's +`sandy2x` implementation. Curiously, however, although the [`sandy2x` +paper][sandy2x] also implements Edwards arithmetic, and cites HWCD08, +it doesn't mention their parallel Edwards formulas. +A 2015 paper by Hernández and López describes an AVX2 implementation +of X25519. Neither the paper nor the code are publicly available, but +it apparently gives only a [slight speedup][avx2trac], suggesting that +it uses a 4-way parallel Montgomery ladder rather than parallel +Edwards formulas. -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][avx2trac]. -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][sandy2x] 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 reason may be that HWCD08 describe their formulas as operating on +four independent processors, which would make a software +implementation impractical: all of the operations are too low-latency +to effectively synchronize. But a closer inspection reveals that the +(more expensive) multiplication and squaring steps are uniform, while +the instruction divergence occurs in the (much cheaper) addition and +subtraction steps. This means that a SIMD implementation can perform +the expensive steps uniformly, and handle divergence in the +inexpensive steps using masking. -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`. +These notes describe modifications to the original parallel formulas +to allow a SIMD implementation, and this module contains an +implementation of the modified formulas using 256-bit AVX2 vector +operations. # Parallel formulas in HWCD'08 @@ -60,218 +66,153 @@ and the unified addition algorithm is presented as follows: | | \\( 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. +Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of +multiplication and squaring of generic field elements, \\(\mathbf D\\) +represents the cost of multiplication by a curve constant (in this +case \\( k = 2d \\)). -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. +Notice that the \\(1\mathbf M\\) and \\(1\mathbf S\\) steps are +uniform. The non-uniform steps are all inexpensive additions or +subtractions, with the exception of the multiplication by the curve +constant \\(k = 2d\\): +$$ +R\_7 \gets 2 d R\_7. +$$ -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 +HWCD suggest parallelising this step by breaking \\(k = 2d\\) 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. +computing \\(k_i R_7 \\) in parallel. This is quite awkward, but if +the curve constant is a ratio \\( d = d\_1/d\_2 \\), then projective +coordinates allow us to instead compute +$$ +(R\_5, R\_6, R\_7, R\_8) \gets (d\_2 R\_5, d\_2 R\_6, 2d\_1 R\_7, d\_2 R\_8). +$$ +This can be performed as a uniform multiplication by a vector of +constants, and if \\(d\_1, d\_2\\) are small, it is relatively +inexpensive. (This trick was suggested by Mike Hamburg). +In the Curve25519 case, we have +$$ +d = \frac{d\_1}{d\_2} = \frac{-121665}{121666}; +$$ +Since \\(2 \cdot 121666 < 2\^{18}\\), all the constants above fit (up +to sign) in 32 bits, so this can be done in parallel as four +multiplications by small constants \\( (121666, 121666, 2\cdot 121665, +2\cdot 121666) \\), followed by a negation to compute \\( - 2\cdot 121665\\). -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) \\). +# Modified parallel formulas -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. +Using the modifications sketched above, we can write SIMD-friendly +versions of the parallel formulas as follows. To avoid confusion with +the original 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 - +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 +(S\_0 &&,&& S\_1 &&,&& S\_2 &&,&& S\_3 ) +&\gets +(Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Y\_2 - X\_2&&,&& Y\_2 + X\_2) +\\\\ +(S\_4 &&,&& S\_5 &&,&& S\_6 &&,&& S\_7 ) +&\gets +(S\_0 \cdot S\_2&&,&& S\_1 \cdot S\_3&&,&& Z\_1 \cdot Z\_2&&,&& T\_1 \cdot T\_2) +\\\\ +(S\_8 &&,&& S\_9 &&,&& S\_{10} &&,&& S\_{11} ) +&\gets +(d\_2 \cdot S\_4 &&,&& d\_2 \cdot S\_5 &&,&& 2 d\_2 \cdot S\_6 &&,&& 2 d\_1 \cdot S\_7 ) +\\\\ +(S\_{12} &&,&& S\_{13} &&,&& S\_{14} &&,&& S\_{15}) +&\gets +(S\_9 - S\_8&&,&& S\_9 + S\_8&&,&& S\_{10} - S\_{11}&&,&& S\_{10} + S\_{11}) +\\\\ +(X\_3&&,&& Y\_3&&,&& Z\_3&&,&& T\_3) +&\gets +(S\_{12} \cdot S\_{14}&&,&& S\_{15} \cdot S\_{13}&&,&& S\_{15} \cdot S\_{14}&&,&& S\_{12} \cdot S\_{13}) \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 \\). +This costs \\( 2\mathbf M + 1 \mathbf D\\). ## Readdition -If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we can precompute - +If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we +can cache the multiplication of the curve constants by computing $$ \begin{aligned} -S\_2 &\gets Y\_2 - X\_2 \\\\ -S\_3 &\gets Y\_2 + X\_2 +(S\_2' &&,&& S\_3' &&,&& Z\_2' &&,&& T\_2' ) +&\gets +(d\_2 \cdot (Y\_2 - X\_2)&&,&& d\_2 \cdot (Y\_1 + X\_1)&&,&& 2d\_2 \cdot Z\_2 &&,&& 2d\_1 \cdot T\_2). \end{aligned} $$ - +This costs \\( 1\mathbf D\\); with \\( (S\_2', S\_3', Z\_2', T\_2')\\) +in hand, the addition formulas above become $$ \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 \\\\ +(S\_0 &&,&& S\_1 &&,&& Z\_1 &&,&& T\_1 ) +&\gets +(Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Z\_1 &&,&& T\_1) +\\\\ +(S\_8 &&,&& S\_9 &&,&& S\_{10} &&,&& S\_{11} ) +&\gets +(S\_0 \cdot S\_2' &&,&& S\_1 \cdot S\_3'&&,&& Z\_1 \cdot Z\_2' &&,&& T\_1 \cdot T\_2') +\\\\ +(S\_{12} &&,&& S\_{13} &&,&& S\_{14} &&,&& S\_{15}) +&\gets +(S\_9 - S\_8&&,&& S\_9 + S\_8&&,&& S\_{10} - S\_{11}&&,&& S\_{10} + S\_{11}) +\\\\ +(X\_3&&,&& Y\_3&&,&& Z\_3&&,&& T\_3) +&\gets +(S\_{12} \cdot S\_{14}&&,&& S\_{15} \cdot S\_{13}&&,&& S\_{15} \cdot S\_{14}&&,&& S\_{12} \cdot S\_{13}) \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 \\). +which costs only \\( 2\mathbf M \\). This precomputation is +essentially similar to the precomputation that HWCD suggest for their +serial formulas. Because the cost of precomputation and then +readdition is the same as addition, it's sufficient to only +implement caching and readdition. ## Doubling +The non-uniform portions of the (re)addition formulas have a fairly +regular structure. Unfortunately, this is not the case for the +doubling formulas, which are much less nice. + 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 +(X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& S\_0) +&\gets +(X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& X\_1 + Y\_1) +\\\\ +(S\_1 &&,&& S\_2 &&,&& S\_3 &&,&& S\_4 ) +&\gets +(X\_1\^2 &&,&& Y\_1\^2&&,&& Z\_1\^2 &&,&& S\_0\^2) +\\\\ +(S\_5 &&,&& S\_6 &&,&& S\_8 &&,&& S\_9 ) +&\gets +(S\_1 + S\_2 &&,&& S\_1 - S\_2 &&,&& S\_1 + 2S\_3 - S\_2 &&,&& S\_1 + S\_2 - S\_4) +\\\\ +(X\_3 &&,&& Y\_3 &&,&& Z\_3 &&,&& T\_3 ) +&\gets +(S\_8 \cdot S\_9 &&,&& S\_5 \cdot S\_6 &&,&& S\_8 \cdot S\_6 &&,&& S\_5 \cdot S\_9) \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. +The intermediate step between the squaring and multiplication requires +a long chain of additions, but with some care and finesse, +described below, it is possible (in our case) to arrange this +computation without requiring an intermediate reduction. + +However, it does mean that the doubling formulas have proportionately +more vectorization overhead than the (re)addition formulas. The +effects of this are discussed in the comparison section below. # Field element representation @@ -329,39 +270,36 @@ 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 +# Avoiding Overflow in Doubling -The non-parallel portion of the doubling formulas is +To analyze the size of the field element coefficients during the +computations, we can parameterize the bounds on the limbs of each +field element by \\( b \in \mathbb R \\) representing the excess bits +above that limb's radix, so that each limb is bounded by either +\\(2\^{25+b} \\) or \\( 2\^{26+b} \\), as appropriate. -$$ -\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 +The multiplication routine requires that its inputs are bounded with \\( 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 +In addition, the multiplication and squaring routines do not +canonically reduce their outputs, but can leave some small uncarried +excesses, so that their reduced outputs are bounded with +\\( b < 0.007 \\). +The non-parallel portion of the doubling formulas is +$$ +\begin{aligned} +(S\_5 &&,&& S\_6 &&,&& S\_8 &&,&& S\_9 ) +&\gets +(S\_1 + S\_2 &&,&& S\_1 - S\_2 &&,&& S\_1 + 2S\_3 - S\_2 &&,&& S\_1 + S\_2 - S\_4) +\end{aligned} +$$ + +Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as $$ \begin{matrix} & S\_1 & S\_1 & S\_1 & S\_1 \\\\ @@ -374,24 +312,22 @@ $$ =& S\_5 & S\_6 & S\_8 & S\_9 \end{matrix} $$ - -results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for +results in bit-excesses \\( < (1.01, 1.60, 2.33, 2.01)\\) 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) +X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.01) \\\\ +Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\ +Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\ +T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 2.01) \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 - +which are too large: it's not possible to arrange the multiplicands so +that one vector has \\(b < 2.5\\) and the other has \\( b < 1.75 \\). +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 \\\\ @@ -404,61 +340,120 @@ $$ =& S\_5 & S\_6 & S\_8 & S\_9 \end{matrix} $$ - -resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for +resulting in bit-excesses \\( < (1.01, 1.60, 2.33, 1.60)\\) 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) +X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.60) \\\\ +Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\ +Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\ +T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 1.60) \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 \\). +whose left-hand sides are all bounded with \\( b < 2.5 \\), +so that we can avoid any intermediate reductions. # 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. +In theory, the parallel Edwards formulas seem to allow a \\(4\\)-way +speedup from parallelism. However, an actual vectorized +implementation has several slowdowns that cut into this speedup. -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 +First, 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. +\times 64 \rightarrow 128\\)-bit (serial) integer multiplier. The +effect of this slowdown is microarchitecture-dependent, since it +requires accounting for the total number of multiplications and +additions and their relative costs. In the future, it will probably +be possible to avoid this slowdown by using the `IFMA52` instructions, +whose parallelism is perfectly suited to these formulas. -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. +Second, the parallel doubling formulas incur both a theoretical and +practical slowdown. The parallel formulas described above work on the +\\( \mathbb P\^3 \\) “extended” coordinates. The \\( \mathbb P\^2 \\) +model introduced earlier by [Bernstein, Birkner, Joye, Lange, and +Peters][bbjlp08] allows slightly faster doublings, so HWCD suggest +mixing coordinate systems while performing scalar multiplication +(attributing the idea to [a 1998 paper][cmo98] by Cohen, Miyagi, and +Ono). The \\( T \\) coordinate is not required for doublings, so when +doublings are followed by doublings, its computation can be skipped. +More details on this approach and the different coordinate systems can +be found in the [`curve_models` module documentation][curve_models]. -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, +Unfortunately, this optimization is not compatible with the parallel +formulas, which cannot save time by skipping a single variable, so the +parallel doubling formulas do slightly more work when counting the +total number of field multiplications and squarings. + +In addition, the parallel doubling formulas have a less regular +pattern of additions and subtractions than the parallel addition +formulas, so the vectorization overhead is proportionately greater. +Both the parallel addition and parallel doubling formulas also require +some shuffling to rearrange data within the vectors, which places more +pressure on the shuffle unit than is desirable. + +This means that the speedup from using a vectorized implementation of +parallel Edwards formulas is likely to be greatest in applications +that do fewer doublings and more additions (like a large multiscalar +multiplication) rather than applications that do fewer additions and +more doublings (like a double-base scalar multiplication). + +Third, current Intel CPUs perform thermal throttling when using wide +vector instructions. A detailed description can be found in §15.26 of +[the Intel Optimization Manual][intel], but using wide vector +instructions prevents the core from operating at higher frequencies. +The core can return to the higher-frequency state after 2 +milliseconds, but this timer is reset every time high-power +instructions are used. + +Any speedup from vectorization therefore has to be weighed against a +slowdown for the next few million 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. +reduce overall performance. This implementation is therefore probably +not suitable for basic applications, like signatures, but is +worthwhile for complex applications, like zero-knowledge proofs, which +do sustained work. + +For this reason, the AVX2 backend is not enabled by default, but can +be selected using the `avx2_backend` feature. + +# Future work + +There are several directions for future improvement: + +* Using the vectorized field arithmetic code to parallelize across + point operations rather than within a single point operation. This + is less flexible, but would give a speedup both from allowing use of + the faster mixed-model arithmetic and from reducing shuffle + pressure. One approach in this direction would be to implement + batched scalar-point operations using vectors of points (AoSoA + layout). This less generally useful but would give a speedup for + Bulletproofs. + +* Extending the implementation to use the full width of AVX512, either + handling the extra parallelism internally to a single point + operation (by using a 2-way parallel implementation of field + arithmetic instead of a wordsliced one), or externally, + parallelizing across point operations. Internal parallelism would + be preferable but might require too much shuffle pressure. + +* Generalizing the implementation to non-AVX2 instructions, + particularly NEON. The current point arithmetic code is written in + terms of field element vectors, which are in turn implemented using + platform SIMD vectors. It should be possible to write an alternate + implementation of the `FieldElement32x4` using NEON without changing + the point arithmetic. NEON has 128-bit vectors rather than 256-bit + vectors, but this may still be worthwhile compared to a serial + implementation. -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). [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 +[curve_models]: https://doc-internal.dalek.rs/curve25519_dalek/curve_models/index.html +[bbjlp08]: https://eprint.iacr.org/2008/013 +[cmo98]: https://link.springer.com/content/pdf/10.1007%2F3-540-49649-1_6.pdf +[intel]: https://software.intel.com/sites/default/files/managed/9e/bc/64-ia-32-architectures-optimization-manual.pdf \ No newline at end of file diff --git a/src/backend/avx2/constants.rs b/src/backend/avx2/constants.rs index e99e312..304a1be 100644 --- a/src/backend/avx2/constants.rs +++ b/src/backend/avx2/constants.rs @@ -16,6 +16,24 @@ use backend::avx2::edwards::{CachedPoint, ExtendedPoint}; use backend::avx2::field::FieldElement32x4; use scalar_mul::window::NafLookupTable8; +/// The identity element as an `ExtendedPoint`. +pub(crate) static EXTENDEDPOINT_IDENTITY: ExtendedPoint = ExtendedPoint(FieldElement32x4([ + u32x8::new(0, 1, 0, 0, 1, 0, 0, 0), + u32x8::splat(0), + u32x8::splat(0), + u32x8::splat(0), + u32x8::splat(0), +])); + +/// The identity element as a `CachedPoint`. +pub(crate) static CACHEDPOINT_IDENTITY: CachedPoint = CachedPoint(FieldElement32x4([ + u32x8::new(121647, 121666, 0, 0, 243332, 67108845, 0, 33554431), + u32x8::new(67108864, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), +])); + /// The low limbs of (2p, 2p, 2p, 2p), so that /// ```no_run /// (2p, 2p, 2p, 2p) = [P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI] @@ -76,14 +94,6 @@ pub(crate) static P_TIMES_16_HI: u32x8 = u32x8::new( 33554431 << 4, ); -pub(crate) static P_TIMES_2_MASKED: FieldElement32x4 = FieldElement32x4([ - u32x8::new(0, 134217690, 0, 67108862, 134217690, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), -]); - /// Odd multiples of the Ed25519 basepoint: pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: NafLookupTable8 = NafLookupTable8([ CachedPoint(FieldElement32x4([ diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index e317b65..43f02f7 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -8,29 +8,53 @@ // - Isis Agora Lovecruft // - Henry de Valence -//! Extended Twisted Edwards for Curve25519, using AVX2. +//! Parallel Edwards Arithmetic for Curve25519. +//! +//! This module currently has two point types: +//! +//! * `ExtendedPoint`: a point stored in vector-friendly format, with +//! vectorized doubling and addition; +//! +//! * `CachedPoint`: used for readdition. +//! +//! Details on the formulas can be found in the documentation for the +//! parent `avx2` module. +//! +//! This API is designed to be safe: vectorized points can only be +//! created from serial points (which do validation on decompression), +//! and operations on valid points return valid points, so invalid +//! point states should be unrepresentable. +//! +//! This design goal is met, with one exception: the `Neg` +//! implementation for the `CachedPoint` performs a lazy negation, so +//! that subtraction can be efficiently implemented as a negation and +//! an addition. Repeatedly negating a `CachedPoint` will cause its +//! coefficients to grow and eventually overflow. Repeatedly negating +//! a point should not be necessary anyways. -// just going to own it -#![allow(bad_style)] +#![allow(non_snake_case)] use core::convert::From; -use core::ops::{Add, Sub, Neg}; +use core::ops::{Add, Neg, Sub}; -use core::simd::{IntoBits, u32x8}; - -use subtle::ConditionallyAssignable; use subtle::Choice; +use subtle::ConditionallyAssignable; use edwards; use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use traits::Identity; -use backend::avx2::field::{D_LANES, Lanes, FieldElement32x4}; +use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle}; +use backend::avx2::constants; -use backend::avx2; - -/// A point on Curve25519, represented in an AVX2-friendly format. +/// A point on Curve25519, using parallel Edwards formulas for curve +/// operations. +/// +/// # Invariant +/// +/// The coefficients of an `ExtendedPoint` are bounded with +/// \\( b < 0.007 \\). #[derive(Copy, Clone, Debug)] pub struct ExtendedPoint(pub(super) FieldElement32x4); @@ -43,7 +67,12 @@ impl From for ExtendedPoint { impl From for edwards::EdwardsPoint { fn from(P: ExtendedPoint) -> edwards::EdwardsPoint { let tmp = P.0.split(); - edwards::EdwardsPoint{X: tmp[0], Y: tmp[1], Z: tmp[2], T: tmp[3]} + edwards::EdwardsPoint { + X: tmp[0], + Y: tmp[1], + Z: tmp[2], + T: tmp[3], + } } } @@ -61,111 +90,71 @@ impl Default for ExtendedPoint { impl Identity for ExtendedPoint { fn identity() -> ExtendedPoint { - ExtendedPoint(FieldElement32x4([ - u32x8::new(0,1,0,0,1,0,0,0), - u32x8::splat(0), - u32x8::splat(0), - u32x8::splat(0), - u32x8::splat(0), - ])) + constants::EXTENDEDPOINT_IDENTITY } } impl ExtendedPoint { + /// Compute the double of this point. pub fn double(&self) -> ExtendedPoint { - unsafe { - use core::arch::x86_64::_mm256_permute2x128_si256; - use core::arch::x86_64::_mm256_permutevar8x32_epi32; - use core::arch::x86_64::_mm256_blend_epi32; - use core::arch::x86_64::_mm256_shuffle_epi32; + // Want to compute (X1 Y1 Z1 X1+Y1). + // Not sure how to do this less expensively than computing + // (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1) + // (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1) + // and then adding. - let P = &self.0; + // Set tmp0 = (X1 Y1 X1 Y1) + let mut tmp0 = self.0.shuffle(Shuffle::ABAB); - let mut t0 = FieldElement32x4::zero(); - let mut t1 = FieldElement32x4::zero(); + // Set tmp1 = (Y1 X1 Y1 X1) + let mut tmp1 = tmp0.shuffle(Shuffle::BADC); - // Want to compute (X1 Y1 Z1 X1+Y1). - // Not sure how to do this less expensively than computing - // (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1) - // (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1) - // and then adding. + // Set tmp0 = (X1 Y1 Z1 X1+Y1) + tmp0 = self.0.blend(tmp0 + tmp1, Lanes::D); - // Set t0 = (X1 Y1 X1 Y1) - t0.0[0] = _mm256_permute2x128_si256(P.0[0].into_bits(), P.0[0].into_bits(), 0b0000_0000).into_bits(); - t0.0[1] = _mm256_permute2x128_si256(P.0[1].into_bits(), P.0[1].into_bits(), 0b0000_0000).into_bits(); - t0.0[2] = _mm256_permute2x128_si256(P.0[2].into_bits(), P.0[2].into_bits(), 0b0000_0000).into_bits(); - t0.0[3] = _mm256_permute2x128_si256(P.0[3].into_bits(), P.0[3].into_bits(), 0b0000_0000).into_bits(); - t0.0[4] = _mm256_permute2x128_si256(P.0[4].into_bits(), P.0[4].into_bits(), 0b0000_0000).into_bits(); + // Set tmp1 = tmp0^2, negating the D values + tmp1 = tmp0.square_and_negate_D(); + // Now tmp1 = (S1 S2 S3 -S4) with b < 0.007 - // Set t1 = (Y1 X1 Y1 X1) - t1.0[0] = _mm256_shuffle_epi32(t0.0[0].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[1] = _mm256_shuffle_epi32(t0.0[1].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[2] = _mm256_shuffle_epi32(t0.0[2].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[3] = _mm256_shuffle_epi32(t0.0[3].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[4] = _mm256_shuffle_epi32(t0.0[4].into_bits(), 0b10_11_00_01).into_bits(); + // See discussion of bounds in the module-level documentation. + // We want to compute + // + // + | S1 | S1 | S1 | S1 | + // + | S2 | | | S2 | + // + | | | S3 | | + // + | | | S3 | | + // + | | | |-S4 | + // + | | 2p | 2p | | + // - | | S2 | S2 | | + // ======================= + // S5 S6 S8 S9 - // Set t0 = (X1+Y1 X1+Y1 X1+Y1 X1+Y1) - t0.0[0] = t0.0[0] + t1.0[0]; - t0.0[1] = t0.0[1] + t1.0[1]; - t0.0[2] = t0.0[2] + t1.0[2]; - t0.0[3] = t0.0[3] + t1.0[3]; - t0.0[4] = t0.0[4] + t1.0[4]; + let zero = FieldElement32x4::zero(); + let S_1 = tmp1.shuffle(Shuffle::AAAA); + let S_2 = tmp1.shuffle(Shuffle::BBBB); - // Set t0 = (X1 Y1 Z1 X1+Y1) - // why does this intrinsic take an i32 for the imm8 ??? - t0.0[0] = _mm256_blend_epi32(P.0[0].into_bits(), t0.0[0].into_bits(), D_LANES as i32).into_bits(); - t0.0[1] = _mm256_blend_epi32(P.0[1].into_bits(), t0.0[1].into_bits(), D_LANES as i32).into_bits(); - t0.0[2] = _mm256_blend_epi32(P.0[2].into_bits(), t0.0[2].into_bits(), D_LANES as i32).into_bits(); - t0.0[3] = _mm256_blend_epi32(P.0[3].into_bits(), t0.0[3].into_bits(), D_LANES as i32).into_bits(); - t0.0[4] = _mm256_blend_epi32(P.0[4].into_bits(), t0.0[4].into_bits(), D_LANES as i32).into_bits(); + tmp0 = zero.blend(tmp1 + tmp1, Lanes::C); + // tmp0 = (0, 0, 2S_3, 0) + tmp0 = tmp0.blend(tmp1, Lanes::D); + // tmp0 = (0, 0, 2S_3, -S_4) + tmp0 = tmp0 + S_1; + // tmp0 = ( S_1, S_1, S_1 + 2S_3, S_1 - S_4) + tmp0 = tmp0 + zero.blend(S_2, Lanes::AD); + // tmp0 = (S_1 + S_2, S_1, S_1 + 2S_3, S_1 + S_2 - S_4) + tmp0 = tmp0 + zero.blend(S_2.negate_lazy(), Lanes::BC); + // tmp0 = (S_1 + S_2, S_1 - S_2, S_1 - S_2 + 2S_3, S_1 + S_2 - S_4) + // b < ( 1.01, 1.6, 2.33, 1.6) + // Now tmp0 = (S_5, S_6, S_8, S_9) - // Set t1 = t0^2, negating the D values - t1 = t0.square_and_negate_D(); + // Set tmp1 = ( S_9, S_6, S_6, S_9) + // b < ( 1.6, 1.6, 1.6, 1.6) + tmp1 = tmp0.shuffle(Shuffle::DBBD); + // Set tmp1 = ( S_8, S_5, S_8, S_5) + // b < (2.33, 1.01, 2.33, 1.01) + tmp0 = tmp0.shuffle(Shuffle::CACA); - // Now t1 = (S1 S2 S3 -S4) - - let c0 = u32x8::new(0,0,2,2,0,0,2,2).into_bits(); // (ABCD) -> (AAAA) - let c1 = u32x8::new(1,1,3,3,1,1,3,3).into_bits(); // (ABCD) -> (BBBB) - - // See discussion of bounds in the module-level documentation. - // - // We want to compute - // - // + | S1 | S1 | S1 | S1 | - // + | S2 | | | S2 | - // + | | | S3 | | - // + | | | S3 | | - // + | | | |-S4 | - // + | | 2p | 2p | | - // - | | S2 | S2 | | - // ======================= - // S5 S6 S8 S9 - // - for i in 0..5 { - let zero = u32x8::splat(0).into_bits(); - let S1: u32x8 = _mm256_permutevar8x32_epi32(t1.0[i].into_bits(), c0).into_bits(); - let S2: u32x8 = _mm256_permutevar8x32_epi32(t1.0[i].into_bits(), c1).into_bits(); - let S3_2: u32x8 = _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into_bits(), 0b01010000).into_bits(); - // tmp0 = (0 0 2*S3 -S4) - let tmp0: u32x8 = _mm256_blend_epi32(S3_2.into_bits(), t1.0[i].into_bits(), 0b10100000).into_bits(); - t0.0[i] = (avx2::constants::P_TIMES_2_MASKED.0[i] + tmp0) + S1; - let S2_pos: u32x8 = _mm256_blend_epi32(zero, S2.into_bits(), 0b10100101).into_bits(); - let S2_neg: u32x8 = _mm256_blend_epi32(S2.into_bits(), zero, 0b10100101).into_bits(); - t0.0[i] = t0.0[i] + S2_pos; - t0.0[i] = t0.0[i] - S2_neg; - } - - let c0 = u32x8::new(4,0,6,2,4,0,6,2).into_bits(); // (ABCD) -> (CACA) - let c1 = u32x8::new(5,1,7,3,1,5,3,7).into_bits(); // (ABCD) -> (DBBD) - - for i in 0..5 { - let tmp = t0.0[i]; - t0.0[i] = _mm256_permutevar8x32_epi32(tmp.into_bits(), c0).into_bits(); - t1.0[i] = _mm256_permutevar8x32_epi32(tmp.into_bits(), c1).into_bits(); - } - - ExtendedPoint(&t0 * &t1) - } + // Bounds on (tmp0, tmp1) are (2.33, 1.6) < (2.5, 1.75). + ExtendedPoint(&tmp0 * &tmp1) } pub fn mul_by_pow_2(&self, k: u32) -> ExtendedPoint { @@ -178,6 +167,15 @@ impl ExtendedPoint { } /// A cached point with some precomputed variables used for readdition. +/// +/// # Warning +/// +/// It is not safe to negate this point more than once. +/// +/// # Invariant +/// +/// As long as the `CachedPoint` is not repeatedly negated, its +/// coefficients will be bounded with \\( b < 1.0 \\). #[derive(Copy, Clone, Debug)] pub struct CachedPoint(pub(super) FieldElement32x4); @@ -185,15 +183,16 @@ impl From for CachedPoint { fn from(P: ExtendedPoint) -> CachedPoint { let mut x = P.0; - // x = (S2 S3 Z2 T2) - x.diff_sum(Lanes::AB); + x = x.blend(x.diff_sum(), Lanes::AB); + // x = (X1 - Y1, X2 + Y2, Z2, T2) = (S2 S3 Z2 T2) + x = x * (121666, 121666, 2*121666, 2*121665); // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) - x.scale_by_curve_constants(); + x = x.blend(-x, Lanes::D); // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) - x.negate_D(); + // The coefficients of the output are bounded with b < 0.007. CachedPoint(x) } } @@ -206,13 +205,7 @@ impl Default for CachedPoint { impl Identity for CachedPoint { fn identity() -> CachedPoint { - CachedPoint(FieldElement32x4([ - u32x8::new(121647, 121666, 0, 0, 243332, 67108845, 0, 33554431), - u32x8::new(67108864, 0, 33554431, 0, 0, 67108863, 0, 33554431), - u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), - u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), - u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), - ])) + constants::CACHEDPOINT_IDENTITY } } @@ -224,52 +217,52 @@ impl ConditionallyAssignable for CachedPoint { impl<'a> Neg for &'a CachedPoint { type Output = CachedPoint; - + /// Lazily negate the point. + /// + /// # Warning + /// + /// Because this method does not perform a reduction, it is not + /// safe to repeatedly negate a point. fn neg(self) -> CachedPoint { - let mut neg = *self; - neg.0.swap_AB(); - neg.0.negate_D_lazy(); - neg + let swapped = self.0.shuffle(Shuffle::BACD); + CachedPoint(swapped.blend(swapped.negate_lazy(), Lanes::D)) } } impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { type Output = ExtendedPoint; - /// Uses a slight tweak of the parallel unified formulas of HWCD'08 + /// Add an `ExtendedPoint` and a `CachedPoint`. fn add(self, other: &'b CachedPoint) -> ExtendedPoint { - unsafe { - use core::arch::x86_64::_mm256_permutevar8x32_epi32; + // The coefficients of an `ExtendedPoint` are reduced after + // every operation. If the `CachedPoint` was negated, its + // coefficients grow by one bit. So on input, `self` is + // bounded with `b < 0.007` and `other` is bounded with + // `b < 1.0`. + + let mut tmp = self.0; - let mut tmp = self.0; + tmp = tmp.blend(tmp.diff_sum(), Lanes::AB); + // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) with b < 1.6 - // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) - tmp.diff_sum(Lanes::AB); + // (tmp, other) bounded with b < (1.6, 1.0) < (2.5, 1.75). + tmp = &tmp * &other.0; + // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) - // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) - tmp = &tmp * &other.0; + tmp = tmp.shuffle(Shuffle::ABDC); + // tmp = (S8 S9 S11 S10) - // tmp = (S8 S9 S11 S10) - tmp.swap_CD(); + tmp = tmp.diff_sum(); + // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - tmp.diff_sum(Lanes::ALL); + let t0 = tmp.shuffle(Shuffle::ADDA); + // t0 = (S12 S15 S15 S12) + let t1 = tmp.shuffle(Shuffle::CBCB); + // t1 = (S14 S13 S14 S13) - let c0 = u32x8::new(0,5,2,7,5,0,7,2); // (ABCD) -> (ADDA) - let c1 = u32x8::new(4,1,6,3,4,1,6,3); // (ABCD) -> (CBCB) - - // set t0 = (S12 S15 S15 S12) - // set t1 = (S14 S13 S14 S13) - let mut t0 = FieldElement32x4::zero(); - let mut t1 = FieldElement32x4::zero(); - for i in 0..5 { - t0.0[i] = _mm256_permutevar8x32_epi32(tmp.0[i].into_bits(), c0.into_bits()).into_bits(); - t1.0[i] = _mm256_permutevar8x32_epi32(tmp.0[i].into_bits(), c1.into_bits()).into_bits(); - } - - // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) - ExtendedPoint(&t0 * &t1) - } + // All coefficients of t0, t1 are bounded with b < 1.6. + // Return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + ExtendedPoint(&t0 * &t1) } } @@ -278,8 +271,9 @@ impl<'a, 'b> Sub<&'b CachedPoint> for &'a ExtendedPoint { /// Implement subtraction by negating the point and adding. /// - /// Empirically, this seems about the same cost as a custom subtraction impl (maybe because the - /// benefit is cancelled by increased code size?) + /// Empirically, this seems about the same cost as a custom + /// subtraction impl (maybe because the benefit is cancelled by + /// increased code size?) fn sub(self, other: &'b CachedPoint) -> ExtendedPoint { self + &(-other) } @@ -290,7 +284,7 @@ impl<'a> From<&'a edwards::EdwardsPoint> for LookupTable { let P = ExtendedPoint::from(*point); let mut points = [CachedPoint::from(P); 8]; for i in 0..7 { - points[i+1] = (&P + &points[i]).into(); + points[i + 1] = (&P + &points[i]).into(); } LookupTable(points) } @@ -335,23 +329,23 @@ mod test { macro_rules! print_var { ($x:ident) => { println!("{} = {:?}", stringify!($x), $x.to_bytes()); - } + }; } - let S0 = &Y1 - &X1; // R1 - let S1 = &Y1 + &X1; // R3 - let S2 = &Y2 - &X2; // R2 - let S3 = &Y2 + &X2; // R4 + let S0 = &Y1 - &X1; // R1 + let S1 = &Y1 + &X1; // R3 + let S2 = &Y2 - &X2; // R2 + let S3 = &Y2 + &X2; // R4 print_var!(S0); print_var!(S1); print_var!(S2); print_var!(S3); println!(""); - let S4 = &S0 * &S2; // R5 = R1 * R2 - let S5 = &S1 * &S3; // R6 = R3 * R4 - let S6 = &Z1 * &Z2; // R8 - let S7 = &T1 * &T2; // R7 + let S4 = &S0 * &S2; // R5 = R1 * R2 + let S5 = &S1 * &S3; // R6 = R3 * R4 + let S6 = &Z1 * &Z2; // R8 + let S7 = &T1 * &T2; // R7 print_var!(S4); print_var!(S5); print_var!(S6); @@ -362,8 +356,8 @@ mod test { let S9 = &S5 * &FieldElement64([ 121666,0,0,0,0]); // R6 let S10 = &S6 * &FieldElement64([2*121666,0,0,0,0]); // R8 let S11 = &S7 * &(-&FieldElement64([2*121665,0,0,0,0])); // R7 - print_var!(S8 ); - print_var!(S9 ); + print_var!(S8); + print_var!(S9); print_var!(S10); print_var!(S11); println!(""); @@ -378,12 +372,17 @@ mod test { print_var!(S15); println!(""); - let X3 = &S12 * &S14; // R1 * R2 - let Y3 = &S15 * &S13; // R3 * R4 - let Z3 = &S15 * &S14; // R2 * R3 - let T3 = &S12 * &S13; // R1 * R4 + let X3 = &S12 * &S14; // R1 * R2 + let Y3 = &S15 * &S13; // R3 * R4 + let Z3 = &S15 * &S14; // R2 * R3 + let T3 = &S12 * &S13; // R1 * R4 - edwards::EdwardsPoint{X: X3, Y: Y3, Z: Z3, T: T3} + edwards::EdwardsPoint { + X: X3, + Y: Y3, + Z: Z3, + T: T3, + } } fn addition_test_helper(P: edwards::EdwardsPoint, Q: edwards::EdwardsPoint) { @@ -437,15 +436,15 @@ mod test { } fn serial_double(P: edwards::EdwardsPoint) -> edwards::EdwardsPoint { - let (X1, Y1, Z1, T1) = (P.X, P.Y, P.Z, P.T); + let (X1, Y1, Z1, _T1) = (P.X, P.Y, P.Z, P.T); macro_rules! print_var { ($x:ident) => { println!("{} = {:?}", stringify!($x), $x.to_bytes()); - } + }; } - let S0 = &X1 + &Y1; // R1 + let S0 = &X1 + &Y1; // R1 print_var!(S0); println!(""); @@ -476,7 +475,12 @@ mod test { let Z3 = &S8 * &S6; let T3 = &S5 * &S9; - edwards::EdwardsPoint{X: X3, Y: Y3, Z: Z3, T: T3} + edwards::EdwardsPoint { + X: X3, + Y: Y3, + Z: Z3, + T: T3, + } } fn doubling_test_helper(P: edwards::EdwardsPoint) { diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 96b4baf..175ff14 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -8,54 +8,143 @@ // - Isis Agora Lovecruft // - Henry de Valence -//! 4-way vectorized 32bit field arithmetic using AVX2. +//! An implementation of 4-way vectorized 32bit field arithmetic using +//! AVX2. +//! +//! The `FieldElement32x4` struct provides a vector of four field +//! elements, implemented using AVX2 operations. Its API is designed +//! to abstract away the platform-dependent details, so that point +//! arithmetic can be implemented only in terms of a vector of field +//! elements. +//! +//! At this level, the API is optimized for speed and not safety. The +//! `FieldElement32x4` does not always perform reductions. The pre- +//! and post-conditions on the bounds of the coefficients are +//! documented for each method, but it is the caller's responsibility +//! to ensure that there are no overflows. -#![allow(bad_style)] +#![allow(non_snake_case)] -pub const A_LANES: u8 = 0b0000_0101; -pub const B_LANES: u8 = 0b0000_1010; -pub const C_LANES: u8 = 0b0101_0000; -pub const D_LANES: u8 = 0b1010_0000; +const A_LANES: u8 = 0b0000_0101; +const B_LANES: u8 = 0b0000_1010; +const C_LANES: u8 = 0b0101_0000; +const D_LANES: u8 = 0b1010_0000; -pub const A_LANES64: u8 = 0b00_00_00_11; -pub const B_LANES64: u8 = 0b00_00_11_00; -pub const C_LANES64: u8 = 0b00_11_00_00; -pub const D_LANES64: u8 = 0b11_00_00_00; +#[allow(unused)] +const A_LANES64: u8 = 0b00_00_00_11; +#[allow(unused)] +const B_LANES64: u8 = 0b00_00_11_00; +#[allow(unused)] +const C_LANES64: u8 = 0b00_11_00_00; +#[allow(unused)] +const D_LANES64: u8 = 0b11_00_00_00; -pub const ALL_LANES: u8 = A_LANES | B_LANES | C_LANES | D_LANES; - -use core::ops::Mul; -use core::simd::{IntoBits, u32x8, i32x8, u64x4}; +use core::ops::{Add, Mul, Neg}; +use core::simd::{i32x8, u32x8, u64x4, IntoBits}; +use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO}; use backend::u64::field::FieldElement64; -use backend::avx2::constants::{P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_16_LO, P_TIMES_16_HI}; -#[derive(Copy, Clone)] -pub enum Lanes { - AB, - CD, - ALL, +/// Unpack 32-bit lanes into 64-bit lanes: +/// ``` +/// (a0, b0, a1, b1, c0, d0, c1, d1) +/// ``` +/// into +/// ``` +/// (a0, 0, b0, 0, c0, 0, d0, 0) +/// (a1, 0, b1, 0, c1, 0, d1, 0) +/// ``` +#[inline(always)] +fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { + let a: u32x8; + let b: u32x8; + let zero = i32x8::new(0, 0, 0, 0, 0, 0, 0, 0); + unsafe { + use core::arch::x86_64::_mm256_unpackhi_epi32; + use core::arch::x86_64::_mm256_unpacklo_epi32; + a = _mm256_unpacklo_epi32(src.into_bits(), zero.into_bits()).into_bits(); + b = _mm256_unpackhi_epi32(src.into_bits(), zero.into_bits()).into_bits(); + } + (a, b) } +/// Repack 64-bit lanes into 32-bit lanes: +/// ``` +/// (a0, 0, b0, 0, c0, 0, d0, 0) +/// (a1, 0, b1, 0, c1, 0, d1, 0) +/// ``` +/// into +/// ``` +/// (a0, b0, a1, b1, c0, d0, c1, d1) +/// ``` #[inline(always)] -fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { +fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { unsafe { use core::arch::x86_64::_mm256_blend_epi32; + use core::arch::x86_64::_mm256_shuffle_epi32; - match control { - Lanes::AB => _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32).into_bits(), - Lanes::CD => _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32).into_bits(), - Lanes::ALL => _mm256_blend_epi32(x.into_bits(), y.into_bits(), ALL_LANES as i32).into_bits(), - } + // Input: x = (a0, 0, b0, 0, c0, 0, d0, 0) + // Input: y = (a1, 0, b1, 0, c1, 0, d1, 0) + + let x_shuffled = _mm256_shuffle_epi32(x.into_bits(), 0b11_01_10_00); + let y_shuffled = _mm256_shuffle_epi32(y.into_bits(), 0b10_00_11_01); + + // x' = (a0, b0, 0, 0, c0, d0, 0, 0) + // y' = ( 0, 0, a1, b1, 0, 0, c1, d1) + + return _mm256_blend_epi32(x_shuffled, y_shuffled, 0b11001100).into_bits(); } } -/// A vector of four `FieldElements`, implemented using AVX2. -#[derive(Clone, Copy, Debug)] -pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]); +/// The `Lanes` enum represents a subset of the lanes `A,B,C,D` of a +/// `FieldElement32x4`. +/// +/// It's used to specify blend operations without +/// having to know details about the data layout of the +/// `FieldElement32x4`. +#[derive(Copy, Clone, Debug)] +pub enum Lanes { + C, + D, + AB, + AC, + CD, + AD, + BC, + ABCD, +} + +/// The `Shuffle` enum represents a shuffle of a `FieldElement32x4`. +/// +/// The enum variants are named by what they do to a vector \\( +/// (A,B,C,D) \\); for instance, `Shuffle::BADC` turns \\( (A, B, C, +/// D) \\) into \\( (B, A, D, C) \\). +#[derive(Copy, Clone, Debug)] +pub enum Shuffle { + AAAA, + BBBB, + CACA, + DBBD, + ADDA, + CBCB, + ABAB, + BADC, + BACD, + ABDC, +} + +/// A vector of four field elements. +/// +/// Each operation on a `FieldElement32x4` has documented effects on +/// the bounds of the coefficients. This API is designed for speed +/// and not safety; it is the caller's responsibility to ensure that +/// the post-conditions of one operation are compatible with the +/// pre-conditions of the next. +#[derive(Clone, Copy, Debug)] +pub struct FieldElement32x4(pub(crate) [u32x8; 5]); -use subtle::ConditionallyAssignable; use subtle::Choice; +use subtle::ConditionallyAssignable; impl ConditionallyAssignable for FieldElement32x4 { fn conditional_assign(&mut self, other: &FieldElement32x4, choice: Choice) { @@ -68,10 +157,11 @@ impl ConditionallyAssignable for FieldElement32x4 { } impl FieldElement32x4 { - pub(crate) fn split(&self) -> [FieldElement64; 4] { + /// Split this vector into an array of four (serial) field + /// elements. + pub fn split(&self) -> [FieldElement64; 4] { let mut out = [FieldElement64::zero(); 4]; for i in 0..5 { - let a_2i = self.0[i].extract(0) as u64; // let b_2i = self.0[i].extract(1) as u64; // let a_2i_1 = self.0[i].extract(2) as u64; // `. @@ -90,14 +180,138 @@ impl FieldElement32x4 { out } + /// Rearrange the elements of this vector according to `control`. + /// + /// The `control` parameter should be a compile-time constant, so + /// that when this function is inlined, LLVM is able to lower the + /// shuffle using an immediate. + #[inline] + pub fn shuffle(&self, control: Shuffle) -> FieldElement32x4 { + #[inline(always)] + fn shuffle_lanes(x: u32x8, control: Shuffle) -> u32x8 { + unsafe { + use core::arch::x86_64::_mm256_permutevar8x32_epi32; + + let c: u32x8 = match control { + Shuffle::AAAA => u32x8::new(0, 0, 2, 2, 0, 0, 2, 2), + Shuffle::BBBB => u32x8::new(1, 1, 3, 3, 1, 1, 3, 3), + Shuffle::CACA => u32x8::new(4, 0, 6, 2, 4, 0, 6, 2), + Shuffle::DBBD => u32x8::new(5, 1, 7, 3, 1, 5, 3, 7), + Shuffle::ADDA => u32x8::new(0, 5, 2, 7, 5, 0, 7, 2), + Shuffle::CBCB => u32x8::new(4, 1, 6, 3, 4, 1, 6, 3), + Shuffle::ABAB => u32x8::new(0, 1, 2, 3, 0, 1, 2, 3), + Shuffle::BADC => u32x8::new(1, 0, 3, 2, 5, 4, 7, 6), + Shuffle::BACD => u32x8::new(1, 0, 3, 2, 4, 5, 6, 7), + Shuffle::ABDC => u32x8::new(0, 1, 2, 3, 5, 4, 7, 6), + }; + // Note that this gets turned into a generic LLVM + // shuffle-by-constants, which can be lowered to a simpler + // instruction than a generic permute. + _mm256_permutevar8x32_epi32(x.into_bits(), c.into_bits()).into_bits() + } + } + + FieldElement32x4([ + shuffle_lanes(self.0[0], control), + shuffle_lanes(self.0[1], control), + shuffle_lanes(self.0[2], control), + shuffle_lanes(self.0[3], control), + shuffle_lanes(self.0[4], control), + ]) + } + + /// Blend `self` with `other`, taking lanes specified in `control` from `other`. + /// + /// The `control` parameter should be a compile-time constant, so + /// that this function can be inlined and LLVM can lower it to a + /// blend instruction using an immediate. + #[inline] + pub fn blend(&self, other: FieldElement32x4, control: Lanes) -> FieldElement32x4 { + #[inline(always)] + fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { + unsafe { + use core::arch::x86_64::_mm256_blend_epi32; + + // This would be much cleaner if we could factor out the match + // statement on the control. Unfortunately, rustc forgets + // constant-info very quickly, so we can't even write + // ``` + // match control { + // Lanes::C => { + // let imm = C_LANES as i32; + // _mm256_blend_epi32(..., imm) + // ``` + // let alone + // ``` + // let imm = match control { + // Lanes::C => C_LANES as i32, + // } + // _mm256_blend_epi32(..., imm) + // ``` + // even though both of these would be constant-folded by LLVM + // at a lower level (as happens in the shuffle implementation, + // which does not require a shuffle immediate but *is* lowered + // to immediate shuffles anyways). + match control { + Lanes::C => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), C_LANES as i32).into_bits() + } + Lanes::D => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), D_LANES as i32).into_bits() + } + Lanes::AD => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | D_LANES) as i32) + .into_bits() + } + Lanes::AB => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32) + .into_bits() + } + Lanes::AC => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | C_LANES) as i32) + .into_bits() + } + Lanes::CD => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32) + .into_bits() + } + Lanes::BC => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (B_LANES | C_LANES) as i32) + .into_bits() + } + Lanes::ABCD => _mm256_blend_epi32( + x.into_bits(), + y.into_bits(), + (A_LANES | B_LANES | C_LANES | D_LANES) as i32, + ).into_bits(), + } + } + } + + FieldElement32x4([ + blend_lanes(self.0[0], other.0[0], control), + blend_lanes(self.0[1], other.0[1], control), + blend_lanes(self.0[2], other.0[2], control), + blend_lanes(self.0[3], other.0[3], control), + blend_lanes(self.0[4], other.0[4], control), + ]) + } + + /// Construct a vector of zeros. pub fn zero() -> FieldElement32x4 { - FieldElement32x4([u32x8::splat(0);5]) + FieldElement32x4([u32x8::splat(0); 5]) } + /// Convenience wrapper around `new(x,x,x,x)`. pub fn splat(x: &FieldElement64) -> FieldElement32x4 { - FieldElement32x4::new(x,x,x,x) + FieldElement32x4::new(x, x, x, x) } + /// Create a `FieldElement32x4` from four `FieldElement64`s. + /// + /// # Postconditions + /// + /// The resulting `FieldElement32x4` is bounded with \\( b < 0.0002 \\). pub fn new( x0: &FieldElement64, x1: &FieldElement64, @@ -119,169 +333,104 @@ impl FieldElement32x4 { buf[i] = u32x8::new(a_2i, b_2i, a_2i_1, b_2i_1, c_2i, d_2i, c_2i_1, d_2i_1); } - let mut out = FieldElement32x4(buf); - out.reduce32(); - return out; + // We don't know that the original `FieldElement64`s were + // fully reduced, so the odd limbs may exceed 2^25. + // Reduce them to be sure. + FieldElement32x4(buf).reduce() } - /// Negate the \\(D\\) variable of \\((A,B,C,D)\\). + /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), without + /// performing a reduction. /// - /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. - pub fn negate_D_lazy(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - self.0[0] = _mm256_blend_epi32(self.0[0].into_bits(), (P_TIMES_2_LO - self.0[0]).into_bits(), D_LANES as i32).into_bits(); - self.0[1] = _mm256_blend_epi32(self.0[1].into_bits(), (P_TIMES_2_HI - self.0[1]).into_bits(), D_LANES as i32).into_bits(); - self.0[2] = _mm256_blend_epi32(self.0[2].into_bits(), (P_TIMES_2_HI - self.0[2]).into_bits(), D_LANES as i32).into_bits(); - self.0[3] = _mm256_blend_epi32(self.0[3].into_bits(), (P_TIMES_2_HI - self.0[3]).into_bits(), D_LANES as i32).into_bits(); - self.0[4] = _mm256_blend_epi32(self.0[4].into_bits(), (P_TIMES_2_HI - self.0[4]).into_bits(), D_LANES as i32).into_bits(); - } - } - - /// Negate the \\(D\\) variable of \\((A,B,C,D)\\). + /// # Preconditions /// - /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. - pub fn negate_D(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - self.0[0] = _mm256_blend_epi32(self.0[0].into_bits(), (P_TIMES_16_LO - self.0[0]).into_bits(), D_LANES as i32).into_bits(); - self.0[1] = _mm256_blend_epi32(self.0[1].into_bits(), (P_TIMES_16_HI - self.0[1]).into_bits(), D_LANES as i32).into_bits(); - self.0[2] = _mm256_blend_epi32(self.0[2].into_bits(), (P_TIMES_16_HI - self.0[2]).into_bits(), D_LANES as i32).into_bits(); - self.0[3] = _mm256_blend_epi32(self.0[3].into_bits(), (P_TIMES_16_HI - self.0[3]).into_bits(), D_LANES as i32).into_bits(); - self.0[4] = _mm256_blend_epi32(self.0[4].into_bits(), (P_TIMES_16_HI - self.0[4]).into_bits(), D_LANES as i32).into_bits(); - } - self.reduce32(); - } - - /// Given `self = (A,B,C,D)`, set `self = (B,A,C,D)` - pub fn swap_AB(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_shuffle_epi32; - use core::arch::x86_64::_mm256_blend_epi32; - for i in 0..5 { - let swapped = _mm256_shuffle_epi32(self.0[i].into_bits(), 0b10_11_00_01); - self.0[i] = _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b00001111).into_bits(); - } - } - } - - /// Given `self = (A,B,C,D)`, set `self = (A,B,D,C)` - pub fn swap_CD(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_shuffle_epi32; - use core::arch::x86_64::_mm256_blend_epi32; - for i in 0..5 { - let swapped = _mm256_shuffle_epi32(self.0[i].into_bits(), 0b10_11_00_01); - self.0[i] = _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b11110000).into_bits(); - } - } - } - - /// Given `self = (A,B,C,D)`, set `self = (B - A, B + A, D - C, D + C)` according to `mask`. + /// The coefficients of `self` must be bounded with \\( b < 0.999 \\). /// - /// This is `#[inline(always)]` because the `mask` parameter should be an immediate. - #[inline(always)] - pub fn diff_sum(&mut self, control: Lanes) { - unsafe { - use core::arch::x86_64::{_mm256_shuffle_epi32, _mm256_blend_epi32}; - - let shuffle = |v: u32x8| -> u32x8 { - _mm256_shuffle_epi32(v.into_bits(), 0b10_11_00_01).into_bits() - }; - - let x01 = self.0[0]; - let x01_shuf = shuffle(x01); - let v1 = (x01_shuf + P_TIMES_2_LO) - x01; - let v2 = x01_shuf + x01; - let diffsum01 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[0] = blend_lanes(x01, diffsum01, control); - - let x23 = self.0[1]; - let x23_shuf = shuffle(x23); - let v1 = (x23_shuf + P_TIMES_2_HI) - x23; - let v2 = x23_shuf + x23; - let diffsum23 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[1] = blend_lanes(x23, diffsum23, control); - - let x45 = self.0[2]; - let x45_shuf = shuffle(x45); - let v1 = (x45_shuf + P_TIMES_2_HI) - x45; - let v2 = x45_shuf + x45; - let diffsum45 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[2] = blend_lanes(x45, diffsum45, control); - - let x67 = self.0[3]; - let x67_shuf = shuffle(x67); - let v1 = (x67_shuf + P_TIMES_2_HI) - x67; - let v2 = x67_shuf + x67; - let diffsum67 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[3] = blend_lanes(x67, diffsum67, control); - - let x89 = self.0[4]; - let x89_shuf = shuffle(x89); - let v1 = (x89_shuf + P_TIMES_2_HI) - x89; - let v2 = x89_shuf + x89; - let diffsum89 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[4] = blend_lanes(x89, diffsum89, control); - } - } - - /// Let `self` \\(= (A, B, C, D) \\). + /// # Postconditions /// - /// Compute - /// $$( 121666A, 121666B, 2\cdot 121666C, 2\cdot 121665 D).$$ - pub fn scale_by_curve_constants(&mut self) { - let mut b = [u64x4::splat(0); 10]; - - let consts = u32x8::new(121666, 0, 121666, 0, 2*121666, 0, 2*121665, 0); - - unsafe { - use core::arch::x86_64::_mm256_mul_epu32; - - let (b0, b1) = unpack_pair(self.0[0]); - b[0] = _mm256_mul_epu32(b0.into_bits(), consts.into_bits()).into_bits(); - b[1] = _mm256_mul_epu32(b1.into_bits(), consts.into_bits()).into_bits(); - - let (b2, b3) = unpack_pair(self.0[1]); - b[2] = _mm256_mul_epu32(b2.into_bits(), consts.into_bits()).into_bits(); - b[3] = _mm256_mul_epu32(b3.into_bits(), consts.into_bits()).into_bits(); - - let (b4, b5) = unpack_pair(self.0[2]); - b[4] = _mm256_mul_epu32(b4.into_bits(), consts.into_bits()).into_bits(); - b[5] = _mm256_mul_epu32(b5.into_bits(), consts.into_bits()).into_bits(); - - let (b6, b7) = unpack_pair(self.0[3]); - b[6] = _mm256_mul_epu32(b6.into_bits(), consts.into_bits()).into_bits(); - b[7] = _mm256_mul_epu32(b7.into_bits(), consts.into_bits()).into_bits(); - - let (b8, b9) = unpack_pair(self.0[4]); - b[8] = _mm256_mul_epu32(b8.into_bits(), consts.into_bits()).into_bits(); - b[9] = _mm256_mul_epu32(b9.into_bits(), consts.into_bits()).into_bits(); - } - - *self = FieldElement32x4::reduce64(b); + /// The coefficients of the result are bounded with \\( b < 1 \\). + #[inline] + pub fn negate_lazy(&self) -> FieldElement32x4 { + // The limbs of self are bounded with b < 0.999, while the + // smallest limb of 2*p is 67108845 > 2^{26+0.9999}, so + // underflows are not possible. + FieldElement32x4([ + P_TIMES_2_LO - self.0[0], + P_TIMES_2_HI - self.0[1], + P_TIMES_2_HI - self.0[2], + P_TIMES_2_HI - self.0[3], + P_TIMES_2_HI - self.0[4], + ]) } - pub fn reduce32(&mut self) { + /// Given `self = (A,B,C,D)`, compute `(B - A, B + A, D - C, D + C)`. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 0.01 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 1.6 \\). + #[inline] + pub fn diff_sum(&self) -> FieldElement32x4 { + // tmp1 = (B, A, D, C) + let tmp1 = self.shuffle(Shuffle::BADC); + // tmp2 = (-A, B, -C, D) + let tmp2 = self.blend(self.negate_lazy(), Lanes::AC); + // (B - A, B + A, D - C, D + C) bounded with b < 1.6 + tmp1 + tmp2 + } - let shifts = i32x8::new(26,26,25,25,26,26,25,25); - let masks = u32x8::new((1<<26)-1, (1<<26)-1, (1<<25)-1, (1<<25)-1, - (1<<26)-1, (1<<26)-1, (1<<25)-1, (1<<25)-1); + /// Reduce this vector of field elements \\(\mathrm{mod} p\\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.0002 \\). + #[inline] + pub fn reduce(&self) -> FieldElement32x4 { + let shifts = i32x8::new(26, 26, 25, 25, 26, 26, 25, 25); + let masks = u32x8::new( + (1 << 26) - 1, + (1 << 26) - 1, + (1 << 25) - 1, + (1 << 25) - 1, + (1 << 26) - 1, + (1 << 26) - 1, + (1 << 25) - 1, + (1 << 25) - 1, + ); - let carry = |v: u32x8| -> u32x8 { + // Let c(x) denote the carryout of the coefficient x. + // + // Given ( x0, y0, x1, y1, z0, w0, z1, w1), + // compute (c(x1), c(y1), c(x0), c(y0), c(z1), c(w1), c(z0), c(w0)). + // + // The carryouts are bounded by 2^(32 - 25) = 2^7. + let rotated_carryout = |v: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_srlv_epi32; - _mm256_srlv_epi32(v.into_bits(), shifts.into_bits()).into_bits() - } - }; - - let swap_lanes = |v: u32x8| -> u32x8 { - unsafe { use core::arch::x86_64::_mm256_shuffle_epi32; - _mm256_shuffle_epi32(v.into_bits(), 0b01_00_11_10).into_bits() + + let c = _mm256_srlv_epi32(v.into_bits(), shifts.into_bits()); + _mm256_shuffle_epi32(c, 0b01_00_11_10).into_bits() } }; + // Combine (lo, lo, lo, lo, lo, lo, lo, lo) + // with (hi, hi, hi, hi, hi, hi, hi, hi) + // to (lo, lo, hi, hi, lo, lo, hi, hi) + // + // This allows combining carryouts, e.g., + // + // lo (c(x1), c(y1), c(x0), c(y0), c(z1), c(w1), c(z0), c(w0)) + // hi (c(x3), c(y3), c(x2), c(y2), c(z3), c(w3), c(z2), c(w2)) + // -> (c(x1), c(y1), c(x2), c(y2), c(z1), c(w1), c(z2), c(w2)) + // + // which is exactly the vector of carryins for + // + // ( x2, y2, x3, y3, z2, w2, z3, w3). + // let combine = |v_lo: u32x8, v_hi: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_blend_epi32; @@ -289,49 +438,79 @@ impl FieldElement32x4 { } }; - let v = &mut self.0; + let mut v = self.0; - let c10 = swap_lanes(carry(v[0])); + let c10 = rotated_carryout(v[0]); v[0] = (v[0] & masks) + combine(u32x8::splat(0), c10); - let c32 = swap_lanes(carry(v[1])); + + let c32 = rotated_carryout(v[1]); v[1] = (v[1] & masks) + combine(c10, c32); - let c54 = swap_lanes(carry(v[2])); + + let c54 = rotated_carryout(v[2]); v[2] = (v[2] & masks) + combine(c32, c54); - let c76 = swap_lanes(carry(v[3])); + + let c76 = rotated_carryout(v[3]); v[3] = (v[3] & masks) + combine(c54, c76); - let c98 = swap_lanes(carry(v[4])); + + let c98 = rotated_carryout(v[4]); v[4] = (v[4] & masks) + combine(c76, c98); - // Still need to account for c9 - // c98 = (c9, c9, c8, c8, c9, c9, c8, c8) - // - let c9_19: u32x8; - unsafe { + let c9_19: u32x8 = unsafe { use core::arch::x86_64::_mm256_mul_epu32; use core::arch::x86_64::_mm256_shuffle_epi32; - let c9_spread = _mm256_shuffle_epi32(c98.into_bits(), 0b11_01_10_00); - let c9_19_spread = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into_bits()); - c9_19 = _mm256_shuffle_epi32(c9_19_spread, 0b11_01_10_00).into_bits(); - } + // Need to rearrange c98, since vpmuludq uses the low + // 32-bits of each 64-bit lane to compute the product: + // + // c98 = (c(x9), c(y9), c(x8), c(y8), c(z9), c(w9), c(z8), c(w8)); + // c9_spread = (c(x9), c(x8), c(y9), c(y8), c(z9), c(z8), c(w9), c(w8)). + let c9_spread = _mm256_shuffle_epi32(c98.into_bits(), 0b11_01_10_00); + + // Since the carryouts are bounded by 2^7, their products with 19 + // are bounded by 2^11.25. This means that + // + // c9_19_spread = (19*c(x9), 0, 19*c(y9), 0, 19*c(z9), 0, 19*c(w9), 0). + let c9_19_spread = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into_bits()); + + // Unshuffle: + // c9_19 = (19*c(x9), 19*c(y9), 0, 0, 19*c(z9), 19*c(w9), 0, 0). + _mm256_shuffle_epi32(c9_19_spread, 0b11_01_10_00).into_bits() + }; + + // Add the final carryin. v[0] = v[0] + c9_19; + + // Each output coefficient has exactly one carryin, which is + // bounded by 2^11.25, so they are bounded as + // + // c_even < 2^26 + 2^11.25 < 26.00006 < 2^{26+b} + // c_odd < 2^25 + 2^11.25 < 25.0001 < 2^{25+b} + // + // where b = 0.0002. + FieldElement32x4(v) } - pub fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { + /// Given an array of wide coefficients, reduce them to a `FieldElement32x4`. + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + #[inline] + fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { // These aren't const because splat isn't a const fn - let LOW_25_BITS: u64x4 = u64x4::splat((1<<25)-1); - let LOW_26_BITS: u64x4 = u64x4::splat((1<<26)-1); + let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1); + let LOW_26_BITS: u64x4 = u64x4::splat((1 << 26) - 1); // Carry the value from limb i = 0..8 to limb i+1 let carry = |z: &mut [u64x4; 10], i: usize| { debug_assert!(i < 9); if i % 2 == 0 { // Even limbs have 26 bits - z[i+1] = z[i+1] + (z[i] >> 26); + z[i + 1] = z[i + 1] + (z[i] >> 26); z[i] = z[i] & LOW_26_BITS; } else { // Odd limbs have 25 bits - z[i+1] = z[i+1] + (z[i] >> 25); + z[i + 1] = z[i + 1] + (z[i] >> 25); z[i] = z[i] & LOW_25_BITS; } }; @@ -370,8 +549,13 @@ impl FieldElement32x4 { z[1] = z[1] + c1; // z1 < 2^25 + 2^17.25 < 2^25.0067 carry(&mut z, 0); // z0 < 2^26, z1 < 2^25.0067 + 2^4.33 = 2^25.007 - // Now repack the [u64x4; 10] into a FieldElement32x4 - + // The output coefficients are bounded with + // + // b = 0.007 for z[1] + // b = 0.0004 for z[5] + // b = 0 for other z[i]. + // + // So the packed result is bounded with b = 0.007. FieldElement32x4([ repack_pair(z[0].into_bits(), z[1].into_bits()), repack_pair(z[2].into_bits(), z[3].into_bits()), @@ -380,65 +564,30 @@ impl FieldElement32x4 { repack_pair(z[8].into_bits(), z[9].into_bits()), ]) } -} -#[inline(always)] -pub fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { - let a: u32x8; - let b: u32x8; - let zero = i32x8::new(0,0,0,0,0,0,0,0); - unsafe { - use core::arch::x86_64::_mm256_unpackhi_epi32; - use core::arch::x86_64::_mm256_unpacklo_epi32; - a = _mm256_unpacklo_epi32(src.into_bits(), zero.into_bits()).into_bits(); - b = _mm256_unpackhi_epi32(src.into_bits(), zero.into_bits()).into_bits(); - } - (a,b) -} - -#[inline(always)] -pub fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { - unsafe { - use core::arch::x86_64::_mm256_shuffle_epi32; - use core::arch::x86_64::_mm256_blend_epi32; - - // Input: x = (a0, 0, b0, 0, c0, 0, d0) - // Input: y = (a1, 0, b1, 0, c1, 0, d1) - - let x_shuffled = _mm256_shuffle_epi32(x.into_bits(), 0b11_01_10_00); - let y_shuffled = _mm256_shuffle_epi32(y.into_bits(), 0b10_00_11_01); - - // x' = (a0, b0, 0, 0, c0, d0, 0, 0) - // y' = ( 0, 0, a1, b1, 0, 0, c1, d1) - - return _mm256_blend_epi32(x_shuffled, y_shuffled, 0b11001100).into_bits(); - } -} - -impl FieldElement32x4 { - /// Square this field element, then conditionally negate according - /// to `neg_mask`. This parameter is hardcoded as `neg_mask = - /// D_LANES64` to negate the \\( D \\) value. + /// Square this field element, and negate the result's \\(D\\) value. /// - /// # Precondition + /// # Preconditions /// - /// Limbs must be bounded by bit-excess \\( b < 2.0 \\). + /// The coefficients of `self` must be bounded with \\( b < 1.5 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). pub fn square_and_negate_D(&self) -> FieldElement32x4 { - let neg_mask = D_LANES64; - #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } #[inline(always)] fn m_lo(x: u32x8, y: u32x8) -> u32x8 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } - let v19 = u32x8::new(19,0,19,0,19,0,19,0); + let v19 = u32x8::new(19, 0, 19, 0, 19, 0, 19, 0); let (x0, x1) = unpack_pair(self.0[0]); let (x2, x3) = unpack_pair(self.0[1]); @@ -514,19 +663,109 @@ impl FieldElement32x4 { } } +impl Neg for FieldElement32x4 { + type Output = FieldElement32x4; + + /// Negate this field element, performing a reduction. + /// + /// If the coefficients are known to be small, use `negate_lazy` + /// to avoid performing a reduction. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 4.0 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.0002 \\). + #[inline] + fn neg(self) -> FieldElement32x4 { + FieldElement32x4([ + P_TIMES_16_LO - self.0[0], + P_TIMES_16_HI - self.0[1], + P_TIMES_16_HI - self.0[2], + P_TIMES_16_HI - self.0[3], + P_TIMES_16_HI - self.0[4], + ]).reduce() + } +} + +impl Add for FieldElement32x4 { + type Output = FieldElement32x4; + /// Add two `FieldElement32x4`s, without performing a reduction. + #[inline] + fn add(self, rhs: FieldElement32x4) -> FieldElement32x4 { + FieldElement32x4([ + self.0[0] + rhs.0[0], + self.0[1] + rhs.0[1], + self.0[2] + rhs.0[2], + self.0[3] + rhs.0[3], + self.0[4] + rhs.0[4], + ]) + } +} + +impl Mul<(u32, u32, u32, u32)> for FieldElement32x4 { + type Output = FieldElement32x4; + /// Perform a multiplication by a vector of small constants. + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + #[inline] + fn mul(self, scalars: (u32, u32, u32, u32)) -> FieldElement32x4 { + unsafe { + use core::arch::x86_64::_mm256_mul_epu32; + + let consts = u32x8::new(scalars.0, 0, scalars.1, 0, scalars.2, 0, scalars.3, 0); + + let (b0, b1) = unpack_pair(self.0[0]); + let (b2, b3) = unpack_pair(self.0[1]); + let (b4, b5) = unpack_pair(self.0[2]); + let (b6, b7) = unpack_pair(self.0[3]); + let (b8, b9) = unpack_pair(self.0[4]); + + FieldElement32x4::reduce64([ + _mm256_mul_epu32(b0.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b1.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b2.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b3.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b4.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b5.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b6.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b7.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b8.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b9.into_bits(), consts.into_bits()).into_bits(), + ]) + } + } +} + impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { type Output = FieldElement32x4; - fn mul(self, _rhs: &'b FieldElement32x4) -> FieldElement32x4 { + /// Multiply `self` by `rhs`. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 2.5 \\). + /// + /// The coefficients of `rhs` must be bounded with \\( b < 1.75 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + /// + fn mul(self, rhs: &'b FieldElement32x4) -> FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } #[inline(always)] fn m_lo(x: u32x8, y: u32x8) -> u32x8 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } let (x0, x1) = unpack_pair(self.0[0]); @@ -535,21 +774,21 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let (x6, x7) = unpack_pair(self.0[3]); let (x8, x9) = unpack_pair(self.0[4]); - let (y0, y1) = unpack_pair(_rhs.0[0]); - let (y2, y3) = unpack_pair(_rhs.0[1]); - let (y4, y5) = unpack_pair(_rhs.0[2]); - let (y6, y7) = unpack_pair(_rhs.0[3]); - let (y8, y9) = unpack_pair(_rhs.0[4]); + let (y0, y1) = unpack_pair(rhs.0[0]); + let (y2, y3) = unpack_pair(rhs.0[1]); + let (y4, y5) = unpack_pair(rhs.0[2]); + let (y6, y7) = unpack_pair(rhs.0[3]); + let (y8, y9) = unpack_pair(rhs.0[4]); - let v19 = u32x8::new(19,0,19,0,19,0,19,0); + let v19 = u32x8::new(19, 0, 19, 0, 19, 0, 19, 0); let y1_19 = m_lo(v19, y1); // This fits in a u32 let y2_19 = m_lo(v19, y2); // iff 26 + b + lg(19) < 32 let y3_19 = m_lo(v19, y3); // if b < 32 - 26 - 4.248 = 1.752 let y4_19 = m_lo(v19, y4); - let y5_19 = m_lo(v19, y5); // below, b<2.5: this is a bottleneck, - let y6_19 = m_lo(v19, y6); // could be avoided by promoting to - let y7_19 = m_lo(v19, y7); // u64 here instead of in m() + let y5_19 = m_lo(v19, y5); + let y6_19 = m_lo(v19, y6); + let y7_19 = m_lo(v19, y7); let y8_19 = m_lo(v19, y8); let y9_19 = m_lo(v19, y9); @@ -570,6 +809,44 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let z8 = m(x0,y8) + m(x1_2,y7) + m(x2,y6) + m(x3_2,y5) + m(x4,y4) + m(x5_2,y3) + m(x6,y2) + m(x7_2,y1) + m(x8,y0) + m(x9_2,y9_19); let z9 = m(x0,y9) + m(x1,y8) + m(x2,y7) + m(x3,y6) + m(x4,y5) + m(x5,y4) + m(x6,y3) + m(x7,y2) + m(x8,y1) + m(x9,y0); + // The bounds on z[i] are the same as in the serial 32-bit code + // and the comment below is copied from there: + + // How big is the contribution to z[i+j] from x[i], y[j]? + // + // Using the bounds above, we get: + // + // i even, j even: x[i]*y[j] < 2^(26+b)*2^(26+b) = 2*2^(51+2*b) + // i odd, j even: x[i]*y[j] < 2^(25+b)*2^(26+b) = 1*2^(51+2*b) + // i even, j odd: x[i]*y[j] < 2^(26+b)*2^(25+b) = 1*2^(51+2*b) + // i odd, j odd: 2*x[i]*y[j] < 2*2^(25+b)*2^(25+b) = 1*2^(51+2*b) + // + // We perform inline reduction mod p by replacing 2^255 by 19 + // (since 2^255 - 19 = 0 mod p). This adds a factor of 19, so + // we get the bounds (z0 is the biggest one, but calculated for + // posterity here in case finer estimation is needed later): + // + // z0 < ( 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 249*2^(51 + 2*b) + // z1 < ( 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 154*2^(51 + 2*b) + // z2 < ( 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 195*2^(51 + 2*b) + // z3 < ( 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 118*2^(51 + 2*b) + // z4 < ( 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 141*2^(51 + 2*b) + // z5 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 82*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 87*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 )*2^(51 + 2b) = 46*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 )*2^(51 + 2b) = 33*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 )*2^(51 + 2b) = 10*2^(51 + 2*b) + // + // So z[0] fits into a u64 if 51 + 2*b + lg(249) < 64 + // if b < 2.5. + + // In fact this bound is slightly sloppy, since it treats both + // inputs x and y as being bounded by the same parameter b, + // while they are in fact bounded by b_x and b_y, and we + // already require that b_y < 1.75 in order to fit the + // multiplications by 19 into a u32. The tighter bound on b_y + // means we could get a tighter bound on the outputs, or a + // looser bound on b_x. FieldElement32x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9]) } } @@ -582,13 +859,14 @@ mod test { #[test] fn scale_by_curve_constants() { let mut x = FieldElement32x4::splat(&FieldElement64::one()); - x.scale_by_curve_constants(); + + x = x * (121666, 121666, 2*121666, 2*121665); let xs = x.split(); - assert_eq!(xs[0], FieldElement64([ 121666,0,0,0,0])); - assert_eq!(xs[1], FieldElement64([ 121666,0,0,0,0])); - assert_eq!(xs[2], FieldElement64([2*121666,0,0,0,0])); - assert_eq!(xs[3], FieldElement64([2*121665,0,0,0,0])); + assert_eq!(xs[0], FieldElement64([121666, 0, 0, 0, 0])); + assert_eq!(xs[1], FieldElement64([121666, 0, 0, 0, 0])); + assert_eq!(xs[2], FieldElement64([2 * 121666, 0, 0, 0, 0])); + assert_eq!(xs[3], FieldElement64([2 * 121665, 0, 0, 0, 0])); } #[test] @@ -598,8 +876,7 @@ mod test { let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]); let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); - let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); - vec.diff_sum(Lanes::ALL); + let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3).diff_sum(); let result = vec.split(); @@ -607,16 +884,6 @@ mod test { assert_eq!(result[1], &x1 + &x0); assert_eq!(result[2], &x3 - &x2); assert_eq!(result[3], &x3 + &x2); - - let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); - vec.diff_sum(Lanes::AB); // leave C,D unchanged - - let result = vec.split(); - - assert_eq!(result[0], &x1 - &x0); - assert_eq!(result[1], &x1 + &x0); - assert_eq!(result[2], x2); - assert_eq!(result[3], x3); } #[test] @@ -636,7 +903,6 @@ mod test { assert_eq!(result[3], -&(&x3 * &x3)); } - #[test] fn multiply_vs_serial() { let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]); @@ -666,7 +932,7 @@ mod test { let src = vec.0[0]; - let (a,b) = unpack_pair(src); + let (a, b) = unpack_pair(src); let expected_a = u32x8::new(10000, 0, 10100, 0, 10200, 0, 10300, 0); let expected_b = u32x8::new(10001, 0, 10101, 0, 10201, 0, 10301, 0); @@ -674,7 +940,7 @@ mod test { assert_eq!(a, expected_a); assert_eq!(b, expected_b); - let expected_src = repack_pair(a,b); + let expected_src = repack_pair(a, b); assert_eq!(src, expected_src); } diff --git a/src/backend/avx2/scalar_mul/straus.rs b/src/backend/avx2/scalar_mul/straus.rs index 52c3079..053713f 100644 --- a/src/backend/avx2/scalar_mul/straus.rs +++ b/src/backend/avx2/scalar_mul/straus.rs @@ -8,6 +8,8 @@ // - Isis Agora Lovecruft // - Henry de Valence +#![allow(non_snake_case)] + use core::borrow::Borrow; use clear_on_drop::ClearOnDrop; diff --git a/src/backend/mod.rs b/src/backend/mod.rs index b095eef..aa44d52 100644 --- a/src/backend/mod.rs +++ b/src/backend/mod.rs @@ -27,6 +27,6 @@ pub mod u32; #[cfg(feature = "u64_backend")] pub mod u64; -#[cfg(all(feature = "avx2_backend", feature = "yolocrypto", target_feature = "avx2"))] +#[cfg(all(feature = "avx2_backend", target_feature = "avx2"))] pub mod avx2; diff --git a/src/lib.rs b/src/lib.rs index 4ebffe7..70f80f5 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -14,7 +14,7 @@ #![cfg_attr(feature = "nightly", feature(cfg_target_feature))] #![cfg_attr(feature = "nightly", feature(external_doc))] -#![cfg_attr(all(feature = "nightly", feature = "yolocrypto"), feature(stdsimd))] +#![cfg_attr(all(feature = "nightly", feature = "avx2_backend"), feature(stdsimd))] // Refuse to compile if documentation is missing, but only on nightly. //