saymrwulf/curve25519-dalek-source
1
0
Fork 0
mirror of https://github.com/saymrwulf/curve25519-dalek-source.git synced 2026-09-10 21:21:13 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'release/0.18.0'

Browse source
This commit is contained in:
Henry de Valence 2018-06-22 12:27:14 -07:00
commit 3a5aef7eea
12 changed files with 1036 additions and 756 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
.travis.yml
View file

@ -10,7 +10,7 @@ env:
# Tests the u64 backend # Tests the u64 backend
- TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std u64_backend' - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std u64_backend'
# Tests the avx2 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 # Tests serde support and default feature selection
- TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='serde' - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='serde'
# Tests building without std. We have to select a backend, so we select the one # Tests building without std. We have to select a backend, so we select the one
@ -21,7 +21,7 @@ matrix:
exclude: exclude:
# Test the avx2 backend only on nightly # Test the avx2 backend only on nightly
- rust: stable - 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. # Test no_std only on nightly.
- rust: stable - rust: stable
env: TEST_COMMAND=build EXTRA_FLAGS=--no-default-features FEATURES='u32_backend' env: TEST_COMMAND=build EXTRA_FLAGS=--no-default-features FEATURES='u32_backend'

6
Cargo.toml
View file

@ -1,6 +1,6 @@
[package] [package]
name = "curve25519-dalek" name = "curve25519-dalek"
version = "0.17.0" version = "0.18.0"
authors = ["Isis Lovecruft <isis@patternsinthevoid.net>", authors = ["Isis Lovecruft <isis@patternsinthevoid.net>",
"Henry de Valence <hdevalence@hdevalence.ca>"] "Henry de Valence <hdevalence@hdevalence.ca>"]
readme = "README.md" readme = "README.md"
@ -41,7 +41,7 @@ harness = false
# match exactly, since the build.rs uses the crate itself as a library. # match exactly, since the build.rs uses the crate itself as a library.
[dependencies] [dependencies]
rand = { version = "0.5.0-pre.2", default-features = false } rand = { version = "0.5.0", default-features = false }
byteorder = { version = "1", default-features = false } byteorder = { version = "1", default-features = false }
digest = "0.7" digest = "0.7"
generic-array = "0.9" generic-array = "0.9"
@ -50,7 +50,7 @@ subtle = { version = "0.6", features = ["generic-impls"], default-features = fal
serde = { version = "1.0", optional = true } serde = { version = "1.0", optional = true }
[build-dependencies] [build-dependencies]
rand = { version = "0.5.0-pre.2", default-features = false } rand = { version = "0.5.0", default-features = false }
byteorder = { version = "1", default-features = false } byteorder = { version = "1", default-features = false }
digest = "0.7" digest = "0.7"
generic-array = "0.9" generic-array = "0.9"

2
Makefile
View file

@ -1,4 +1,4 @@
FEATURES := nightly yolocrypto FEATURES := nightly yolocrypto avx2_backend
doc: doc:
cargo rustdoc --features "$(FEATURES)" -- --html-in-header docs/assets/rustdoc-include-katex-header.html cargo rustdoc --features "$(FEATURES)" -- --html-in-header docs/assets/rustdoc-include-katex-header.html

15
README.md
View file

@ -50,7 +50,7 @@ make doc-internal
To import `curve25519-dalek`, add the following to the dependencies section of To import `curve25519-dalek`, add the following to the dependencies section of
your project's `Cargo.toml`: your project's `Cargo.toml`:
```toml ```toml
curve25519-dalek = "^0.17" curve25519-dalek = "^0.18"
``` ```
Then import the crate as: Then import the crate as:
```rust,no_run ```rust,no_run
@ -71,15 +71,17 @@ Curve arithmetic is implemented using one of the following backends:
* a `u32` backend using `u64` products; * a `u32` backend using `u64` products;
* a `u64` backend using `u128` products; * a `u64` backend using `u128` products;
* an experimental AVX2 backend, available using the `yolocrypto` feature when * an `avx2` backend using parallel formulas, available when compiling for a
compiling for a target with `target_feature=+avx2`. target with `target_feature=+avx2`.
By default the `u64` backend is selected. To select a specific backend, use: By default the `u64` backend is selected. To select a specific backend, use:
```sh ```sh
cargo build --no-default-features --features "std u32_backend" 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 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]: Benchmarks are run using [`criterion.rs`][criterion]:
@ -88,7 +90,7 @@ Benchmarks are run using [`criterion.rs`][criterion]:
export RUSTFLAGS="-C target_cpu=native" export RUSTFLAGS="-C target_cpu=native"
cargo bench --no-default-features --features "std u32_backend" 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 u64_backend"
cargo bench --no-default-features --features "std avx2_backend yolocrypto" cargo bench --no-default-features --features "std avx2_backend"
``` ```
# Contributing # Contributing
@ -117,7 +119,8 @@ to the Dalek race.*
Portions of this library were originally a port of [Adam Langley's Portions of this library were originally a port of [Adam Langley's
Golang ed25519 library](https://github.com/agl/ed25519), which was in 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 fast `u32` and `u64` scalar arithmetic was implemented by Andrew Moon, and
the addition chain for scalar inversion was provided by Brian Smith. the addition chain for scalar inversion was provided by Brian Smith.

2
build.rs
View file

@ -1,5 +1,5 @@
#![cfg_attr(feature = "nightly", feature(cfg_target_feature))] #![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(unused_variables)]
#![allow(non_snake_case)] #![allow(non_snake_case)]
#![allow(dead_code)] #![allow(dead_code)]

563
docs/avx2-notes.md
View file

@ -1,42 +1,48 @@
An implementation of group operations on the twisted Edwards form of A vectorized implementation of group operations on the twisted Edwards
Curve25519, using AVX2 to implement the 4-way parallel formulas of form of Curve25519, using a modification of the 4-way parallel
Hisil, Wong, Carter, and Dawson (HWCD). formulas of Hisil, Wong, Carter, and Dawson.
Their 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08], which
introduced the extended coordinates used in other parts of `-dalek`,
also describes 4-way parallel formulas for point addition and
doubling:
* a unified addition algorithm taking an effective \\(2\mathbf M + # Overview
1\mathbf D\\);
* a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf The 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08] by Hisil,
S\\); 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 However, they also describe 4-way parallel formulas for point addition
an effective \\(2 \mathbf M \\). 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 Unlike their serial formulas, which are used widely, their parallel
multiplication and squaring of generic field elements and \\(\mathbf formulas do not seem to have been implemented in software before. The
D\\) represents the cost of multiplication by a curve constant. 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. The reason may be that HWCD08 describe their formulas as operating on
A 2015 paper by Hernández and López mentions using AVX2 for the X25519 four independent processors, which would make a software
Montgomery ladder, but neither the paper nor the code are publicly implementation impractical: all of the operations are too low-latency
available, and it apparently gives only a [slight speedup][avx2trac]. to effectively synchronize. But a closer inspection reveals that the
The 2008 HWCD paper also describes and analyzes a 2-wide variant of the (more expensive) multiplication and squaring steps are uniform, while
Montgomery ladder (for comparison with parallel Edwards formulas); this the instruction divergence occurs in the (much cheaper) addition and
strategy was used in 2015 by Tung Chou's `sandy2x` implementation, which subtraction steps. This means that a SIMD implementation can perform
used a 2-wide field implementation in 128-bit vector registers. the expensive steps uniformly, and handle divergence in the
Curiously, however, although the [`sandy2x` paper][sandy2x] also inexpensive steps using masking.
implements Edwards arithmetic, and cites the HWCD paper, it doesn't
mention the parallel formulas from HWCD, suggesting that they have been
overlooked for software implementations.
The notes below describe a tweak to the \\( 2\mathbf M + 1\mathbf D \\) These notes describe modifications to the original parallel formulas
unified addition formulas to give \\( 2\mathbf M \\) readdition with to allow a SIMD implementation, and this module contains an
\\(1\mathbf D\\) precomputation, and a tweak to the doubling formulas to implementation of the modified formulas using 256-bit AVX2 vector
avoid an extra reduction. These tweaked formulas are the ones used by operations.
the `avx2` backend of `curve25519-dalek`.
# Parallel formulas in HWCD'08 # 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 \\) | | | \\( 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 \\) | | \\(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 Notice that the \\(1\mathbf M\\) and \\(1\mathbf S\\) steps are
low-latency to parallelize across threads. However, the main cost uniform. The non-uniform steps are all inexpensive additions or
is in the multiplication and squaring steps, which are uniform, while subtractions, with the exception of the multiplication by the curve
the divergent steps involve inexpensive additions and subtractions. constant \\(k = 2d\\):
$$
R\_7 \gets 2 d R\_7.
$$
This means we can use SIMD to implement the expensive portions in HWCD suggest parallelising this step by breaking \\(k = 2d\\) into four
parallel, and handle the instruction divergence on the inexpensive parts
using masking.
The remaining obstacle to parallelism is the multiplication by the curve
constant \\(k = 2d\\). In the Curve25519 case, this is
$$ k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. $$
HWCD suggest parallelising this step by breaking \\(k\\) into four
parts as \\(k = k_0 + 2\^n k_1 + 2\^{2n} k_2 + 2\^{3n} k_3 \\) and 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 computing \\(k_i R_7 \\) in parallel. This is quite awkward, but if
somewhat awkward in our case, since we would normally represent the curve constant is a ratio \\( d = d\_1/d\_2 \\), then projective
\\(k\\) as \\( 10 \\) 32-bit limbs, and \\(10 \\) is not divisible coordinates allow us to instead compute
by \\(4\\), so we would need a specialized routine to perform a $$
vectorized multiplication by 64-bit constants. (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 # Modified parallel formulas
\\(R_7\\) by \\( -2\cdot 121665 \\) and multiply the other three
variables by \\(121666\\). This trick was suggested by Mike
Hamburg. Ignoring the sign for the moment, since
\\(2 \cdot 121666 < 2\^{18}\\), all these constants fit in 32 bits,
so (up to sign) this can be done in parallel as four multiplications
by small constants \\( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \\).
How do we handle the sign? Using the modifications sketched above, we can write SIMD-friendly
Since we're primarily interested in Ristretto performance, not versions of the parallel formulas as follows. To avoid confusion with
Curve25519 performance, we could alternately work on the the original formulas, temporary variables are named \\(S\\) instead
\\(4\\)-isogenous "IsoEd25519" curve, which has \\(d = 121665\\). of \\(R\\) and are in static single-assignment form.
However, this would only save the negation step, since multiplying
one field element by a 32-bit constant is not much easier than
multiplying four field elements by 32-bit constants, and it would
prevent accelerating Curve25519, so we don't make this choice.
Instead, we just negate one lane, and move the \\(1 \mathbf D\\)
into precomputation (see below).
# Tweaked formulas
After tweaking the formulas as described above, we obtain the
following. To avoid confusion with the original HWCD formulas,
temporary variables are named \\(S\\) instead of \\(R\\) and are in
static single-assignment form.
## Addition ## Addition
This implementation only implements readdition, but the tweaked addition To add points
formulas are described first. To add points \\(P_1 = (X_1 : Y_1 : Z_1 : \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\)
T_1) \\) and \\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\), we compute and
\\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\),
we compute
$$ $$
\begin{aligned} \begin{aligned}
S\_0 &\gets Y\_1 - X\_1 \\\\ (S\_0 &&,&& S\_1 &&,&& S\_2 &&,&& S\_3 )
S\_1 &\gets Y\_1 + X\_1 \\\\ &\gets
S\_2 &\gets Y\_2 - X\_2 \\\\ (Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Y\_2 - X\_2&&,&& Y\_2 + X\_2)
S\_3 &\gets 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} \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 \\). 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 ## 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} \begin{aligned}
S\_2 &\gets Y\_2 - X\_2 \\\\ (S\_2' &&,&& S\_3' &&,&& Z\_2' &&,&& T\_2' )
S\_3 &\gets Y\_2 + X\_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} \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} \begin{aligned}
S\_2' &\gets S\_2 \cdot 121666 \\\\ (S\_0 &&,&& S\_1 &&,&& Z\_1 &&,&& T\_1 )
S\_3' &\gets S\_3 \cdot 121666 \\\\ &\gets
Z\_2' &\gets Z\_2 \cdot 2 \cdot 121666 \\\\ (Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Z\_1 &&,&& T\_1)
T\_2' &\gets T\_2 \cdot -2 \cdot 121665 \\\\ \\\\
(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} \end{aligned}
$$ $$
which costs only \\( 2\mathbf M \\). This precomputation is
to obtain the `CachedPoint` \\( (S\_2', S\_3', Z\_2', T\_2') \\). essentially similar to the precomputation that HWCD suggest for their
This precomputation is essentially the same as that suggested in serial formulas. Because the cost of precomputation and then
§3.1 of HWCD, with the difference that the multiplication by the curve readdition is the same as addition, it's sufficient to only
constant \\( -121665 / 121666 \\) is spread over all four implement caching and readdition.
coordinates, to allow a vectorized computation of four
multiplications of small constants instead of a serial computation
of multiplication by a large constant.
To perform readdition of \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and
\\(P_2 = (S\_2', S\_3', Z\_2', T\_2') \\), we compute
$$
\begin{aligned}
S\_0 &\gets Y\_1 - X\_1 \\\\
S\_1 &\gets Y\_1 + X\_1
\end{aligned}
$$
$$
\begin{aligned}
S\_8 &\gets S\_0 S\_2' \\\\
S\_9 &\gets S\_1 S\_3' \\\\
S\_{10} &\gets Z\_1 Z\_2' \\\\
S\_{11} &\gets T\_1 T\_2'
\end{aligned}
$$
$$
\begin{aligned}
S\_{12} &\gets S\_9 - S\_8 \\\\
S\_{13} &\gets S\_9 + S\_8 \\\\
S\_{14} &\gets S\_{10} - S\_{11} \\\\
S\_{15} &\gets S\_{10} + S\_{11}
\end{aligned}
$$
$$
\begin{aligned}
X\_3 &\gets S\_{12} S\_{14} \\\\
Y\_3 &\gets S\_{15} S\_{13} \\\\
Z\_3 &\gets S\_{15} S\_{14} \\\\
T\_3 &\gets S\_{12} S\_{13}
\end{aligned}
$$
to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\).
Compared to the addition formulas above, this saves \\( 1\mathbf D \\).
## Doubling ## 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 To double a point \\( P = (X\_1 : Y\_1 : Z\_1 : T\_1) \\), we compute
$$ S\_0 \gets X\_1 + Y\_1 $$
$$ $$
\begin{aligned} \begin{aligned}
S\_1 &\gets X\_1\^2 \\\\ (X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& S\_0)
S\_2 &\gets Y\_1\^2 \\\\ &\gets
S\_3 &\gets Z\_1\^2 \\\\ (X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& X\_1 + Y\_1)
S\_4 &\gets S\_0\^2 \\\\
(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} \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 \\). 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 The intermediate step between the squaring and multiplication requires
are less nice. However, with some careful bounds-juggling, it is a long chain of additions, but with some care and finesse,
possible to implement them without inserting extra carry chains, as described below, it is possible (in our case) to arrange this
described below. 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 # 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 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. 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.
$$ The multiplication routine requires that its inputs are bounded with
\begin{aligned}
S\_5 &\gets S\_1 + S\_2 \\\\
S\_6 &\gets S\_1 - S\_2 \\\\
S\_7 &\gets 2S\_3 \\\\
S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\
S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4
\end{aligned}
$$
Performing too many intermediate additions and subtractions grows
the bounds beyond what is allowed as input to multiplication,
forcing an extra carry pass. However, it is just possible to avoid
this by rearranging signs.
Assume that the bounds on the limbs of each field element are
parameterized by \\( b \in \mathbb R \\) representing the excess
bits, so that each limb is bounded by either
\\( 2\^{25+b} \\) or \\( 2\^{26+b} \\).
The multiplication routine requires that its inputs are bounded by
\\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\) \\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\)
into 32 bits. Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\) into 32 bits. Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\)
when \\( x < 2\^{27.75} = 2\^{26 + 1.75} \\). However, this is only 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 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} \begin{matrix}
& S\_1 & S\_1 & S\_1 & S\_1 \\\\ & S\_1 & S\_1 & S\_1 & S\_1 \\\\
@ -374,24 +312,22 @@ $$
=& S\_5 & S\_6 & S\_8 & S\_9 =& S\_5 & S\_6 & S\_8 & S\_9
\end{matrix} \end{matrix}
$$ $$
results in bit-excesses \\( < (1.01, 1.60, 2.33, 2.01)\\) for
results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for
\\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute
are then are then
$$ $$
\begin{aligned} \begin{aligned}
X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.00) \\\\ X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.01) \\\\
Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\
Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\
T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 2.00) T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 2.01)
\end{aligned} \end{aligned}
$$ $$
which are too large: it's not possible to arrange the multiplicands so
which are too large. However, if we flip the sign of \\( S\_4 = that one vector has \\(b < 2.5\\) and the other has \\( b < 1.75 \\).
S\_0\^2 \\) during squaring, so that we output \\(S\_4' = -S\_4 However, if we flip the sign of \\( S\_4 = S\_0\^2 \\) during
\pmod p\\), then we can compute squaring, so that we output \\(S\_4' = -S\_4 \pmod p\\), then we can
compute
$$ $$
\begin{matrix} \begin{matrix}
& S\_1 & S\_1 & S\_1 & S\_1 \\\\ & S\_1 & S\_1 & S\_1 & S\_1 \\\\
@ -404,61 +340,120 @@ $$
=& S\_5 & S\_6 & S\_8 & S\_9 =& S\_5 & S\_6 & S\_8 & S\_9
\end{matrix} \end{matrix}
$$ $$
resulting in bit-excesses \\( < (1.01, 1.60, 2.33, 1.60)\\) for
resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for
\\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute
are then are then
$$ $$
\begin{aligned} \begin{aligned}
X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.59) \\\\ X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.60) \\\\
Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\
Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\
T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 1.59) T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 1.60)
\end{aligned} \end{aligned}
$$ $$
whose right-hand sides are all bounded with \\( b < 1.75 \\) and 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 # Comparison to non-vectorized formulas
HWCD also suggest using a mixed representation, passing between \\( In theory, the parallel Edwards formulas seem to allow a \\(4\\)-way
\mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\) speedup from parallelism. However, an actual vectorized
"projective" coordinates, where doubling is slightly cheaper (saving implementation has several slowdowns that cut into this speedup.
about \\(\mathbf 1M\\). This approach is used for the
non-vectorized `u32` and `u64` backends, and more
details on the different coordinate systems can be found in the
`curve_models` module documentation.
This optimization is not compatible with the parallel formulas, which are First, the parallel formulas can only use a \\( 32 \times 32
therefore slightly less efficient when counting the total number of
field multiplications and squarings. In particular, vectorized doublings
are less efficient than serial doublings.
In addition, the parallel formulas can only use a \\( 32 \times 32
\rightarrow 64 \\)-bit integer multiplier, so the speedup from \rightarrow 64 \\)-bit integer multiplier, so the speedup from
vectorization must overcome the disadvantage of losing the \\( 64 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 Second, the parallel doubling formulas incur both a theoretical and
`ymm16..ymm31` registers to reduce register pressure, and avoid practical slowdown. The parallel formulas described above work on the
spills during field multiplication. This gives a small but \\( \mathbb P\^3 \\) “extended” coordinates. The \\( \mathbb P\^2 \\)
noticeable speedup. 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 Unfortunately, this optimization is not compatible with the parallel
(particularly Skylake and Skylake-X microarchitectures) perform thermal formulas, which cannot save time by skipping a single variable, so the
throttling when using wide vector instructions. For a mixed workload, 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 where point operations are interspersed with other tasks, this can
reduce overall performance. This probably means that this reduce overall performance. This implementation is therefore probably
implementation is not suitable for basic applications, like signatures, not suitable for basic applications, like signatures, but is
but could still be worthwhile for complex applications, like worthwhile for complex applications, like zero-knowledge proofs, which
zero-knowledge proofs, which do enough work to make it worthwhile. 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 [sandy2x]: https://eprint.iacr.org/2015/943.pdf
[avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28 [avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28
[hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf [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

26
src/backend/avx2/constants.rs
View file

@ -16,6 +16,24 @@ use backend::avx2::edwards::{CachedPoint, ExtendedPoint};
use backend::avx2::field::FieldElement32x4; use backend::avx2::field::FieldElement32x4;
use scalar_mul::window::NafLookupTable8; 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 /// The low limbs of (2p, 2p, 2p, 2p), so that
/// ```no_run /// ```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] /// (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, 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: /// Odd multiples of the Ed25519 basepoint:
pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: NafLookupTable8<CachedPoint> = NafLookupTable8([ pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: NafLookupTable8<CachedPoint> = NafLookupTable8([
CachedPoint(FieldElement32x4([ CachedPoint(FieldElement32x4([

266
src/backend/avx2/edwards.rs
View file

@ -8,29 +8,53 @@
// - Isis Agora Lovecruft <isis@patternsinthevoid.net> // - Isis Agora Lovecruft <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca> // - Henry de Valence <hdevalence@hdevalence.ca>
//! 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(non_snake_case)]
#![allow(bad_style)]
use core::convert::From; 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::Choice;
use subtle::ConditionallyAssignable;
use edwards; use edwards;
use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8};
use traits::Identity; 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, using parallel Edwards formulas for curve
/// operations.
/// A point on Curve25519, represented in an AVX2-friendly format. ///
/// # Invariant
///
/// The coefficients of an `ExtendedPoint` are bounded with
/// \\( b < 0.007 \\).
#[derive(Copy, Clone, Debug)] #[derive(Copy, Clone, Debug)]
pub struct ExtendedPoint(pub(super) FieldElement32x4); pub struct ExtendedPoint(pub(super) FieldElement32x4);
@ -43,7 +67,12 @@ impl From<edwards::EdwardsPoint> for ExtendedPoint {
impl From<ExtendedPoint> for edwards::EdwardsPoint { impl From<ExtendedPoint> for edwards::EdwardsPoint {
fn from(P: ExtendedPoint) -> edwards::EdwardsPoint { fn from(P: ExtendedPoint) -> edwards::EdwardsPoint {
let tmp = P.0.split(); 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,74 +90,33 @@ impl Default for ExtendedPoint {
impl Identity for ExtendedPoint { impl Identity for ExtendedPoint {
fn identity() -> ExtendedPoint { fn identity() -> ExtendedPoint {
ExtendedPoint(FieldElement32x4([ constants::EXTENDEDPOINT_IDENTITY
u32x8::new(0,1,0,0,1,0,0,0),
u32x8::splat(0),
u32x8::splat(0),
u32x8::splat(0),
u32x8::splat(0),
]))
} }
} }
impl ExtendedPoint { impl ExtendedPoint {
/// Compute the double of this point.
pub fn double(&self) -> ExtendedPoint { 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;
let P = &self.0;
let mut t0 = FieldElement32x4::zero();
let mut t1 = FieldElement32x4::zero();
// Want to compute (X1 Y1 Z1 X1+Y1). // Want to compute (X1 Y1 Z1 X1+Y1).
// Not sure how to do this less expensively than computing // Not sure how to do this less expensively than computing
// (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1) // (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1)
// (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1) // (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1)
// and then adding. // and then adding.
// Set t0 = (X1 Y1 X1 Y1) // Set tmp0 = (X1 Y1 X1 Y1)
t0.0[0] = _mm256_permute2x128_si256(P.0[0].into_bits(), P.0[0].into_bits(), 0b0000_0000).into_bits(); let mut tmp0 = self.0.shuffle(Shuffle::ABAB);
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 t1 = (Y1 X1 Y1 X1) // Set tmp1 = (Y1 X1 Y1 X1)
t1.0[0] = _mm256_shuffle_epi32(t0.0[0].into_bits(), 0b10_11_00_01).into_bits(); let mut tmp1 = tmp0.shuffle(Shuffle::BADC);
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();
// Set t0 = (X1+Y1 X1+Y1 X1+Y1 X1+Y1) // Set tmp0 = (X1 Y1 Z1 X1+Y1)
t0.0[0] = t0.0[0] + t1.0[0]; tmp0 = self.0.blend(tmp0 + tmp1, Lanes::D);
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];
// Set t0 = (X1 Y1 Z1 X1+Y1) // Set tmp1 = tmp0^2, negating the D values
// why does this intrinsic take an i32 for the imm8 ??? tmp1 = tmp0.square_and_negate_D();
t0.0[0] = _mm256_blend_epi32(P.0[0].into_bits(), t0.0[0].into_bits(), D_LANES as i32).into_bits(); // Now tmp1 = (S1 S2 S3 -S4) with b < 0.007
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();
// Set t1 = t0^2, negating the D values
t1 = t0.square_and_negate_D();
// 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. // See discussion of bounds in the module-level documentation.
//
// We want to compute // We want to compute
// //
// + | S1 | S1 | S1 | S1 | // + | S1 | S1 | S1 | S1 |
@ -140,32 +128,33 @@ impl ExtendedPoint {
// - | | S2 | S2 | | // - | | S2 | S2 | |
// ======================= // =======================
// S5 S6 S8 S9 // 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 zero = FieldElement32x4::zero();
let c1 = u32x8::new(5,1,7,3,1,5,3,7).into_bits(); // (ABCD) -> (DBBD) let S_1 = tmp1.shuffle(Shuffle::AAAA);
let S_2 = tmp1.shuffle(Shuffle::BBBB);
for i in 0..5 { tmp0 = zero.blend(tmp1 + tmp1, Lanes::C);
let tmp = t0.0[i]; // tmp0 = (0, 0, 2S_3, 0)
t0.0[i] = _mm256_permutevar8x32_epi32(tmp.into_bits(), c0).into_bits(); tmp0 = tmp0.blend(tmp1, Lanes::D);
t1.0[i] = _mm256_permutevar8x32_epi32(tmp.into_bits(), c1).into_bits(); // 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)
ExtendedPoint(&t0 * &t1) // 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);
// 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 { 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. /// 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)] #[derive(Copy, Clone, Debug)]
pub struct CachedPoint(pub(super) FieldElement32x4); pub struct CachedPoint(pub(super) FieldElement32x4);
@ -185,15 +183,16 @@ impl From<ExtendedPoint> for CachedPoint {
fn from(P: ExtendedPoint) -> CachedPoint { fn from(P: ExtendedPoint) -> CachedPoint {
let mut x = P.0; let mut x = P.0;
// x = (S2 S3 Z2 T2) x = x.blend(x.diff_sum(), Lanes::AB);
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 = (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 = (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) CachedPoint(x)
} }
} }
@ -206,13 +205,7 @@ impl Default for CachedPoint {
impl Identity for CachedPoint { impl Identity for CachedPoint {
fn identity() -> CachedPoint { fn identity() -> CachedPoint {
CachedPoint(FieldElement32x4([ constants::CACHEDPOINT_IDENTITY
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),
]))
} }
} }
@ -224,62 +217,63 @@ impl ConditionallyAssignable for CachedPoint {
impl<'a> Neg for &'a CachedPoint { impl<'a> Neg for &'a CachedPoint {
type Output = 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 { fn neg(self) -> CachedPoint {
let mut neg = *self; let swapped = self.0.shuffle(Shuffle::BACD);
neg.0.swap_AB(); CachedPoint(swapped.blend(swapped.negate_lazy(), Lanes::D))
neg.0.negate_D_lazy();
neg
} }
} }
impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint {
type Output = 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 { fn add(self, other: &'b CachedPoint) -> ExtendedPoint {
unsafe { // The coefficients of an `ExtendedPoint` are reduced after
use core::arch::x86_64::_mm256_permutevar8x32_epi32; // 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 = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) tmp = tmp.blend(tmp.diff_sum(), Lanes::AB);
tmp.diff_sum(Lanes::AB); // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) with b < 1.6
// tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) // (tmp, other) bounded with b < (1.6, 1.0) < (2.5, 1.75).
tmp = &tmp * &other.0; tmp = &tmp * &other.0;
// tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11)
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 c0 = u32x8::new(0,5,2,7,5,0,7,2); // (ABCD) -> (ADDA) let t0 = tmp.shuffle(Shuffle::ADDA);
let c1 = u32x8::new(4,1,6,3,4,1,6,3); // (ABCD) -> (CBCB) // t0 = (S12 S15 S15 S12)
let t1 = tmp.shuffle(Shuffle::CBCB);
// t1 = (S14 S13 S14 S13)
// set t0 = (S12 S15 S15 S12) // All coefficients of t0, t1 are bounded with b < 1.6.
// set t1 = (S14 S13 S14 S13) // Return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3)
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) ExtendedPoint(&t0 * &t1)
} }
} }
}
impl<'a, 'b> Sub<&'b CachedPoint> for &'a ExtendedPoint { impl<'a, 'b> Sub<&'b CachedPoint> for &'a ExtendedPoint {
type Output = ExtendedPoint; type Output = ExtendedPoint;
/// Implement subtraction by negating the point and adding. /// Implement subtraction by negating the point and adding.
/// ///
/// Empirically, this seems about the same cost as a custom subtraction impl (maybe because the /// Empirically, this seems about the same cost as a custom
/// benefit is cancelled by increased code size?) /// subtraction impl (maybe because the benefit is cancelled by
/// increased code size?)
fn sub(self, other: &'b CachedPoint) -> ExtendedPoint { fn sub(self, other: &'b CachedPoint) -> ExtendedPoint {
self + &(-other) self + &(-other)
} }
@ -335,7 +329,7 @@ mod test {
macro_rules! print_var { macro_rules! print_var {
($x:ident) => { ($x:ident) => {
println!("{} = {:?}", stringify!($x), $x.to_bytes()); println!("{} = {:?}", stringify!($x), $x.to_bytes());
} };
} }
let S0 = &Y1 - &X1; // R1 let S0 = &Y1 - &X1; // R1
@ -383,7 +377,12 @@ mod test {
let Z3 = &S15 * &S14; // R2 * R3 let Z3 = &S15 * &S14; // R2 * R3
let T3 = &S12 * &S13; // R1 * R4 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) { fn addition_test_helper(P: edwards::EdwardsPoint, Q: edwards::EdwardsPoint) {
@ -437,12 +436,12 @@ mod test {
} }
fn serial_double(P: edwards::EdwardsPoint) -> edwards::EdwardsPoint { 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 { macro_rules! print_var {
($x:ident) => { ($x:ident) => {
println!("{} = {:?}", stringify!($x), $x.to_bytes()); println!("{} = {:?}", stringify!($x), $x.to_bytes());
} };
} }
let S0 = &X1 + &Y1; // R1 let S0 = &X1 + &Y1; // R1
@ -476,7 +475,12 @@ mod test {
let Z3 = &S8 * &S6; let Z3 = &S8 * &S6;
let T3 = &S5 * &S9; 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) { fn doubling_test_helper(P: edwards::EdwardsPoint) {

780
src/backend/avx2/field.rs
View file

@ -8,54 +8,143 @@
// - Isis Agora Lovecruft <isis@patternsinthevoid.net> // - Isis Agora Lovecruft <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca> // - Henry de Valence <hdevalence@hdevalence.ca>
//! 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; const A_LANES: u8 = 0b0000_0101;
pub const B_LANES: u8 = 0b0000_1010; const B_LANES: u8 = 0b0000_1010;
pub const C_LANES: u8 = 0b0101_0000; const C_LANES: u8 = 0b0101_0000;
pub const D_LANES: u8 = 0b1010_0000; const D_LANES: u8 = 0b1010_0000;
pub const A_LANES64: u8 = 0b00_00_00_11; #[allow(unused)]
pub const B_LANES64: u8 = 0b00_00_11_00; const A_LANES64: u8 = 0b00_00_00_11;
pub const C_LANES64: u8 = 0b00_11_00_00; #[allow(unused)]
pub const D_LANES64: u8 = 0b11_00_00_00; 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::{Add, Mul, Neg};
use core::simd::{i32x8, u32x8, u64x4, IntoBits};
use core::ops::Mul;
use core::simd::{IntoBits, u32x8, i32x8, u64x4};
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::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)] /// Unpack 32-bit lanes into 64-bit lanes:
pub enum Lanes { /// ```
AB, /// (a0, b0, a1, b1, c0, d0, c1, d1)
CD, /// ```
ALL, /// 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)] #[inline(always)]
fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { fn repack_pair(x: u32x8, y: u32x8) -> u32x8 {
unsafe { unsafe {
use core::arch::x86_64::_mm256_blend_epi32; use core::arch::x86_64::_mm256_blend_epi32;
use core::arch::x86_64::_mm256_shuffle_epi32;
match control { // Input: x = (a0, 0, b0, 0, c0, 0, d0, 0)
Lanes::AB => _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32).into_bits(), // Input: y = (a1, 0, b1, 0, c1, 0, d1, 0)
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(), 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. /// 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)] #[derive(Clone, Copy, Debug)]
pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]); pub struct FieldElement32x4(pub(crate) [u32x8; 5]);
use subtle::ConditionallyAssignable;
use subtle::Choice; use subtle::Choice;
use subtle::ConditionallyAssignable;
impl ConditionallyAssignable for FieldElement32x4 { impl ConditionallyAssignable for FieldElement32x4 {
fn conditional_assign(&mut self, other: &FieldElement32x4, choice: Choice) { fn conditional_assign(&mut self, other: &FieldElement32x4, choice: Choice) {
@ -68,10 +157,11 @@ impl ConditionallyAssignable for FieldElement32x4 {
} }
impl 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]; let mut out = [FieldElement64::zero(); 4];
for i in 0..5 { for i in 0..5 {
let a_2i = self.0[i].extract(0) as u64; // let a_2i = self.0[i].extract(0) as u64; //
let b_2i = self.0[i].extract(1) as u64; // let b_2i = self.0[i].extract(1) as u64; //
let a_2i_1 = self.0[i].extract(2) as u64; // `. let a_2i_1 = self.0[i].extract(2) as u64; // `.
@ -90,14 +180,138 @@ impl FieldElement32x4 {
out 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 { 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 { 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( pub fn new(
x0: &FieldElement64, x0: &FieldElement64,
x1: &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); 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); // We don't know that the original `FieldElement64`s were
out.reduce32(); // fully reduced, so the odd limbs may exceed 2^25.
return out; // 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. /// # Preconditions
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)\\).
/// ///
/// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. /// The coefficients of `self` must be bounded with \\( b < 0.999 \\).
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`.
/// ///
/// This is `#[inline(always)]` because the `mask` parameter should be an immediate. /// # Postconditions
#[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) \\).
/// ///
/// Compute /// The coefficients of the result are bounded with \\( b < 1 \\).
/// $$( 121666A, 121666B, 2\cdot 121666C, 2\cdot 121665 D).$$ #[inline]
pub fn scale_by_curve_constants(&mut self) { pub fn negate_lazy(&self) -> FieldElement32x4 {
let mut b = [u64x4::splat(0); 10]; // The limbs of self are bounded with b < 0.999, while the
// smallest limb of 2*p is 67108845 > 2^{26+0.9999}, so
let consts = u32x8::new(121666, 0, 121666, 0, 2*121666, 0, 2*121665, 0); // underflows are not possible.
FieldElement32x4([
unsafe { P_TIMES_2_LO - self.0[0],
use core::arch::x86_64::_mm256_mul_epu32; P_TIMES_2_HI - self.0[1],
P_TIMES_2_HI - self.0[2],
let (b0, b1) = unpack_pair(self.0[0]); P_TIMES_2_HI - self.0[3],
b[0] = _mm256_mul_epu32(b0.into_bits(), consts.into_bits()).into_bits(); P_TIMES_2_HI - self.0[4],
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); /// 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
} }
pub fn reduce32(&mut self) { /// 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 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, 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,
(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 { unsafe {
use core::arch::x86_64::_mm256_srlv_epi32; 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; 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 { let combine = |v_lo: u32x8, v_hi: u32x8| -> u32x8 {
unsafe { unsafe {
use core::arch::x86_64::_mm256_blend_epi32; use core::arch::x86_64::_mm256_blend_epi32;
@ -289,35 +438,65 @@ 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); 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); 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); 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); 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); v[4] = (v[4] & masks) + combine(c76, c98);
// Still need to account for c9 let c9_19: u32x8 = unsafe {
// c98 = (c9, c9, c8, c8, c9, c9, c8, c8)
//
let c9_19: u32x8;
unsafe {
use core::arch::x86_64::_mm256_mul_epu32; use core::arch::x86_64::_mm256_mul_epu32;
use core::arch::x86_64::_mm256_shuffle_epi32; use core::arch::x86_64::_mm256_shuffle_epi32;
// 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); 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()); 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();
}
// 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; 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 // These aren't const because splat isn't a const fn
let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1); let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1);
let LOW_26_BITS: u64x4 = u64x4::splat((1 << 26) - 1); let LOW_26_BITS: u64x4 = u64x4::splat((1 << 26) - 1);
@ -370,8 +549,13 @@ impl FieldElement32x4 {
z[1] = z[1] + c1; // z1 < 2^25 + 2^17.25 < 2^25.0067 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 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([ FieldElement32x4([
repack_pair(z[0].into_bits(), z[1].into_bits()), repack_pair(z[0].into_bits(), z[1].into_bits()),
repack_pair(z[2].into_bits(), z[3].into_bits()), repack_pair(z[2].into_bits(), z[3].into_bits()),
@ -380,52 +564,17 @@ impl FieldElement32x4 {
repack_pair(z[8].into_bits(), z[9].into_bits()), repack_pair(z[8].into_bits(), z[9].into_bits()),
]) ])
} }
}
#[inline(always)] /// Square this field element, and negate the result's \\(D\\) value.
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.
/// ///
/// # 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 { pub fn square_and_negate_D(&self) -> FieldElement32x4 {
let neg_mask = D_LANES64;
#[inline(always)] #[inline(always)]
fn m(x: u32x8, y: u32x8) -> u64x4 { fn m(x: u32x8, y: u32x8) -> u64x4 {
use core::arch::x86_64::_mm256_mul_epu32; use core::arch::x86_64::_mm256_mul_epu32;
@ -514,9 +663,99 @@ 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<FieldElement32x4> 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 { impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 {
type Output = 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)] #[inline(always)]
fn m(x: u32x8, y: u32x8) -> u64x4 { fn m(x: u32x8, y: u32x8) -> u64x4 {
use core::arch::x86_64::_mm256_mul_epu32; use core::arch::x86_64::_mm256_mul_epu32;
@ -535,11 +774,11 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 {
let (x6, x7) = unpack_pair(self.0[3]); let (x6, x7) = unpack_pair(self.0[3]);
let (x8, x9) = unpack_pair(self.0[4]); let (x8, x9) = unpack_pair(self.0[4]);
let (y0, y1) = unpack_pair(_rhs.0[0]); let (y0, y1) = unpack_pair(rhs.0[0]);
let (y2, y3) = unpack_pair(_rhs.0[1]); let (y2, y3) = unpack_pair(rhs.0[1]);
let (y4, y5) = unpack_pair(_rhs.0[2]); let (y4, y5) = unpack_pair(rhs.0[2]);
let (y6, y7) = unpack_pair(_rhs.0[3]); let (y6, y7) = unpack_pair(rhs.0[3]);
let (y8, y9) = unpack_pair(_rhs.0[4]); 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);
@ -547,9 +786,9 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 {
let y2_19 = m_lo(v19, y2); // iff 26 + b + lg(19) < 32 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 y3_19 = m_lo(v19, y3); // if b < 32 - 26 - 4.248 = 1.752
let y4_19 = m_lo(v19, y4); let y4_19 = m_lo(v19, y4);
let y5_19 = m_lo(v19, y5); // below, b<2.5: this is a bottleneck, let y5_19 = m_lo(v19, y5);
let y6_19 = m_lo(v19, y6); // could be avoided by promoting to let y6_19 = m_lo(v19, y6);
let y7_19 = m_lo(v19, y7); // u64 here instead of in m() let y7_19 = m_lo(v19, y7);
let y8_19 = m_lo(v19, y8); let y8_19 = m_lo(v19, y8);
let y9_19 = m_lo(v19, y9); 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 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); 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]) FieldElement32x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9])
} }
} }
@ -582,7 +859,8 @@ mod test {
#[test] #[test]
fn scale_by_curve_constants() { fn scale_by_curve_constants() {
let mut x = FieldElement32x4::splat(&FieldElement64::one()); let mut x = FieldElement32x4::splat(&FieldElement64::one());
x.scale_by_curve_constants();
x = x * (121666, 121666, 2*121666, 2*121665);
let xs = x.split(); let xs = x.split();
assert_eq!(xs[0], FieldElement64([121666, 0, 0, 0, 0])); assert_eq!(xs[0], FieldElement64([121666, 0, 0, 0, 0]));
@ -598,8 +876,7 @@ mod test {
let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]); let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]);
let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]);
let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3).diff_sum();
vec.diff_sum(Lanes::ALL);
let result = vec.split(); let result = vec.split();
@ -607,16 +884,6 @@ mod test {
assert_eq!(result[1], &x1 + &x0); assert_eq!(result[1], &x1 + &x0);
assert_eq!(result[2], &x3 - &x2); assert_eq!(result[2], &x3 - &x2);
assert_eq!(result[3], &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] #[test]
@ -636,7 +903,6 @@ mod test {
assert_eq!(result[3], -&(&x3 * &x3)); assert_eq!(result[3], -&(&x3 * &x3));
} }
#[test] #[test]
fn multiply_vs_serial() { fn multiply_vs_serial() {
let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]); let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]);

2
src/backend/avx2/scalar_mul/straus.rs
View file

@ -8,6 +8,8 @@
// - Isis Agora Lovecruft <isis@patternsinthevoid.net> // - Isis Agora Lovecruft <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca> // - Henry de Valence <hdevalence@hdevalence.ca>
#![allow(non_snake_case)]
use core::borrow::Borrow; use core::borrow::Borrow;
use clear_on_drop::ClearOnDrop; use clear_on_drop::ClearOnDrop;

2
src/backend/mod.rs
View file

@ -27,6 +27,6 @@ pub mod u32;
#[cfg(feature = "u64_backend")] #[cfg(feature = "u64_backend")]
pub mod u64; pub mod u64;
#[cfg(all(feature = "avx2_backend", feature = "yolocrypto", target_feature = "avx2"))] #[cfg(all(feature = "avx2_backend", target_feature = "avx2"))]
pub mod avx2; pub mod avx2;

2
src/lib.rs
View file

@ -14,7 +14,7 @@
#![cfg_attr(feature = "nightly", feature(cfg_target_feature))] #![cfg_attr(feature = "nightly", feature(cfg_target_feature))]
#![cfg_attr(feature = "nightly", feature(external_doc))] #![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. // Refuse to compile if documentation is missing, but only on nightly.
// //