From 7e0ddf6b98ea1b357d1653eae31e797789501849 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 4 Apr 2018 21:13:09 -0700 Subject: [PATCH 1/9] Rewrite NAF code to work with more window sizes Change Scalar::non_adjacent_form() to take a width parameter. This rewrite also makes it faster, although it's probably a ways off from optimal. I don't know how much it matters. TODO: write up description of why this computes the same thing. Thanks to @oleganza for pointing out an error reading bits across words in an earlier version of this code. --- Cargo.toml | 2 + build.rs | 1 + .../avx2/scalar_mul/vartime_double_base.rs | 4 +- src/backend/avx2/scalar_mul/vartime_straus.rs | 2 +- src/lib.rs | 2 + src/scalar.rs | 78 +++++++++++-------- src/scalar_mul/vartime_double_base.rs | 4 +- src/scalar_mul/vartime_straus.rs | 2 +- 8 files changed, 58 insertions(+), 37 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 9cb39f1..749700b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -41,6 +41,7 @@ harness = false # match exactly, since the build.rs uses the crate itself as a library. [dependencies] +byteorder = "1" digest = "0.7" generic-array = "0.9" clear_on_drop = "=0.2.3" @@ -49,6 +50,7 @@ serde = { version = "1.0", optional = true } rand = { version = "0.4", optional = true } [build-dependencies] +byteorder = "1" digest = "0.7" generic-array = "0.9" clear_on_drop = "=0.2.3" diff --git a/build.rs b/build.rs index 8331688..595d369 100644 --- a/build.rs +++ b/build.rs @@ -5,6 +5,7 @@ #![allow(non_snake_case)] #![allow(dead_code)] +extern crate byteorder; extern crate clear_on_drop; extern crate core; extern crate digest; diff --git a/src/backend/avx2/scalar_mul/vartime_double_base.rs b/src/backend/avx2/scalar_mul/vartime_double_base.rs index 1ba09f1..aa2aa97 100644 --- a/src/backend/avx2/scalar_mul/vartime_double_base.rs +++ b/src/backend/avx2/scalar_mul/vartime_double_base.rs @@ -18,8 +18,8 @@ use backend::avx2::constants::BASEPOINT_ODD_LOOKUP_TABLE; /// Compute \\(aA + bB\\) in variable time, where \\(B\\) is the Ed25519 basepoint. pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { - let a_naf = a.non_adjacent_form(); - let b_naf = b.non_adjacent_form(); + let a_naf = a.non_adjacent_form(5); + let b_naf = b.non_adjacent_form(5); // Find starting index let mut i: usize = 255; diff --git a/src/backend/avx2/scalar_mul/vartime_straus.rs b/src/backend/avx2/scalar_mul/vartime_straus.rs index 224713d..f05754b 100644 --- a/src/backend/avx2/scalar_mul/vartime_straus.rs +++ b/src/backend/avx2/scalar_mul/vartime_straus.rs @@ -27,7 +27,7 @@ where { let nafs: Vec<_> = scalars .into_iter() - .map(|c| c.borrow().non_adjacent_form()) + .map(|c| c.borrow().non_adjacent_form(5)) .collect(); let lookup_tables: Vec<_> = points .into_iter() diff --git a/src/lib.rs b/src/lib.rs index 0a8a1a3..3a3cb84 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -41,6 +41,8 @@ extern crate alloc; extern crate clear_on_drop; +extern crate byteorder; + // The `Digest` trait is implemented using `generic_array`, so we need it // too. Hopefully we can eliminate `generic_array` from `Digest` once const // generics land. diff --git a/src/scalar.rs b/src/scalar.rs index 7dae1e1..f4a86dd 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -493,7 +493,7 @@ impl Scalar { bits } - /// Compute a width-5 "Non-Adjacent Form" of this scalar. + /// Compute a width-\\(w\\) "Non-Adjacent Form" of this scalar. /// /// A width-\\(w\\) NAF of a positive integer \\(k\\) is an expression /// $$ @@ -507,37 +507,53 @@ impl Scalar { /// Intuitively, this is like a binary expansion, except that we /// allow some coefficients to grow up to \\(2\^{w-1}\\) so that the /// nonzero coefficients are as sparse as possible. - pub(crate) fn non_adjacent_form(&self) -> [i8; 256] { - // Step 1: write out bits of the scalar - let mut naf = self.bits(); + pub(crate) fn non_adjacent_form(&self, w: usize) -> [i8; 256] { + use byteorder::{ByteOrder, LittleEndian}; - // Step 2: zero coefficients by carrying them upwards or downwards - 'bits: for i in 0..256 { - if naf[i] == 0 { continue 'bits; } - 'window: for b in 1..6 { - if i+b >= 256 { break 'window; } - if naf[i+b] == 0 { continue 'window; } - let potential_carry = naf[i+b] << b; - if naf[i+b] + potential_carry <= 15 { - // Eliminate naf[i+b] by carrying its value onto naf[i] - naf[i] += potential_carry; - naf[i+b] = 0; - } else if naf[i+b] - potential_carry >= -15 { - // Eliminate naf[i+b] by carrying its value upwards. - naf[i] -= potential_carry; // Subtract 2^(i+b) - 'carry: for k in i+b..256 { - if naf[k] != 0 { - // Since naf[k] = 0 or 1 for k > i, naf[k] == 1. - naf[k] = 0; // Subtract 2^k - } else { - // By now we have subtracted 2^k = - // 2^(i+b) + 2^(i+b) + 2^(i+b+1) + ... + 2^(k-1). - naf[k] = 1; // Add back 2^k. - break 'carry; - } - } - } + let mut naf = [0i8; 256]; + + let mut x_u64 = [0u64; 5]; + LittleEndian::read_u64_into(&self.bytes, &mut x_u64[0..4]); + + let width = 1 << w; + let window_mask = width - 1; + + let mut pos = 0; + let mut carry = 0; + while pos < 256 { + // Construct a buffer of bits of the scalar, starting at bit `pos` + let u64_idx = pos / 64; + let bit_idx = pos % 64; + let bit_buf: u64; + if bit_idx < 64 - w { + // This window's bits are contained in a single u64 + bit_buf = x_u64[u64_idx] >> bit_idx; + } else { + // Combine the current u64's bits with the bits from the next u64 + bit_buf = (x_u64[u64_idx] >> bit_idx) | (x_u64[1+u64_idx] << (64 - bit_idx)); } + + // Add the carry into the current window + let window = carry + (bit_buf & window_mask); + + if window & 1 == 0 { + // If the window value is even, preserve the carry and continue. + // Why is the carry preserved? + // If carry == 0 and window & 1 == 0, then the next carry should be 0 + // If carry == 1 and window & 1 == 0, then bit_buf & 1 == 1 so the next carry should be 1 + pos += 1; + continue; + } + + if window < width/2 { + carry = 0; + naf[pos] = window as i8; + } else { + carry = 1; + naf[pos] = (window as i8) - (width as i8); + } + + pos += w; } naf @@ -789,7 +805,7 @@ mod test { #[test] fn non_adjacent_form() { - let naf = A_SCALAR.non_adjacent_form(); + let naf = A_SCALAR.non_adjacent_form(5); for i in 0..256 { assert_eq!(naf[i], A_NAF[i]); } diff --git a/src/scalar_mul/vartime_double_base.rs b/src/scalar_mul/vartime_double_base.rs index 60c51e6..9fcd267 100644 --- a/src/scalar_mul/vartime_double_base.rs +++ b/src/scalar_mul/vartime_double_base.rs @@ -18,8 +18,8 @@ use scalar_mul::window::OddLookupTable; /// Compute \\(aA + bB\\) in variable time, where \\(B\\) is the Ed25519 basepoint. pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { - let a_naf = a.non_adjacent_form(); - let b_naf = b.non_adjacent_form(); + let a_naf = a.non_adjacent_form(5); + let b_naf = b.non_adjacent_form(5); // Find starting index let mut i: usize = 255; diff --git a/src/scalar_mul/vartime_straus.rs b/src/scalar_mul/vartime_straus.rs index c55ce21..b419c0a 100644 --- a/src/scalar_mul/vartime_straus.rs +++ b/src/scalar_mul/vartime_straus.rs @@ -27,7 +27,7 @@ where { let nafs: Vec<_> = scalars .into_iter() - .map(|c| c.borrow().non_adjacent_form()) + .map(|c| c.borrow().non_adjacent_form(5)) .collect(); let lookup_tables: Vec<_> = points .into_iter() From 753a0292ded69b34bc69687e6c8e25c5a907bd9c Mon Sep 17 00:00:00 2001 From: Henry & Isis Date: Thu, 5 Apr 2018 04:48:28 +0000 Subject: [PATCH 2/9] Rename OddLookupTable to NafLookupTable5. An OddLookupTable corresponds to a non-adjacent form of width 5. --- build.rs | 8 ++++---- src/backend/avx2/constants.rs | 6 +++--- src/backend/avx2/edwards.rs | 6 +++--- .../avx2/scalar_mul/vartime_double_base.rs | 4 ++-- src/backend/avx2/scalar_mul/vartime_straus.rs | 4 ++-- src/scalar_mul/vartime_double_base.rs | 4 ++-- src/scalar_mul/vartime_straus.rs | 4 ++-- src/scalar_mul/window.rs | 16 ++++++++-------- 8 files changed, 26 insertions(+), 26 deletions(-) diff --git a/build.rs b/build.rs index 595d369..7d06199 100644 --- a/build.rs +++ b/build.rs @@ -59,7 +59,7 @@ mod scalar_mul; use edwards::EdwardsBasepointTable; use curve_models::AffineNielsPoint; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; fn main() { // Enable the "precomputed_tables" feature in the main build stage @@ -86,7 +86,7 @@ use edwards::EdwardsBasepointTable; use curve_models::AffineNielsPoint; use scalar_mul::window::LookupTable; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; /// Table containing precomputed multiples of the Ed25519 basepoint \\\\(B = (x, 4/5)\\\\). pub const ED25519_BASEPOINT_TABLE: EdwardsBasepointTable = ED25519_BASEPOINT_TABLE_INNER_DOC_HIDDEN; @@ -101,13 +101,13 @@ pub const ED25519_BASEPOINT_TABLE_INNER_DOC_HIDDEN: EdwardsBasepointTable = {:?} // Now generate AFFINE_ODD_MULTIPLES_OF_BASEPOINT let B = &constants::ED25519_BASEPOINT_POINT; - let odd_multiples = OddLookupTable::::from(B); + let odd_multiples = NafLookupTable5::::from(B); f.write_all( format!( "\n /// Odd multiples of the basepoint `[B, 3B, 5B, 7B, 9B, 11B, 13B, 15B]`. -pub(crate) const AFFINE_ODD_MULTIPLES_OF_BASEPOINT: OddLookupTable = {:?}; +pub(crate) const AFFINE_ODD_MULTIPLES_OF_BASEPOINT: NafLookupTable5 = {:?}; \n\n", &odd_multiples ).as_bytes(), diff --git a/src/backend/avx2/constants.rs b/src/backend/avx2/constants.rs index e0479db..0ef2300 100644 --- a/src/backend/avx2/constants.rs +++ b/src/backend/avx2/constants.rs @@ -12,7 +12,7 @@ use core::simd::u32x8; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; use backend::avx2::field::FieldElement32x4; use backend::avx2::edwards::{ExtendedPoint, CachedPoint}; @@ -53,8 +53,8 @@ pub(crate) static P_TIMES_2_MASKED: FieldElement32x4 = FieldElement32x4([ ]); /// Odd multiples of the Ed25519 basepoint: -pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: OddLookupTable = - OddLookupTable([ +pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: NafLookupTable5 = + NafLookupTable5([ CachedPoint(FieldElement32x4([ u32x8::new(3571425, 10045002, 19036563, 1096096, 243332, 65897020, 0, 28963681), u32x8::new(30896895, 63055514, 1614915, 5095970, 0, 53791688, 0, 31258312), diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index c0ad9ac..3c52bdb 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -22,7 +22,7 @@ use subtle::ConditionallyAssignable; use subtle::Choice; use edwards; -use scalar_mul::window::{LookupTable, OddLookupTable}; +use scalar_mul::window::{LookupTable, NafLookupTable5}; use traits::Identity; @@ -294,7 +294,7 @@ impl<'a> From<&'a edwards::EdwardsPoint> for LookupTable { } } -impl<'a> From<&'a edwards::EdwardsPoint> for OddLookupTable { +impl<'a> From<&'a edwards::EdwardsPoint> for NafLookupTable5 { fn from(point: &'a edwards::EdwardsPoint) -> Self { let A = ExtendedPoint::from(*point); let mut Ai = [CachedPoint::from(A); 8]; @@ -303,7 +303,7 @@ impl<'a> From<&'a edwards::EdwardsPoint> for OddLookupTable { Ai[i + 1] = (&A2 + &Ai[i]).into(); } // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A] - OddLookupTable(Ai) + NafLookupTable5(Ai) } } diff --git a/src/backend/avx2/scalar_mul/vartime_double_base.rs b/src/backend/avx2/scalar_mul/vartime_double_base.rs index aa2aa97..17d21f2 100644 --- a/src/backend/avx2/scalar_mul/vartime_double_base.rs +++ b/src/backend/avx2/scalar_mul/vartime_double_base.rs @@ -12,7 +12,7 @@ use traits::Identity; use scalar::Scalar; use edwards::EdwardsPoint; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; use backend::avx2::edwards::{CachedPoint, ExtendedPoint}; use backend::avx2::constants::BASEPOINT_ODD_LOOKUP_TABLE; @@ -30,7 +30,7 @@ pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { } } - let table_A = OddLookupTable::::from(A); + let table_A = NafLookupTable5::::from(A); let table_B = &BASEPOINT_ODD_LOOKUP_TABLE; let mut Q = ExtendedPoint::identity(); diff --git a/src/backend/avx2/scalar_mul/vartime_straus.rs b/src/backend/avx2/scalar_mul/vartime_straus.rs index f05754b..30e68d0 100644 --- a/src/backend/avx2/scalar_mul/vartime_straus.rs +++ b/src/backend/avx2/scalar_mul/vartime_straus.rs @@ -14,7 +14,7 @@ use core::borrow::Borrow; use traits::Identity; use scalar::Scalar; use edwards::EdwardsPoint; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; use backend::avx2::edwards::{CachedPoint, ExtendedPoint}; /// Perform variable-time, variable-base scalar multiplication. @@ -31,7 +31,7 @@ where .collect(); let lookup_tables: Vec<_> = points .into_iter() - .map(|point| OddLookupTable::::from(point.borrow())) + .map(|point| NafLookupTable5::::from(point.borrow())) .collect(); let mut Q = ExtendedPoint::identity(); diff --git a/src/scalar_mul/vartime_double_base.rs b/src/scalar_mul/vartime_double_base.rs index 9fcd267..bdb7005 100644 --- a/src/scalar_mul/vartime_double_base.rs +++ b/src/scalar_mul/vartime_double_base.rs @@ -14,7 +14,7 @@ use traits::Identity; use scalar::Scalar; use edwards::EdwardsPoint; use curve_models::{ProjectiveNielsPoint, ProjectivePoint}; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; /// Compute \\(aA + bB\\) in variable time, where \\(B\\) is the Ed25519 basepoint. pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { @@ -30,7 +30,7 @@ pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { } } - let table_A = OddLookupTable::::from(A); + let table_A = NafLookupTable5::::from(A); let table_B = &constants::AFFINE_ODD_MULTIPLES_OF_BASEPOINT; let mut r = ProjectivePoint::identity(); diff --git a/src/scalar_mul/vartime_straus.rs b/src/scalar_mul/vartime_straus.rs index b419c0a..03c1f75 100644 --- a/src/scalar_mul/vartime_straus.rs +++ b/src/scalar_mul/vartime_straus.rs @@ -15,7 +15,7 @@ use traits::Identity; use scalar::Scalar; use edwards::EdwardsPoint; use curve_models::{CompletedPoint, ProjectivePoint, ProjectiveNielsPoint}; -use scalar_mul::window::OddLookupTable; +use scalar_mul::window::NafLookupTable5; /// Perform variable-time, variable-base scalar multiplication. pub(crate) fn multiscalar_mul(scalars: I, points: J) -> EdwardsPoint @@ -31,7 +31,7 @@ where .collect(); let lookup_tables: Vec<_> = points .into_iter() - .map(|P| OddLookupTable::::from(P.borrow())) + .map(|P| NafLookupTable5::::from(P.borrow())) .collect(); let mut r = ProjectivePoint::identity(); diff --git a/src/scalar_mul/window.rs b/src/scalar_mul/window.rs index 9b25917..e0b69c2 100644 --- a/src/scalar_mul/window.rs +++ b/src/scalar_mul/window.rs @@ -122,9 +122,9 @@ impl<'a> From<&'a EdwardsPoint> for LookupTable { /// Holds odd multiples 1A, 3A, ..., 15A of a point A. #[derive(Copy, Clone)] -pub(crate) struct OddLookupTable(pub(crate) [T; 8]); +pub(crate) struct NafLookupTable5(pub(crate) [T; 8]); -impl OddLookupTable { +impl NafLookupTable5 { /// Given public, odd \\( x \\) with \\( 0 < x < 2^4 \\), return \\(xA\\). pub fn select(&self, x: usize) -> T { debug_assert_eq!(x & 1, 1); @@ -134,13 +134,13 @@ impl OddLookupTable { } } -impl Debug for OddLookupTable { +impl Debug for NafLookupTable5 { fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { - write!(f, "OddLookupTable({:?})", self.0) + write!(f, "NafLookupTable5({:?})", self.0) } } -impl<'a> From<&'a EdwardsPoint> for OddLookupTable { +impl<'a> From<&'a EdwardsPoint> for NafLookupTable5 { fn from(A: &'a EdwardsPoint) -> Self { let mut Ai = [A.to_projective_niels(); 8]; let A2 = A.double(); @@ -148,11 +148,11 @@ impl<'a> From<&'a EdwardsPoint> for OddLookupTable { Ai[i + 1] = (&A2 + &Ai[i]).to_extended().to_projective_niels(); } // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A] - OddLookupTable(Ai) + NafLookupTable5(Ai) } } -impl<'a> From<&'a EdwardsPoint> for OddLookupTable { +impl<'a> From<&'a EdwardsPoint> for NafLookupTable5 { fn from(A: &'a EdwardsPoint) -> Self { let mut Ai = [A.to_affine_niels(); 8]; let A2 = A.double(); @@ -160,6 +160,6 @@ impl<'a> From<&'a EdwardsPoint> for OddLookupTable { Ai[i + 1] = (&A2 + &Ai[i]).to_extended().to_affine_niels(); } // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A] - OddLookupTable(Ai) + NafLookupTable5(Ai) } } From 62d43752df7fbfd14ad775d5e605536e9cfaa0ba Mon Sep 17 00:00:00 2001 From: Henry & Isis Date: Thu, 5 Apr 2018 05:12:13 +0000 Subject: [PATCH 3/9] Add NafLookupTable8 and use for pre-computed basepoint table generation. --- build.rs | 10 +++--- src/scalar_mul/vartime_double_base.rs | 2 +- src/scalar_mul/window.rs | 47 +++++++++++++++++++++++++++ 3 files changed, 53 insertions(+), 6 deletions(-) diff --git a/build.rs b/build.rs index 7d06199..389e7bb 100644 --- a/build.rs +++ b/build.rs @@ -59,7 +59,7 @@ mod scalar_mul; use edwards::EdwardsBasepointTable; use curve_models::AffineNielsPoint; -use scalar_mul::window::NafLookupTable5; +use scalar_mul::window::NafLookupTable8; fn main() { // Enable the "precomputed_tables" feature in the main build stage @@ -86,7 +86,7 @@ use edwards::EdwardsBasepointTable; use curve_models::AffineNielsPoint; use scalar_mul::window::LookupTable; -use scalar_mul::window::NafLookupTable5; +use scalar_mul::window::NafLookupTable8; /// Table containing precomputed multiples of the Ed25519 basepoint \\\\(B = (x, 4/5)\\\\). pub const ED25519_BASEPOINT_TABLE: EdwardsBasepointTable = ED25519_BASEPOINT_TABLE_INNER_DOC_HIDDEN; @@ -101,13 +101,13 @@ pub const ED25519_BASEPOINT_TABLE_INNER_DOC_HIDDEN: EdwardsBasepointTable = {:?} // Now generate AFFINE_ODD_MULTIPLES_OF_BASEPOINT let B = &constants::ED25519_BASEPOINT_POINT; - let odd_multiples = NafLookupTable5::::from(B); + let odd_multiples = NafLookupTable8::::from(B); f.write_all( format!( "\n -/// Odd multiples of the basepoint `[B, 3B, 5B, 7B, 9B, 11B, 13B, 15B]`. -pub(crate) const AFFINE_ODD_MULTIPLES_OF_BASEPOINT: NafLookupTable5 = {:?}; +/// Odd multiples of the basepoint `[B, 3B, 5B, 7B, 9B, 11B, 13B, 15B, ..., 127B]`. +pub(crate) const AFFINE_ODD_MULTIPLES_OF_BASEPOINT: NafLookupTable8 = {:?}; \n\n", &odd_multiples ).as_bytes(), diff --git a/src/scalar_mul/vartime_double_base.rs b/src/scalar_mul/vartime_double_base.rs index bdb7005..86022d2 100644 --- a/src/scalar_mul/vartime_double_base.rs +++ b/src/scalar_mul/vartime_double_base.rs @@ -19,7 +19,7 @@ use scalar_mul::window::NafLookupTable5; /// Compute \\(aA + bB\\) in variable time, where \\(B\\) is the Ed25519 basepoint. pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { let a_naf = a.non_adjacent_form(5); - let b_naf = b.non_adjacent_form(5); + let b_naf = b.non_adjacent_form(8); // Find starting index let mut i: usize = 255; diff --git a/src/scalar_mul/window.rs b/src/scalar_mul/window.rs index e0b69c2..91ebb65 100644 --- a/src/scalar_mul/window.rs +++ b/src/scalar_mul/window.rs @@ -163,3 +163,50 @@ impl<'a> From<&'a EdwardsPoint> for NafLookupTable5 { NafLookupTable5(Ai) } } + +/// Holds stuff up to 8. +#[derive(Copy, Clone)] +pub(crate) struct NafLookupTable8(pub(crate) [T; 64]); + +impl NafLookupTable8 { + pub fn select(&self, x: usize) -> T { + debug_assert_eq!(x & 1, 1); + debug_assert!(x < 256); + + self.0[x / 2] + } +} + +impl Debug for NafLookupTable8 { + fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { + write!(f, "NafLookupTable8([\n")?; + for i in 0..64 { + write!(f, "\t{:?},\n", &self.0[i])?; + } + write!(f, "])") + } +} + +impl<'a> From<&'a EdwardsPoint> for NafLookupTable8 { + fn from(A: &'a EdwardsPoint) -> Self { + let mut Ai = [A.to_projective_niels(); 64]; + let A2 = A.double(); + for i in 0..63 { + Ai[i + 1] = (&A2 + &Ai[i]).to_extended().to_projective_niels(); + } + // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A, ..., 127A] + NafLookupTable8(Ai) + } +} + +impl<'a> From<&'a EdwardsPoint> for NafLookupTable8 { + fn from(A: &'a EdwardsPoint) -> Self { + let mut Ai = [A.to_affine_niels(); 64]; + let A2 = A.double(); + for i in 0..63 { + Ai[i + 1] = (&A2 + &Ai[i]).to_extended().to_affine_niels(); + } + // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A, ..., 127A] + NafLookupTable8(Ai) + } +} From 6b768c2a1ae981e9f984f40664a307e2b94c0fca Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 5 Apr 2018 16:20:47 -0700 Subject: [PATCH 4/9] Change AVX2 backend to use width-8 tables --- src/backend/avx2/constants.rs | 3450 ++++++++++++++++- src/backend/avx2/edwards.rs | 15 +- .../avx2/scalar_mul/vartime_double_base.rs | 2 +- 3 files changed, 3392 insertions(+), 75 deletions(-) diff --git a/src/backend/avx2/constants.rs b/src/backend/avx2/constants.rs index 0ef2300..e99e312 100644 --- a/src/backend/avx2/constants.rs +++ b/src/backend/avx2/constants.rs @@ -12,102 +12,3406 @@ use core::simd::u32x8; -use scalar_mul::window::NafLookupTable5; +use backend::avx2::edwards::{CachedPoint, ExtendedPoint}; use backend::avx2::field::FieldElement32x4; -use backend::avx2::edwards::{ExtendedPoint, CachedPoint}; +use scalar_mul::window::NafLookupTable8; /// 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] /// ``` -pub(crate) static P_TIMES_2_LO: u32x8 = - u32x8::new(67108845 << 1, 67108845 << 1, 33554431 << 1, 33554431 << 1, 67108845 << 1, 67108845 << 1, 33554431 << 1, 33554431 << 1); +pub(crate) static P_TIMES_2_LO: u32x8 = u32x8::new( + 67108845 << 1, + 67108845 << 1, + 33554431 << 1, + 33554431 << 1, + 67108845 << 1, + 67108845 << 1, + 33554431 << 1, + 33554431 << 1, +); /// The high 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] /// ``` -pub(crate) static P_TIMES_2_HI: u32x8 = - u32x8::new(67108863 << 1, 67108863 << 1, 33554431 << 1, 33554431 << 1, 67108863 << 1, 67108863 << 1, 33554431 << 1, 33554431 << 1); +pub(crate) static P_TIMES_2_HI: u32x8 = u32x8::new( + 67108863 << 1, + 67108863 << 1, + 33554431 << 1, + 33554431 << 1, + 67108863 << 1, + 67108863 << 1, + 33554431 << 1, + 33554431 << 1, +); /// The low limbs of (16p, 16p, 16p, 16p), so that /// ```no_run /// (16p, 16p, 16p, 16p) = [P_TIMES_16_LO, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI] /// ``` -pub(crate) static P_TIMES_16_LO: u32x8 = - u32x8::new(67108845 << 4, 67108845 << 4, 33554431 << 4, 33554431 << 4, 67108845 << 4, 67108845 << 4, 33554431 << 4, 33554431 << 4); +pub(crate) static P_TIMES_16_LO: u32x8 = u32x8::new( + 67108845 << 4, + 67108845 << 4, + 33554431 << 4, + 33554431 << 4, + 67108845 << 4, + 67108845 << 4, + 33554431 << 4, + 33554431 << 4, +); /// The high limbs of (16p, 16p, 16p, 16p), so that /// ```no_run /// (16p, 16p, 16p, 16p) = [P_TIMES_16_LO, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI] /// ``` -pub(crate) static P_TIMES_16_HI: u32x8 = - u32x8::new(67108863 << 4, 67108863 << 4, 33554431 << 4, 33554431 << 4, 67108863 << 4, 67108863 << 4, 33554431 << 4, 33554431 << 4); +pub(crate) static P_TIMES_16_HI: u32x8 = u32x8::new( + 67108863 << 4, + 67108863 << 4, + 33554431 << 4, + 33554431 << 4, + 67108863 << 4, + 67108863 << 4, + 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) + 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: NafLookupTable5 = - NafLookupTable5([ - CachedPoint(FieldElement32x4([ - u32x8::new(3571425, 10045002, 19036563, 1096096, 243332, 65897020, 0, 28963681), - u32x8::new(30896895, 63055514, 1614915, 5095970, 0, 53791688, 0, 31258312), - u32x8::new(13347627, 40339464, 2236269, 11185503, 0, 22520087, 0, 8659512), - u32x8::new(11125413, 29139905, 32037254, 28360723, 0, 64556417, 0, 9635759), - u32x8::new(33268144, 47262491, 4336918, 15795740, 0, 22027545, 0, 4846528), - ])), - CachedPoint(FieldElement32x4([ - u32x8::new(47099681, 31447946, 29365447, 24740513, 42991046, 18317844, 16051644, 21404226), - u32x8::new(31708133, 28909527, 2366091, 13703791, 469246, 54159622, 2601402, 32988002), - u32x8::new(63432457, 30251794, 15163516, 18491340, 28144087, 35605455, 13682295, 18474872), - u32x8::new(12221607, 4967598, 26061980, 26008006, 20226147, 9726961, 17410, 18051083), - u32x8::new(60569645, 62487085, 11911242, 21920922, 4092105, 38186967, 22431483, 31366585), - ])), - CachedPoint(FieldElement32x4([ - u32x8::new(18147205, 62587998, 2554617, 536692, 11924528, 26674131, 17645433, 24341419), - u32x8::new(11573357, 27579485, 31491870, 29000885, 10800976, 51902791, 28076395, 20464029), - u32x8::new(56031649, 10856669, 11791193, 26769430, 25306956, 5922200, 6630685, 9385098), - u32x8::new(31319348, 23906711, 16290213, 32142166, 61106354, 17181823, 3548308, 12022566), - u32x8::new(5904298, 50218605, 11826440, 5492249, 10379071, 3472255, 172742, 31948344), - ])), - CachedPoint(FieldElement32x4([ - u32x8::new(10625852, 15193821, 22918394, 23676410, 53695416, 54987793, 10067515, 11747680), - u32x8::new(65013325, 1309652, 29616320, 28922974, 60360891, 19621771, 9938982, 30406429), - u32x8::new(54967954, 65931918, 5595602, 25719523, 64909864, 30566415, 15945272, 8495317), - u32x8::new(1167157, 55265018, 11507029, 31641054, 43497904, 2367338, 12937761, 27517066), - u32x8::new(656704, 2544994, 13006713, 480979, 38471594, 62541240, 25353597, 11531760), - ])), - CachedPoint(FieldElement32x4([ - u32x8::new(22176662, 3984313, 27495285, 4110608, 2909584, 30594106, 15677919, 2549183), - u32x8::new(33979105, 62269905, 2071511, 6894756, 53189950, 47232857, 6408191, 6123225), - u32x8::new(32553873, 63948030, 12612401, 3633166, 24054373, 37626618, 14481327, 8520484), - u32x8::new(56552486, 10749438, 12034813, 28811946, 1445640, 36755601, 12104575, 10257833), - u32x8::new(22795808, 48761311, 1136056, 9380768, 1411523, 5341811, 27318329, 9686767)])), - CachedPoint(FieldElement32x4([ - u32x8::new(21157200, 39156966, 20473176, 4934657, 61478183, 45121537, 5429856, 13035023), - u32x8::new(7954529, 58789246, 31440083, 7054221, 38438565, 36856107, 1364112, 14548122), - u32x8::new(26120083, 36321360, 4919997, 31687496, 33757765, 36237559, 15243054, 32163861), - u32x8::new(25878307, 46544824, 19455951, 2414935, 16844726, 56521560, 32680554, 26660660), - u32x8::new(48360220, 43407178, 12187042, 24925816, 7423722, 25746484, 12814654, 17395963), - ])), - CachedPoint(FieldElement32x4([ - u32x8::new(63153652, 32195955, 4087908, 8431689, 30392384, 47203165, 8986649, 9053039), - u32x8::new(63659241, 47988767, 2931872, 19953600, 11747107, 51610101, 20952181, 13364887), - u32x8::new(3659197, 58790649, 5930099, 2605312, 28477896, 580728, 20579735, 2610622), - u32x8::new(41781607, 17161358, 10690531, 24368015, 47027031, 36742339, 5414694, 13156365), - u32x8::new(13237853, 51182423, 8954802, 29006542, 22643989, 56896541, 22830593, 10289708), - ])), - CachedPoint(FieldElement32x4([ - u32x8::new(1401265, 58846825, 30911620, 32239180, 15391552, 15200821, 6339309, 16403588), - u32x8::new(55913797, 29541724, 1664461, 21709410, 38470488, 47097092, 17674945, 32666066), - u32x8::new(22844482, 10797709, 27548106, 31638735, 34500968, 26611503, 19727211, 13160873), - u32x8::new(31485204, 14496164, 13981208, 10276888, 5748808, 35024436, 2740987, 7479021), - u32x8::new(58541207, 14866135, 32344041, 545930, 62661488, 6941250, 27940205, 11976112), - ])), - ]); +pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: NafLookupTable8 = NafLookupTable8([ + CachedPoint(FieldElement32x4([ + u32x8::new( + 3571425, + 10045002, + 19036563, + 1096096, + 243332, + 65897020, + 0, + 28963681, + ), + u32x8::new( + 30896895, + 63055514, + 1614915, + 5095970, + 0, + 53791688, + 0, + 31258312, + ), + u32x8::new( + 13347627, + 40339464, + 2236269, + 11185503, + 0, + 22520087, + 0, + 8659512, + ), + u32x8::new( + 11125413, + 29139905, + 32037254, + 28360723, + 0, + 64556417, + 0, + 9635759, + ), + u32x8::new( + 33268144, + 47262491, + 4336918, + 15795740, + 0, + 22027545, + 0, + 4846528, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 47099681, + 31447946, + 29365447, + 24740513, + 42991046, + 18317844, + 16051644, + 21404226, + ), + u32x8::new( + 31708133, + 28909527, + 2366091, + 13703791, + 469246, + 54159622, + 2601402, + 32988002, + ), + u32x8::new( + 63432457, + 30251794, + 15163516, + 18491340, + 28144087, + 35605455, + 13682295, + 18474872, + ), + u32x8::new( + 12221607, + 4967598, + 26061980, + 26008006, + 20226147, + 9726961, + 17410, + 18051083, + ), + u32x8::new( + 60569645, + 62487085, + 11911242, + 21920922, + 4092105, + 38186967, + 22431483, + 31366585, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 18147205, + 62587998, + 2554617, + 536692, + 11924528, + 26674131, + 17645433, + 24341419, + ), + u32x8::new( + 11573357, + 27579485, + 31491870, + 29000885, + 10800976, + 51902791, + 28076395, + 20464029, + ), + u32x8::new( + 56031649, + 10856669, + 11791193, + 26769430, + 25306956, + 5922200, + 6630685, + 9385098, + ), + u32x8::new( + 31319348, + 23906711, + 16290213, + 32142166, + 61106354, + 17181823, + 3548308, + 12022566, + ), + u32x8::new( + 5904298, + 50218605, + 11826440, + 5492249, + 10379071, + 3472255, + 172742, + 31948344, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 10625852, + 15193821, + 22918394, + 23676410, + 53695416, + 54987793, + 10067515, + 11747680, + ), + u32x8::new( + 65013325, + 1309652, + 29616320, + 28922974, + 60360891, + 19621771, + 9938982, + 30406429, + ), + u32x8::new( + 54967954, + 65931918, + 5595602, + 25719523, + 64909864, + 30566415, + 15945272, + 8495317, + ), + u32x8::new( + 1167157, + 55265018, + 11507029, + 31641054, + 43497904, + 2367338, + 12937761, + 27517066, + ), + u32x8::new( + 656704, + 2544994, + 13006713, + 480979, + 38471594, + 62541240, + 25353597, + 11531760, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 22176662, + 3984313, + 27495285, + 4110608, + 2909584, + 30594106, + 15677919, + 2549183, + ), + u32x8::new( + 33979105, + 62269905, + 2071511, + 6894756, + 53189950, + 47232857, + 6408191, + 6123225, + ), + u32x8::new( + 32553873, + 63948030, + 12612401, + 3633166, + 24054373, + 37626618, + 14481327, + 8520484, + ), + u32x8::new( + 56552486, + 10749438, + 12034813, + 28811946, + 1445640, + 36755601, + 12104575, + 10257833, + ), + u32x8::new( + 22795808, + 48761311, + 1136056, + 9380768, + 1411523, + 5341811, + 27318329, + 9686767, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 21157200, + 39156966, + 20473176, + 4934657, + 61478183, + 45121537, + 5429856, + 13035023, + ), + u32x8::new( + 7954529, + 58789246, + 31440083, + 7054221, + 38438565, + 36856107, + 1364112, + 14548122, + ), + u32x8::new( + 26120083, + 36321360, + 4919997, + 31687496, + 33757765, + 36237559, + 15243054, + 32163861, + ), + u32x8::new( + 25878307, + 46544824, + 19455951, + 2414935, + 16844726, + 56521560, + 32680554, + 26660660, + ), + u32x8::new( + 48360220, + 43407178, + 12187042, + 24925816, + 7423722, + 25746484, + 12814654, + 17395963, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 63153652, + 32195955, + 4087908, + 8431689, + 30392384, + 47203165, + 8986649, + 9053039, + ), + u32x8::new( + 63659241, + 47988767, + 2931872, + 19953600, + 11747107, + 51610101, + 20952181, + 13364887, + ), + u32x8::new( + 3659197, + 58790649, + 5930099, + 2605312, + 28477896, + 580728, + 20579735, + 2610622, + ), + u32x8::new( + 41781607, + 17161358, + 10690531, + 24368015, + 47027031, + 36742339, + 5414694, + 13156365, + ), + u32x8::new( + 13237853, + 51182423, + 8954802, + 29006542, + 22643989, + 56896541, + 22830593, + 10289708, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 1401265, + 58846825, + 30911620, + 32239180, + 15391552, + 15200821, + 6339309, + 16403588, + ), + u32x8::new( + 55913797, + 29541724, + 1664461, + 21709410, + 38470488, + 47097092, + 17674945, + 32666066, + ), + u32x8::new( + 22844482, + 10797709, + 27548106, + 31638735, + 34500968, + 26611503, + 19727211, + 13160873, + ), + u32x8::new( + 31485204, + 14496164, + 13981208, + 10276888, + 5748808, + 35024436, + 2740987, + 7479021, + ), + u32x8::new( + 58541207, + 14866135, + 32344041, + 545930, + 62661488, + 6941250, + 27940205, + 11976112, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 39849808, + 44781685, + 15697329, + 24387845, + 12501486, + 50260092, + 23199481, + 31929024, + ), + u32x8::new( + 24823070, + 27956017, + 27034296, + 10316465, + 47664045, + 11152446, + 15719183, + 30181617, + ), + u32x8::new( + 20771189, + 19969144, + 31433937, + 19185213, + 27565920, + 10384445, + 2893359, + 9255362, + ), + u32x8::new( + 42894974, + 11925545, + 32134441, + 32738810, + 55916336, + 32479272, + 19563550, + 5511385, + ), + u32x8::new( + 17857161, + 47809169, + 14564114, + 27997751, + 33024640, + 38669671, + 31956536, + 27313245, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 58237774, + 15917425, + 18872208, + 19394230, + 17374297, + 6101419, + 4839741, + 6596900, + ), + u32x8::new( + 66947393, + 15744215, + 18368993, + 17750160, + 41006525, + 9205497, + 2629667, + 32170865, + ), + u32x8::new( + 66481381, + 1919414, + 28338762, + 7372967, + 33819153, + 4156199, + 27126309, + 12739816, + ), + u32x8::new( + 44117158, + 58545296, + 22521371, + 11809712, + 28998792, + 50731010, + 30215699, + 25748377, + ), + u32x8::new( + 23561284, + 4160244, + 9035405, + 24895184, + 39761639, + 59253416, + 8684759, + 22487864, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 12671134, + 56419053, + 16092401, + 30038207, + 4002647, + 47822606, + 7151311, + 28430768, + ), + u32x8::new( + 61041684, + 35765374, + 30598048, + 19666539, + 44150175, + 40140037, + 290469, + 28442674, + ), + u32x8::new( + 18847796, + 1371617, + 33316881, + 13199936, + 43646578, + 17068881, + 12074900, + 1537415, + ), + u32x8::new( + 10052225, + 38316070, + 27469797, + 5297537, + 50725570, + 20435349, + 10339121, + 2779737, + ), + u32x8::new( + 18372189, + 15466385, + 24762130, + 22217964, + 23503887, + 47844464, + 10415034, + 2606889, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 55082775, + 45300503, + 16032654, + 5964396, + 17743504, + 24634761, + 19493066, + 5184611, + ), + u32x8::new( + 50172633, + 35093294, + 10040575, + 23616256, + 4543900, + 61852191, + 4049821, + 7423669, + ), + u32x8::new( + 20295398, + 40009376, + 10487190, + 15670429, + 51972856, + 58649552, + 20436392, + 3432497, + ), + u32x8::new( + 35189420, + 54117751, + 12825868, + 6283038, + 27540739, + 30648758, + 22658912, + 9466689, + ), + u32x8::new( + 51737549, + 40725785, + 17409814, + 25201086, + 21156239, + 34176168, + 26814520, + 5956424, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 8211442, + 8014184, + 6260823, + 22108096, + 32182620, + 51844847, + 2466270, + 28582231, + ), + u32x8::new( + 27199739, + 3848333, + 31738017, + 10892045, + 4963982, + 65391770, + 32551997, + 28906469, + ), + u32x8::new( + 16606846, + 32207068, + 26404535, + 7614129, + 45416902, + 65584718, + 13821785, + 2646060, + ), + u32x8::new( + 36090634, + 57981287, + 32247670, + 22837502, + 31003861, + 55448117, + 6062915, + 20369975, + ), + u32x8::new( + 27381403, + 50578107, + 522631, + 29521058, + 31137497, + 40220737, + 27628049, + 1824195, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 59402443, + 17056879, + 29262689, + 6131785, + 52551472, + 43367471, + 29423199, + 18899208, + ), + u32x8::new( + 5749414, + 43514612, + 11365899, + 21514624, + 65591890, + 60945892, + 19841732, + 5628567, + ), + u32x8::new( + 19334369, + 52500268, + 12307673, + 5267367, + 3212103, + 9035822, + 29142161, + 30520954, + ), + u32x8::new( + 57261330, + 6819646, + 22089161, + 9800373, + 55155453, + 62250856, + 13766735, + 25244545, + ), + u32x8::new( + 54370226, + 61888301, + 24496089, + 2540581, + 65637506, + 60274355, + 18154273, + 11687259, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 12521903, + 26014045, + 13995625, + 33360175, + 23605474, + 7376434, + 27229267, + 17195036, + ), + u32x8::new( + 59482891, + 10074423, + 574357, + 3857753, + 61377787, + 50306685, + 5241065, + 20234396, + ), + u32x8::new( + 23674717, + 6997172, + 20771841, + 16858511, + 40565304, + 29973136, + 7049812, + 14585010, + ), + u32x8::new( + 1427477, + 13295732, + 31762066, + 31499740, + 60419925, + 54666164, + 22009424, + 8089609, + ), + u32x8::new( + 58154031, + 41593020, + 15342328, + 957047, + 38937260, + 37037498, + 24871992, + 32973409, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 30654745, + 51286025, + 21206982, + 2433562, + 12780105, + 31732574, + 33087964, + 33081189, + ), + u32x8::new( + 66640017, + 42720009, + 16567620, + 15300745, + 1530367, + 33001123, + 20930247, + 21042661, + ), + u32x8::new( + 15003356, + 5294119, + 22985605, + 18928772, + 32628461, + 18230172, + 14773298, + 27193722, + ), + u32x8::new( + 27555, + 65346287, + 17017174, + 7837720, + 21499787, + 42855613, + 22474984, + 13675085, + ), + u32x8::new( + 24164369, + 50130116, + 5973149, + 24152073, + 1577334, + 25400030, + 18648484, + 32228854, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 49518649, + 59119280, + 31670678, + 20396561, + 61728330, + 651402, + 176032, + 9529498, + ), + u32x8::new( + 61765532, + 9082232, + 32794568, + 15526956, + 48543100, + 32614212, + 19001206, + 25680229, + ), + u32x8::new( + 32086091, + 10373081, + 8996131, + 31822823, + 35788988, + 49973190, + 30542040, + 17858455, + ), + u32x8::new( + 48130197, + 58121889, + 27753291, + 29923268, + 54448075, + 43300790, + 9336565, + 15770022, + ), + u32x8::new( + 57725546, + 20557498, + 9366233, + 16023566, + 16189031, + 2837363, + 24315301, + 27003505, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 28286608, + 10767548, + 18220739, + 5413236, + 48253387, + 58255702, + 11864864, + 28527159, + ), + u32x8::new( + 45038176, + 58655197, + 25648758, + 10951484, + 42564382, + 34542843, + 23146954, + 22234334, + ), + u32x8::new( + 14858710, + 24978793, + 15040559, + 4379220, + 47621477, + 40271440, + 15650420, + 1998736, + ), + u32x8::new( + 24106391, + 9626149, + 344505, + 25253814, + 34579800, + 59687089, + 25718289, + 25904133, + ), + u32x8::new( + 1981195, + 37751302, + 26132048, + 1764722, + 13288231, + 28808622, + 12531301, + 18292949, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 13869851, + 31448904, + 14963539, + 7581293, + 20536485, + 35021083, + 21257574, + 33356609, + ), + u32x8::new( + 36903364, + 18429241, + 11097857, + 5943856, + 60583077, + 40015815, + 30509523, + 31915271, + ), + u32x8::new( + 49161801, + 40681915, + 67892, + 25454357, + 22779677, + 25798439, + 15964829, + 5863227, + ), + u32x8::new( + 60810637, + 4496471, + 5217137, + 14095116, + 50942411, + 50712663, + 2507380, + 26844507, + ), + u32x8::new( + 34579752, + 53519385, + 10859797, + 18816024, + 42552864, + 39478521, + 6783896, + 17277037, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 43287109, + 27900723, + 33182187, + 2766754, + 17041989, + 1018260, + 33392790, + 4830032, + ), + u32x8::new( + 60194178, + 30788903, + 24728888, + 14513195, + 20897010, + 28843233, + 20111980, + 17475240, + ), + u32x8::new( + 46042274, + 19257042, + 4628173, + 31649727, + 27388316, + 66631493, + 11541886, + 6408028, + ), + u32x8::new( + 57024680, + 49536568, + 32050358, + 31321917, + 17437691, + 49672356, + 2884755, + 20493991, + ), + u32x8::new( + 59553007, + 46782643, + 29001173, + 1814088, + 21930692, + 51319706, + 14965872, + 30748046, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 16441817, + 36111849, + 6900424, + 602234, + 46522199, + 16441484, + 8135070, + 21726541, + ), + u32x8::new( + 37711225, + 32701959, + 11679112, + 13125533, + 32154135, + 9407918, + 26554289, + 620848, + ), + u32x8::new( + 19233407, + 30086864, + 14679568, + 2797374, + 4892806, + 7993077, + 247658, + 5632804, + ), + u32x8::new( + 37427262, + 26675495, + 27125659, + 13496131, + 50718473, + 40115609, + 28505351, + 27837393, + ), + u32x8::new( + 196819, + 18410429, + 7070012, + 21691388, + 29763371, + 24754123, + 9727048, + 10930179, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 28319289, + 40734650, + 16225680, + 24739184, + 64272368, + 35356897, + 7866648, + 13635853, + ), + u32x8::new( + 34165295, + 48328447, + 27041670, + 23643655, + 48949950, + 52963288, + 30411133, + 6045174, + ), + u32x8::new( + 18583559, + 41649834, + 9813585, + 26098520, + 25682734, + 26733526, + 19276490, + 10654728, + ), + u32x8::new( + 34867476, + 52715968, + 5694571, + 13380978, + 15134994, + 1831255, + 8608001, + 17266401, + ), + u32x8::new( + 59925903, + 44282172, + 27802465, + 1855069, + 14234749, + 36635487, + 11302294, + 10938429, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 8373273, + 49064494, + 4932071, + 32997499, + 38472880, + 29335908, + 14504412, + 22460029, + ), + u32x8::new( + 31795930, + 50785923, + 25835990, + 25790073, + 65669841, + 11360450, + 9969157, + 9008164, + ), + u32x8::new( + 50262498, + 45869261, + 16124434, + 15336007, + 882762, + 42522623, + 11277198, + 26296377, + ), + u32x8::new( + 42332732, + 59129236, + 14452816, + 567985, + 208061, + 34722729, + 32008143, + 14828749, + ), + u32x8::new( + 17937794, + 36846032, + 32102665, + 4442466, + 19745435, + 31633451, + 7146411, + 15812027, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 30741269, + 38648744, + 12562645, + 30092623, + 25073992, + 28730659, + 27911745, + 30000958, + ), + u32x8::new( + 2859794, + 25991700, + 17776078, + 27091930, + 2328322, + 60061146, + 18581824, + 18039008, + ), + u32x8::new( + 58206333, + 17917354, + 1972306, + 11853766, + 2655376, + 60543390, + 18416710, + 13287440, + ), + u32x8::new( + 62746330, + 61423885, + 21246577, + 2266675, + 60099139, + 14804707, + 14772234, + 20679434, + ), + u32x8::new( + 26987698, + 15488817, + 715616, + 2339565, + 51980752, + 17333865, + 21965103, + 10839820, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 18672548, + 57660959, + 16042910, + 19519287, + 62865851, + 17580961, + 26628347, + 23774759, + ), + u32x8::new( + 368070, + 3464471, + 25888304, + 30370559, + 52396053, + 45426828, + 28745251, + 9246829, + ), + u32x8::new( + 29090099, + 57950037, + 23104657, + 4903923, + 10987778, + 56163684, + 23621539, + 10332760, + ), + u32x8::new( + 53338235, + 44851161, + 21606845, + 31069622, + 4243630, + 34464392, + 11286454, + 5802022, + ), + u32x8::new( + 46710757, + 63389067, + 11642865, + 1980986, + 12967337, + 28162061, + 3854192, + 30432268, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 12179834, + 41005450, + 12809619, + 33525228, + 4624405, + 46957889, + 16968743, + 11827816, + ), + u32x8::new( + 51521162, + 12466775, + 31791271, + 15303651, + 49798465, + 62714504, + 6509600, + 12918560, + ), + u32x8::new( + 20445559, + 1756449, + 28848701, + 7920171, + 9835040, + 5900071, + 28757409, + 12376688, + ), + u32x8::new( + 18259496, + 14281012, + 21767026, + 10232236, + 20000226, + 12400540, + 4104902, + 23570543, + ), + u32x8::new( + 3687440, + 26546648, + 13328821, + 26841081, + 49822734, + 22334054, + 244496, + 24862543, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 59523541, + 62195428, + 3853227, + 13954801, + 12387708, + 47627615, + 27221350, + 17899572, + ), + u32x8::new( + 63193587, + 36343307, + 14595132, + 6880795, + 1364792, + 37648434, + 3259017, + 20536046, + ), + u32x8::new( + 30362834, + 10440372, + 9574624, + 11729232, + 63861613, + 21748389, + 5530846, + 2721586, + ), + u32x8::new( + 18339760, + 1550632, + 17170271, + 25732971, + 28459263, + 63142237, + 21642345, + 31557672, + ), + u32x8::new( + 10611282, + 5204623, + 18049257, + 214175, + 19432723, + 49809070, + 26010406, + 27449522, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 19770733, + 26478685, + 9464541, + 29158041, + 28604307, + 45196604, + 7586524, + 6641859, + ), + u32x8::new( + 65654484, + 52230498, + 30886612, + 19112823, + 47271809, + 38942611, + 16020035, + 10773481, + ), + u32x8::new( + 27464323, + 54451016, + 20646645, + 17732915, + 23008717, + 53626684, + 3253189, + 15614410, + ), + u32x8::new( + 52381752, + 40693008, + 7063024, + 28469981, + 51159478, + 44543211, + 19941777, + 5985451, + ), + u32x8::new( + 13553668, + 35524849, + 14788737, + 1883845, + 12385775, + 47958835, + 29135466, + 1776722, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 36719806, + 20827965, + 23175373, + 32996806, + 42041892, + 65708790, + 5467143, + 20884008, + ), + u32x8::new( + 43256281, + 40770646, + 17244063, + 31959819, + 64366384, + 43544617, + 25057754, + 12628720, + ), + u32x8::new( + 17337782, + 58472057, + 27906934, + 15305274, + 30292418, + 39284317, + 16946773, + 24806712, + ), + u32x8::new( + 6485126, + 32447403, + 16261486, + 13561940, + 49439635, + 10738368, + 16419889, + 8897231, + ), + u32x8::new( + 44812203, + 40122262, + 25496058, + 2759794, + 25295304, + 52178368, + 24154195, + 29334408, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 42307254, + 57217102, + 1088936, + 3832827, + 33905401, + 23130334, + 6958056, + 12622851, + ), + u32x8::new( + 3881189, + 14870059, + 19712830, + 6071598, + 38147944, + 60776394, + 3427938, + 13765703, + ), + u32x8::new( + 7666911, + 24227591, + 17077136, + 22967588, + 6874639, + 30915523, + 11451695, + 24292224, + ), + u32x8::new( + 13659529, + 31984463, + 28764736, + 20506164, + 64729627, + 49321636, + 28284636, + 25472371, + ), + u32x8::new( + 39360308, + 42281399, + 9446504, + 868960, + 49227724, + 21351115, + 30561851, + 11292096, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 7071115, + 46444090, + 5387916, + 15432877, + 27226682, + 41506862, + 2398278, + 3978240, + ), + u32x8::new( + 51009614, + 54216973, + 24368938, + 31392616, + 38456150, + 62313644, + 6729154, + 99724, + ), + u32x8::new( + 17474332, + 62857913, + 2619930, + 30659308, + 18268181, + 32809239, + 22826292, + 24561895, + ), + u32x8::new( + 38187020, + 67003092, + 14118280, + 16500577, + 18808560, + 64983716, + 25712929, + 32518261, + ), + u32x8::new( + 25735813, + 62284262, + 10824872, + 20558596, + 48149681, + 31162667, + 22608274, + 26285185, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 963440, + 63742255, + 10230323, + 25515008, + 32506414, + 6105697, + 25980317, + 24645129, + ), + u32x8::new( + 7162189, + 8101249, + 14679265, + 33443386, + 2002396, + 8541405, + 19442276, + 4795881, + ), + u32x8::new( + 8116694, + 51463069, + 4415528, + 25599140, + 55805721, + 39582709, + 6719436, + 30033839, + ), + u32x8::new( + 14468202, + 42181869, + 25188826, + 9639755, + 47546189, + 62711146, + 32762447, + 18338064, + ), + u32x8::new( + 33880058, + 32810909, + 8969931, + 13095238, + 38360605, + 40138517, + 9246134, + 4928058, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 63655588, + 17883670, + 9410246, + 26162761, + 5000571, + 7349225, + 23785252, + 32751089, + ), + u32x8::new( + 28568737, + 10733123, + 9342397, + 21570673, + 54096560, + 32467591, + 20494687, + 21511513, + ), + u32x8::new( + 47675157, + 47932807, + 29250946, + 15672208, + 59760469, + 9945465, + 14939287, + 18437405, + ), + u32x8::new( + 37985267, + 8609815, + 31573002, + 3373596, + 47828883, + 20834216, + 13248616, + 24154292, + ), + u32x8::new( + 5543543, + 29553242, + 3386453, + 30501150, + 25058089, + 15236571, + 8814395, + 32462955, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 39158670, + 15322548, + 20495103, + 3312736, + 14557171, + 12985179, + 8044741, + 3176899, + ), + u32x8::new( + 24673290, + 29693310, + 21412266, + 18324699, + 2154518, + 40329021, + 17500543, + 3954277, + ), + u32x8::new( + 36758685, + 38738957, + 165513, + 14691866, + 3070475, + 10424235, + 17096536, + 16896898, + ), + u32x8::new( + 59790459, + 43094586, + 8720681, + 10423589, + 1122030, + 31545615, + 4463786, + 31811293, + ), + u32x8::new( + 49778992, + 60881044, + 20509974, + 5832494, + 64155961, + 31483358, + 4511231, + 20307815, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 2863373, + 40876242, + 26865913, + 24067353, + 15726407, + 40919070, + 12953902, + 9931535, + ), + u32x8::new( + 60934877, + 42512204, + 21649141, + 21945190, + 52211954, + 60984193, + 7046207, + 5363493, + ), + u32x8::new( + 4205971, + 64068464, + 18197273, + 7327176, + 51527794, + 21166920, + 20669933, + 11828242, + ), + u32x8::new( + 59782815, + 49617225, + 15379924, + 457923, + 9320508, + 21498914, + 3242540, + 31563182, + ), + u32x8::new( + 27714753, + 8664670, + 3366162, + 26338598, + 56775518, + 25796006, + 13129151, + 21388876, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 59276548, + 49972346, + 16795002, + 33455915, + 48430097, + 53857205, + 18627071, + 32474471, + ), + u32x8::new( + 42160315, + 50705892, + 13530540, + 28012698, + 19833221, + 55886870, + 20191784, + 9644313, + ), + u32x8::new( + 20372416, + 28414713, + 24084234, + 31804096, + 33815377, + 36131001, + 17251241, + 18291088, + ), + u32x8::new( + 56234667, + 14920441, + 2033267, + 29572003, + 1724043, + 45519699, + 17873735, + 501988, + ), + u32x8::new( + 50031659, + 31517850, + 15697583, + 1016845, + 43104661, + 54769582, + 8008601, + 27257051, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 52951491, + 66542164, + 14853573, + 30444631, + 12045973, + 24321813, + 16545674, + 18160646, + ), + u32x8::new( + 60107911, + 1126003, + 5947677, + 19486116, + 41119984, + 30860440, + 7935395, + 13354438, + ), + u32x8::new( + 17841328, + 11063269, + 1664538, + 26687568, + 6268968, + 22280371, + 17275484, + 4523163, + ), + u32x8::new( + 15886041, + 56799482, + 15446552, + 21712778, + 1005290, + 17827215, + 4978741, + 6854882, + ), + u32x8::new( + 34319277, + 47731002, + 20321804, + 28544575, + 29591814, + 63376351, + 24754545, + 26001714, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 66783087, + 5234346, + 46102, + 8566476, + 19947339, + 20180418, + 25398238, + 3726678, + ), + u32x8::new( + 63890180, + 46380965, + 20674069, + 5366544, + 59661487, + 48406612, + 31533614, + 7071217, + ), + u32x8::new( + 13104676, + 1406631, + 24326736, + 19854367, + 61039528, + 11019904, + 31967425, + 19219275, + ), + u32x8::new( + 39003597, + 30143957, + 15351834, + 8639435, + 57309582, + 61436794, + 15830475, + 10090318, + ), + u32x8::new( + 45923044, + 6700175, + 99413, + 21263025, + 23762647, + 53905481, + 6063914, + 10065424, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 42822326, + 57678669, + 4052879, + 25452667, + 54049411, + 2373092, + 22337016, + 7701046, + ), + u32x8::new( + 44382355, + 43307377, + 16761537, + 30373573, + 49790216, + 23230748, + 25655306, + 10519391, + ), + u32x8::new( + 919475, + 59371245, + 1273450, + 25558666, + 9724711, + 8556709, + 25755845, + 10887647, + ), + u32x8::new( + 25465699, + 44651158, + 17658392, + 11257418, + 29735193, + 22885150, + 7094716, + 26828565, + ), + u32x8::new( + 48237389, + 47661599, + 27054393, + 7328070, + 27280193, + 65616691, + 23062005, + 4170709, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 26535281, + 60238317, + 30343788, + 25790743, + 37993933, + 24614372, + 9523840, + 10401918, + ), + u32x8::new( + 2783987, + 29468958, + 4697011, + 19804475, + 37246678, + 46797720, + 10261254, + 18942252, + ), + u32x8::new( + 58135580, + 60247753, + 25301938, + 6844561, + 20949454, + 39844754, + 4552026, + 919057, + ), + u32x8::new( + 6694071, + 44126261, + 32285330, + 31370180, + 24603698, + 53328179, + 13971149, + 5325636, + ), + u32x8::new( + 64879487, + 582094, + 17982081, + 19190425, + 24951286, + 26923842, + 29077174, + 33286062, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 54863941, + 67016431, + 1224043, + 23371240, + 62940074, + 52101083, + 13523637, + 30366406, + ), + u32x8::new( + 36324581, + 25407485, + 18258623, + 4698602, + 50300544, + 2658516, + 26300935, + 2611030, + ), + u32x8::new( + 27183975, + 21791014, + 18105064, + 9875199, + 58118912, + 54198635, + 6400311, + 14767984, + ), + u32x8::new( + 33918318, + 42937962, + 14809334, + 22136592, + 10636588, + 29082337, + 29829692, + 28549776, + ), + u32x8::new( + 61080905, + 854212, + 12202487, + 20004503, + 9256495, + 6903981, + 20567109, + 347423, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 41391822, + 34336880, + 22362564, + 14247996, + 12115604, + 41583344, + 7639288, + 28910945, + ), + u32x8::new( + 62066617, + 59758859, + 26665947, + 11614812, + 65737664, + 45704543, + 30324810, + 12868376, + ), + u32x8::new( + 17491771, + 43589814, + 9454919, + 26047850, + 52629282, + 39304244, + 3868968, + 19296062, + ), + u32x8::new( + 17826638, + 30413590, + 32534225, + 32741469, + 15012391, + 14365713, + 33039233, + 14791399, + ), + u32x8::new( + 64115596, + 59197067, + 32739005, + 23275744, + 32954320, + 22241406, + 20788442, + 4942942, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 31956192, + 59570132, + 2784352, + 4237732, + 47222312, + 4860927, + 18658867, + 15279314, + ), + u32x8::new( + 63240583, + 28160478, + 23524941, + 13390861, + 66437406, + 57718120, + 33345312, + 28896298, + ), + u32x8::new( + 39026193, + 46239965, + 21440243, + 25070488, + 64012383, + 60999016, + 16517060, + 29565907, + ), + u32x8::new( + 18118181, + 60161496, + 4212092, + 23976240, + 36277753, + 62363144, + 5816868, + 16964362, + ), + u32x8::new( + 18196138, + 62490693, + 281468, + 7934713, + 56027312, + 62015725, + 4837237, + 32932252, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 29885826, + 51028067, + 30418143, + 33438769, + 62542283, + 39442528, + 31535876, + 143299, + ), + u32x8::new( + 17143063, + 56709783, + 14451852, + 15782104, + 32762665, + 14047066, + 26295037, + 5432487, + ), + u32x8::new( + 75151, + 533606, + 7539077, + 30926189, + 38410914, + 23771680, + 4872443, + 29199566, + ), + u32x8::new( + 61522396, + 48934708, + 16223126, + 207380, + 11171993, + 47975147, + 14164574, + 352966, + ), + u32x8::new( + 15449006, + 56530757, + 26796528, + 12045834, + 63738697, + 40667227, + 33001582, + 9101885, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 43331297, + 18431341, + 25801195, + 17267698, + 19365485, + 57295202, + 22218985, + 21284590, + ), + u32x8::new( + 2429849, + 19152559, + 10762172, + 22564684, + 21880390, + 66866426, + 20357935, + 22641906, + ), + u32x8::new( + 19771185, + 31652693, + 3666117, + 28136958, + 23624283, + 55101502, + 6313920, + 6783662, + ), + u32x8::new( + 3487137, + 7092443, + 11001876, + 26196524, + 47319246, + 44542068, + 17594073, + 15027760, + ), + u32x8::new( + 49563607, + 32191113, + 4991283, + 25400512, + 46539152, + 4155103, + 32368171, + 201203, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 20548943, + 14334571, + 4073874, + 6368588, + 53208883, + 56484515, + 15970071, + 25561889, + ), + u32x8::new( + 49915097, + 44030795, + 11202344, + 29284344, + 60258023, + 66225712, + 8075764, + 12383512, + ), + u32x8::new( + 45248912, + 4933668, + 9592153, + 5819559, + 31030983, + 38174071, + 32435814, + 7442522, + ), + u32x8::new( + 62688129, + 48218381, + 22089545, + 12897361, + 21050881, + 34278889, + 7569163, + 3225449, + ), + u32x8::new( + 19050183, + 51089071, + 32935757, + 22640195, + 66122318, + 47144608, + 18743677, + 25177079, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 41186817, + 46681702, + 31819867, + 32997133, + 38559207, + 27147015, + 30293819, + 16762988, + ), + u32x8::new( + 24154689, + 51762873, + 23883879, + 13510519, + 55338250, + 61224161, + 11663149, + 30803960, + ), + u32x8::new( + 18104238, + 14117824, + 11724021, + 21362053, + 65704761, + 35530242, + 13498058, + 33522849, + ), + u32x8::new( + 63812888, + 23995539, + 28920539, + 24005193, + 26412223, + 36582218, + 4251418, + 26160309, + ), + u32x8::new( + 16822053, + 66064082, + 3482145, + 31979593, + 45937188, + 54475379, + 612917, + 7976478, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 46509314, + 55327128, + 8944536, + 274914, + 26432930, + 53829300, + 21192572, + 3569894, + ), + u32x8::new( + 20919764, + 64356651, + 30642344, + 17215170, + 20335124, + 11203745, + 18663316, + 19024174, + ), + u32x8::new( + 59297055, + 53842463, + 3680204, + 9806710, + 54004169, + 51484914, + 29807998, + 20134199, + ), + u32x8::new( + 14781592, + 22628010, + 26877930, + 25880359, + 30434803, + 190607, + 30184292, + 8991040, + ), + u32x8::new( + 64400983, + 64591751, + 854562, + 28216111, + 20010398, + 50414793, + 9803872, + 22687008, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 15091184, + 32550863, + 8818643, + 4244752, + 43123513, + 64565526, + 408838, + 13206998, + ), + u32x8::new( + 16405061, + 60379639, + 31489017, + 20949281, + 27568751, + 38734986, + 8364264, + 12451020, + ), + u32x8::new( + 16005217, + 58008076, + 1406778, + 26546927, + 39571784, + 56365493, + 31274296, + 8918790, + ), + u32x8::new( + 23271122, + 19453469, + 27718201, + 32742670, + 234332, + 36785342, + 22601675, + 14331046, + ), + u32x8::new( + 40636025, + 22442705, + 22115403, + 23745859, + 41164945, + 61012, + 12499614, + 542137, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 62776018, + 32835413, + 17373246, + 17187309, + 54469193, + 21770290, + 15923753, + 28996575, + ), + u32x8::new( + 59385210, + 63082298, + 12568449, + 8509004, + 9483342, + 16105238, + 5756054, + 26890758, + ), + u32x8::new( + 53987996, + 38201748, + 5521661, + 19060159, + 18663191, + 9093637, + 27786835, + 31189196, + ), + u32x8::new( + 65872678, + 43635130, + 27903055, + 25020300, + 65772737, + 38110437, + 5213502, + 21909342, + ), + u32x8::new( + 4438979, + 9680838, + 10212446, + 4764184, + 13235684, + 58245995, + 20264570, + 21024049, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 60835961, + 48209103, + 31049052, + 4688268, + 12426713, + 59829045, + 22302488, + 29008521, + ), + u32x8::new( + 50401667, + 29716596, + 23531224, + 7581281, + 49071895, + 6952617, + 14934683, + 8218256, + ), + u32x8::new( + 1601446, + 36631413, + 31774811, + 29625330, + 56786114, + 8331539, + 23129509, + 19783344, + ), + u32x8::new( + 59514327, + 64513110, + 1772300, + 5701338, + 5737511, + 16147555, + 9461515, + 5703271, + ), + u32x8::new( + 33072974, + 54300426, + 11940114, + 1308663, + 15627555, + 4931627, + 28443714, + 20924342, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 18135013, + 20358426, + 4922557, + 10015355, + 65729669, + 34786528, + 26248549, + 29194359, + ), + u32x8::new( + 797666, + 34997544, + 24316856, + 25107230, + 24612576, + 4761401, + 15307321, + 32404252, + ), + u32x8::new( + 16501152, + 60565831, + 9487105, + 9316022, + 24986054, + 31917592, + 3962024, + 2501883, + ), + u32x8::new( + 63356796, + 50432342, + 18044926, + 30566881, + 42032028, + 31415202, + 13524600, + 16119907, + ), + u32x8::new( + 3927286, + 57022374, + 9265437, + 21620772, + 19481940, + 3806938, + 24836192, + 14572399, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 10785787, + 46564798, + 368445, + 33181384, + 5319843, + 52687136, + 30347110, + 29837357, + ), + u32x8::new( + 56436732, + 47859251, + 24141084, + 22250712, + 59046084, + 4963427, + 33463413, + 17168859, + ), + u32x8::new( + 15512044, + 6366740, + 4737504, + 27644548, + 30307977, + 25037929, + 14593903, + 12836490, + ), + u32x8::new( + 63878897, + 34013023, + 5860752, + 7244096, + 3689461, + 57012135, + 18389096, + 11589351, + ), + u32x8::new( + 4682110, + 36302830, + 653422, + 22316819, + 14081831, + 5657024, + 11088376, + 24110612, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 39907267, + 45940262, + 24887471, + 18342609, + 878445, + 40456159, + 12019082, + 345107, + ), + u32x8::new( + 12794982, + 28893944, + 9447505, + 11387200, + 16961963, + 13916996, + 10893728, + 25898006, + ), + u32x8::new( + 44934162, + 53465865, + 3583620, + 1102334, + 53917811, + 63478576, + 2426066, + 10389549, + ), + u32x8::new( + 45096036, + 37595344, + 19367718, + 20257175, + 10280866, + 41653449, + 27665642, + 375926, + ), + u32x8::new( + 45847901, + 24064074, + 32494820, + 32204556, + 10720704, + 51079060, + 1297436, + 29853825, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 66303987, + 36060363, + 16494578, + 24962147, + 11971403, + 49538586, + 25060560, + 1964341, + ), + u32x8::new( + 25988481, + 27641502, + 24909517, + 27237087, + 66646363, + 52777626, + 16360849, + 10459972, + ), + u32x8::new( + 43930529, + 34374176, + 31225968, + 8807030, + 10394758, + 35904854, + 25325589, + 19335583, + ), + u32x8::new( + 25094697, + 34380951, + 20051185, + 32287161, + 11739332, + 53887441, + 30517319, + 26601892, + ), + u32x8::new( + 8868546, + 35635502, + 32513071, + 28248087, + 51946989, + 14222744, + 19198839, + 23261841, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 51218008, + 5070126, + 11046681, + 5320810, + 61212079, + 34104447, + 23895089, + 6460727, + ), + u32x8::new( + 39843528, + 46278671, + 10426120, + 25624792, + 66658766, + 37140083, + 28933107, + 12969597, + ), + u32x8::new( + 59635793, + 40220191, + 5751421, + 173680, + 58321825, + 740337, + 1412847, + 7682623, + ), + u32x8::new( + 975962, + 56440763, + 20812276, + 22631115, + 49095824, + 19883130, + 2419746, + 31043648, + ), + u32x8::new( + 66208703, + 39669328, + 22525915, + 3748897, + 65994776, + 34533552, + 8126286, + 18326047, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 64176557, + 3912400, + 19351673, + 30068471, + 31190055, + 24221683, + 33142424, + 28698542, + ), + u32x8::new( + 34784792, + 4109933, + 3867193, + 19557314, + 2112512, + 32715890, + 24550117, + 16595976, + ), + u32x8::new( + 35542761, + 48024875, + 10925431, + 31526577, + 66577735, + 23189821, + 13375709, + 1735095, + ), + u32x8::new( + 59699254, + 43854093, + 29783239, + 24777271, + 19600372, + 39924461, + 2896720, + 1472185, + ), + u32x8::new( + 56389656, + 35980854, + 33172342, + 1370336, + 23707480, + 57654949, + 7850973, + 12655016, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 38372660, + 57101970, + 7044964, + 12732710, + 57535705, + 6043201, + 30858914, + 10946592, + ), + u32x8::new( + 21023468, + 6946992, + 26403324, + 23901823, + 35695559, + 23440687, + 4763891, + 6514074, + ), + u32x8::new( + 28662273, + 30933699, + 9352242, + 26354829, + 37402243, + 3145176, + 8770289, + 525937, + ), + u32x8::new( + 54933102, + 36695832, + 3281859, + 4755022, + 23043294, + 32794379, + 15618886, + 23602412, + ), + u32x8::new( + 9931565, + 29897140, + 2480737, + 24193701, + 7833615, + 2284939, + 893926, + 13421882, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 22917795, + 22088359, + 28978099, + 19794863, + 60542318, + 29878494, + 31053731, + 9080720, + ), + u32x8::new( + 23679072, + 52547035, + 28424916, + 20647332, + 4008761, + 28267029, + 12961289, + 1589095, + ), + u32x8::new( + 55616194, + 26678929, + 14998265, + 23274397, + 54625466, + 46244264, + 28627706, + 33030665, + ), + u32x8::new( + 11527330, + 6449415, + 26531607, + 3472938, + 41541592, + 62607682, + 19862690, + 20564723, + ), + u32x8::new( + 32843805, + 49066843, + 28425824, + 19521495, + 48792073, + 48242878, + 27392443, + 13175986, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 16185025, + 61537525, + 2961305, + 1492442, + 25123147, + 3095034, + 31896958, + 33089615, + ), + u32x8::new( + 64748157, + 18336595, + 16522231, + 25426312, + 65718949, + 35485695, + 30554083, + 10205918, + ), + u32x8::new( + 39626934, + 39271045, + 16420458, + 9826240, + 56483981, + 27128085, + 3783403, + 13360006, + ), + u32x8::new( + 30793778, + 66771960, + 17241420, + 6564573, + 61102581, + 29974476, + 32385512, + 9011754, + ), + u32x8::new( + 28068166, + 11862220, + 14323567, + 12380617, + 52090465, + 16029056, + 24495309, + 21409233, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 59411973, + 57437124, + 11695483, + 17586857, + 16108987, + 43449109, + 31098002, + 6248476, + ), + u32x8::new( + 42258047, + 61595931, + 29308533, + 11742653, + 43042345, + 27373650, + 30165249, + 21929989, + ), + u32x8::new( + 49907221, + 9620337, + 21888081, + 20981082, + 56288861, + 61562203, + 33223566, + 3582446, + ), + u32x8::new( + 57535017, + 41003416, + 22080416, + 14463796, + 65518565, + 18127889, + 24370863, + 33332664, + ), + u32x8::new( + 66655380, + 6430175, + 471782, + 11947673, + 30596400, + 18898659, + 15930721, + 4211851, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 6757410, + 65455566, + 13584784, + 11362173, + 10797127, + 24451471, + 19541370, + 29309435, + ), + u32x8::new( + 40360156, + 17685025, + 18326181, + 3846903, + 13693365, + 63049479, + 31900359, + 23385063, + ), + u32x8::new( + 52455038, + 57513503, + 22163311, + 27095042, + 48610726, + 66454160, + 12085341, + 26357004, + ), + u32x8::new( + 22097042, + 14063840, + 6705778, + 14342902, + 66139825, + 20702105, + 31279090, + 7495745, + ), + u32x8::new( + 27360710, + 49314837, + 18774847, + 7146436, + 37066216, + 42004961, + 22409916, + 10524446, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 1497507, + 33054449, + 11839906, + 2960428, + 40538463, + 18884538, + 25018820, + 4073970, + ), + u32x8::new( + 54484385, + 43640735, + 2808257, + 20710708, + 39840730, + 27222424, + 21783544, + 11848522, + ), + u32x8::new( + 45765237, + 48200555, + 9299019, + 9393151, + 34818188, + 56098995, + 13575233, + 21012731, + ), + u32x8::new( + 4265428, + 49627650, + 24960282, + 9425650, + 47883651, + 2797524, + 11853190, + 22877329, + ), + u32x8::new( + 25008173, + 64199503, + 380047, + 12107343, + 12329448, + 11914399, + 764281, + 29687002, + ), + ])), + CachedPoint(FieldElement32x4([ + u32x8::new( + 35889734, + 23047226, + 4022841, + 7017445, + 7274086, + 53316179, + 25100176, + 15310676, + ), + u32x8::new( + 42409427, + 30270106, + 6823853, + 31551384, + 40645017, + 66489807, + 18021817, + 32669351, + ), + u32x8::new( + 39827134, + 43680850, + 28297996, + 20258133, + 26058742, + 52643238, + 22238331, + 21690533, + ), + u32x8::new( + 60808002, + 17499995, + 30042246, + 29310584, + 48219954, + 29389518, + 8680514, + 17844709, + ), + u32x8::new( + 6452896, + 50116553, + 9532047, + 26821214, + 44524351, + 50428429, + 21904953, + 12608048, + ), + ])), +]); diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 3c52bdb..55f1567 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -22,7 +22,7 @@ use subtle::ConditionallyAssignable; use subtle::Choice; use edwards; -use scalar_mul::window::{LookupTable, NafLookupTable5}; +use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use traits::Identity; @@ -307,6 +307,19 @@ impl<'a> From<&'a edwards::EdwardsPoint> for NafLookupTable5 { } } +impl<'a> From<&'a edwards::EdwardsPoint> for NafLookupTable8 { + fn from(point: &'a edwards::EdwardsPoint) -> Self { + let A = ExtendedPoint::from(*point); + let mut Ai = [CachedPoint::from(A); 64]; + let A2 = A.double(); + for i in 0..63 { + Ai[i + 1] = (&A2 + &Ai[i]).into(); + } + // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A, ..., 127A] + NafLookupTable8(Ai) + } +} + #[cfg(test)] mod test { use super::*; diff --git a/src/backend/avx2/scalar_mul/vartime_double_base.rs b/src/backend/avx2/scalar_mul/vartime_double_base.rs index 17d21f2..b3d6333 100644 --- a/src/backend/avx2/scalar_mul/vartime_double_base.rs +++ b/src/backend/avx2/scalar_mul/vartime_double_base.rs @@ -19,7 +19,7 @@ use backend::avx2::constants::BASEPOINT_ODD_LOOKUP_TABLE; /// Compute \\(aA + bB\\) in variable time, where \\(B\\) is the Ed25519 basepoint. pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint { let a_naf = a.non_adjacent_form(5); - let b_naf = b.non_adjacent_form(5); + let b_naf = b.non_adjacent_form(8); // Find starting index let mut i: usize = 255; From 68bbd1bd03cd80520ce2edc7b1a433655dbaa3fe Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 6 Apr 2018 12:11:59 -0700 Subject: [PATCH 5/9] Document algorithm for NAFs --- src/scalar.rs | 75 ++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 69 insertions(+), 6 deletions(-) diff --git a/src/scalar.rs b/src/scalar.rs index f4a86dd..c6b9917 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -80,6 +80,7 @@ pub struct Scalar { /// /// The integer representing this scalar must be bounded above by \\(2\^{255}\\), or equivalently the high bit of `bytes[31]` must be zero. /// + /// This ensures that there is room for a carry bit when computing a NAF representation. // XXX This is pub(crate) so we can write literal constants. If const fns were stable, we could make the Scalar constructors const fns and use those instead. pub(crate) bytes: [u8; 32], } @@ -497,17 +498,79 @@ impl Scalar { /// /// A width-\\(w\\) NAF of a positive integer \\(k\\) is an expression /// $$ - /// k = \sum_{i=0}\^n k\_i 2\^i, + /// k = \sum_{i=0}\^m n\_i 2\^i, /// $$ /// where each nonzero - /// coefficient \\(k\_i\\) is odd and bounded by \\(|k\_i| < 2\^{w-1}\\), - /// \\(k\_{n-1}\\) is nonzero, and at most one of any \\(w\\) consecutive + /// coefficient \\(n\_i\\) is odd and bounded by \\(|n\_i| < 2\^{w-1}\\), + /// \\(n\_{m-1}\\) is nonzero, and at most one of any \\(w\\) consecutive /// coefficients is nonzero. (Hankerson, Menezes, Vanstone; def 3.32). /// + /// The length of the NAF is at most one more than the length of + /// the binary representation of \\(k\\). This is why the + /// `Scalar` type maintains an invariant that the top bit is + /// \\(0\\), so that the NAF of a scalar has at most 256 digits. + /// /// Intuitively, this is like a binary expansion, except that we - /// allow some coefficients to grow up to \\(2\^{w-1}\\) so that the - /// nonzero coefficients are as sparse as possible. + /// allow some coefficients to grow in magnitude up to + /// \\(2\^{w-1}\\) so that the nonzero coefficients are as sparse + /// as possible. + /// + /// When doing scalar multiplication, we can then use a lookup + /// table of precomputed multiples of a point to add the nonzero + /// terms \\( k_i P \\). Using signed digits cuts the table size + /// in half, and using odd digits cuts the table size in half + /// again. + /// + /// To compute a \\(w\\)-NAF, we use a modification of Algorithm 3.35 of HMV: + /// + /// 1. \\( i \gets 0 \\) + /// 2. While \\( k \ge 1 \\): + /// 1. If \\(k\\) is odd, \\( n_i \gets k \operatorname{mods} 2^w \\), \\( k \gets k - n_i \\). + /// 2. If \\(k\\) is even, \\( n_i \gets 0 \\). + /// 3. \\( k \gets k / 2 \\), \\( i \gets i + 1 \\). + /// 3. Return \\( n_0, n_1, ... , \\) + /// + /// Here \\( \bar x = x \operatorname{mods} 2^w \\) means the + /// \\( \bar x \\) with \\( \bar x \equiv x \pmod{2^w} \\) and + /// \\( -2^{w-1} \leq \bar x < 2^w \\). + /// + /// We implement this by scanning across the bits of \\(k\\) from LSB to MSB. + /// Write the bits of \\(k\\) as + /// $$ + /// k = \sum\_{i=0}\^m k\_i 2^i, + /// $$ + /// and split the sum as + /// $$ + /// k = \sum\_{i=0}^{w-1} k\_i 2^i + 2^w \sum\_{i=0} k\_{i+w} 2^i + /// $$ + /// where the first part is \\( k \mod 2^w \\). + /// + /// If \\( k \mod 2^w\\) is odd, and \\( k \mod 2^w < 2^{w-1} \\), then we emit + /// \\( n_0 = k \mod 2^w \\). Instead of computing + /// \\( k - n_0 \\), we just advance \\(w\\) bits and reindex. + /// + /// If \\( k \mod 2^w\\) is odd, and \\( k \mod 2^w \ge 2^{w-1} \\), then + /// \\( n_0 = k \operatorname{mods} 2^w = k \mod 2^w - 2^w \\). + /// The quantity \\( k - n_0 \\) is + /// $$ + /// \begin{aligned} + /// k - n_0 &= \sum\_{i=0}^{w-1} k\_i 2^i + 2^w \sum\_{i=0} k\_{i+w} 2^i + /// - \sum\_{i=0}^{w-1} k\_i 2^i + 2^w \\\\ + /// &= 2^w + 2^w \sum\_{i=0} k\_{i+w} 2^i + /// \end{aligned} + /// $$ + /// so instead of computing the subtraction, we can set a carry + /// bit, advance \\(w\\) bits, and reindex. + /// + /// If \\( k \mod 2^w\\) is even, we emit \\(0\\), advance 1 bit + /// and reindex. In fact, by setting all digits to \\(0\\) + /// initially, we don't need to emit anything. pub(crate) fn non_adjacent_form(&self, w: usize) -> [i8; 256] { + // required by the NAF definition + debug_assert!( w >= 2 ); + // required so that the NAF digits fit in i8 + debug_assert!( w <= 8 ); + use byteorder::{ByteOrder, LittleEndian}; let mut naf = [0i8; 256]; @@ -564,7 +627,7 @@ impl Scalar { /// $$ /// a = a\_0 + a\_1 16\^1 + \cdots + a_{63} 16\^{63}, /// $$ - /// with \\(-8 \leq a_i < 8\\) for \\(0 \leq i < 63\\) and \\(-8 \leq a_63 \leq 8\\). + /// with \\(-8 \leq a_i < 8\\) for \\(0 \leq i < 63\\) and \\(-8 \leq a_{63} \leq 8\\). pub(crate) fn to_radix_16(&self) -> [i8; 64] { debug_assert!(self[31] <= 127); let mut output = [0i8; 64]; From 9fc5602ce79103525faa5c4d0c95b758b8c70540 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Sun, 8 Apr 2018 15:15:23 -0700 Subject: [PATCH 6/9] Split AVX2 docs into markdown file --- docs/avx2-notes.md | 464 +++++++++++++++++++++++++++++++++++++++ src/backend/avx2/mod.rs | 467 +--------------------------------------- src/ristretto.rs | 2 + 3 files changed, 468 insertions(+), 465 deletions(-) create mode 100644 docs/avx2-notes.md diff --git a/docs/avx2-notes.md b/docs/avx2-notes.md new file mode 100644 index 0000000..2e70d82 --- /dev/null +++ b/docs/avx2-notes.md @@ -0,0 +1,464 @@ +An implementation of group operations on the twisted Edwards form of +Curve25519, using AVX2 to implement the 4-way parallel formulas of +Hisil, Wong, Carter, and Dawson (HWCD). + +Their 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08], which +introduced the extended coordinates used in other parts of `-dalek`, +also describes 4-way parallel formulas for point addition and +doubling: + +* a unified addition algorithm taking an effective \\(2\mathbf M + +1\mathbf D\\); + +* a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf +S\\); + +* a dedicated (i.e., for distinct points) addition algorithm taking +an effective \\(2 \mathbf M \\). + +Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of +multiplication and squaring of generic field elements and \\(\mathbf +D\\) represents the cost of multiplication by a curve constant. + +Currently, this implementation uses only the first two algorithms. + +# Parallel formulas + +The doubling formula is presented in the HWCD paper as follows: + +| Cost | Processor 1 | Processor 2 | Processor 3 | Processor 4 | +|------------------|--------------------------------|--------------------------------|--------------------------------|--------------------------------| +| | idle | idle | idle | \\( R\_1 \gets X\_1 + Y\_1 \\) | +| \\(1\mathbf S\\) | \\( R\_2 \gets X\_1\^2 \\) | \\( R\_3 \gets Y\_1\^2 \\) | \\( R\_4 \gets Z\_1\^2 \\) | \\( R\_5 \gets R\_1\^2 \\) | +| | \\( R\_6 \gets R\_2 + R\_3 \\) | \\( R\_7 \gets R\_2 - R\_3 \\) | \\( R\_4 \gets 2 R\_4 \\) | idle | +| | idle | \\( R\_1 \gets R\_4 + R\_7 \\) | idle | \\( R\_2 \gets R\_6 - R\_5 \\) | +| \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_6 R\_7 \\) | \\( T\_3 \gets R\_2 R\_6 \\) | \\( Z\_3 \gets R\_1 R\_7 \\) | + +and the unified addition algorithm is presented as follows: + +| Cost | Processor 1 | Processor 2 | Processor 3 | Processor 4 | +|------------------|--------------------------------|--------------------------------|--------------------------------|--------------------------------| +| | \\( R\_1 \gets Y\_1 - X\_1 \\) | \\( R\_2 \gets Y\_2 - X\_2 \\) | \\( R\_3 \gets Y\_1 + X\_1 \\) | \\( R\_4 \gets Y\_2 + X\_2 \\) | +| \\(1\mathbf M\\) | \\( R\_5 \gets R\_1 R\_2 \\) | \\( R\_6 \gets R\_3 R\_4 \\) | \\( R\_7 \gets T\_1 T\_2 \\) | \\( R\_8 \gets Z\_1 Z\_2 \\) | +| \\(1\mathbf D\\) | idle | idle | \\( R\_7 \gets k R\_7 \\) | \\( R\_8 \gets 2 R\_8 \\) | +| | \\( R\_1 \gets R\_6 - R\_5 \\) | \\( R\_2 \gets R\_8 - R\_7 \\) | \\( R\_3 \gets R\_8 + R\_7 \\) | \\( R\_4 \gets R\_6 + R\_5 \\) | +| \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_3 R\_4 \\) | \\( T\_3 \gets R\_1 R\_4 \\) | \\( Z\_3 \gets R\_2 R\_3 \\) | + +Here \\( k = 2d \\) is a curve constant. + +# Implementation strategy + +For a software implementation, each "processor"'s operations are too +low-latency to parallelize across threads. However, the main cost +is in the multiplication and squaring steps, which share a single +instruction. + +Our strategy is to implement 4-wide multiplication and squaring +using one 64-bit AVX2 lane for each field element. Field elements +are represented in the usual way as 10 `u32` limbs in radix +\\(25.5\\) (i.e., alternating between \\(2\^{26}\\) for even limbs +and \\(2\^{25}\\) for odd limbs). This has the effect that passing +between the parallel 32-bit AVX2 representation and the serial +64-bit representation amounts to regrouping digits. + +The addition and subtraction steps are done largely serially, using +masking to handle the instruction divergence. The remaining +obstacle to parallelism is the multiplication by the curve constant +\\(k = 2d\\). In the Curve25519 case, this is + +$$ k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. $$ + +HWCD suggest parallelising this step by breaking \\(k\\) into four +parts as \\(k = k_0 + 2\^n k_1 + 2\^{2n} k_2 + 2\^{3n} k_3 \\) and +computing \\(k_i R_7 \\) in parallel. However, this would be +somewhat awkward in our case, since we would normally represent +\\(k\\) as \\( 10 \\) 32-bit limbs, and \\(10 \\) is not divisible +by \\(4\\), so we would need a specialized routine to perform a +vectorized multiplication by 64-bit constants. + +Instead, since we are working projectively, we can multiply +\\(R_7\\) by \\( -2\cdot 121665 \\) and multiply the other three +variables by \\(121666\\). This trick was suggested by Mike +Hamburg. Ignoring the sign for the moment, since +\\(2 \cdot 121666 < 2\^{18}\\), all these constants fit in 32 bits, +so (up to sign) this can be done in parallel as four multiplications +by small constants \\( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \\). + +How do we handle the sign? +Since we're primarily interested in Ristretto performance, not +Curve25519 performance, we could alternately work on the +\\(4\\)-isogenous "IsoEd25519" curve, which has \\(d = 121665\\). +However, this would only save the negation step, since multiplying +one field element by a 32-bit constant is not much easier than +multiplying four field elements by 32-bit constants, and it would +prevent accelerating Curve25519, so we don't make this choice. +Instead, we just negate one lane, and move the \\(1 \mathbf D\\) +into precomputation (see below). + +The 4-wide formulas of the HWCD paper do not seem to have been +implemented using SIMD before. The HWCD paper also describes and +analyzes a 2-wide variant of the Montgomery ladder (for comparison +with parallel Edwards formulas); this strategy was used in 2015 by +Tung Chou's `sandy2x` implementation, which used a 2-wide field +implementation in 128-bit vector registers. + +Curiously, however, although the [`sandy2x` paper][sandy2x] also +implements Edwards arithmetic, and cites the HWCD paper, it doesn't +mention or discuss the parallel formulas from HWCD, or that the +2-wide Montgomery formulas it uses were previously published there. +There is also a 2015 paper by Hernández and López on using AVX2 for +the X25519 Montgomery ladder, but neither the paper nor the code are +publicly available, and it apparently gives only a [slight +speedup][avx2trac], suggesting that it also overlooked the +HWCD formulas. + +HWCD also suggest using a mixed representation, passing between \\( +\mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\) +"projective" coordinates, where doubling is slightly cheaper (saving +about \\(\mathbf 1M\\). This approach is used for the +non-vectorized `u32` and `u64` backends, and more +details on the different coordinate systems can be found in the +`curve_models` module documentation. + +This optimization is not compatible with the parallel formulas, which are +therefore slightly less efficient when counting the total number of +field multiplications and squarings. In particular, vectorized doublings +are less efficient than serial doublings. +In addition, the parallel formulas can only use a \\( 32 \times 32 +\rightarrow 64 \\)-bit integer multiplier, so the speedup from +vectorization must overcome the disadvantage of losing the \\( 64 +\times 64 \rightarrow 128\\)-bit (serial) integer multiplier. + +# Tweaked formulas + +After tweaking the formulas as described above, we obtain the +following. To avoid confusion with the original HWCD formulas, +temporary variables are named \\(S\\) instead of \\(R\\) and are in +static single-assignment (SSA) form. + +## Addition + +To add points \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and \\(P_2 = (X_2 +: Y_2 : Z_2 : T_2 ) \\), we compute + +$$ +\begin{aligned} +S\_0 &\gets Y\_1 - X\_1 \\\\ +S\_1 &\gets Y\_1 + X\_1 \\\\ +S\_2 &\gets Y\_2 - X\_2 \\\\ +S\_3 &\gets Y\_2 + X\_2 +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_4 &\gets S\_0 S\_2 \\\\ +S\_5 &\gets S\_1 S\_3 \\\\ +S\_6 &\gets Z\_1 Z\_2 \\\\ +S\_7 &\gets T\_1 T\_2 +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_8 &\gets S\_4 \cdot 121666 \\\\ +S\_9 &\gets S\_5 \cdot 121666 \\\\ +S\_{10} &\gets S\_6 \cdot 2 \cdot 121666 \\\\ +S\_{11} &\gets S\_7 \cdot -2 \cdot 121665 +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_{12} &\gets S\_9 - S\_8 \\\\ +S\_{13} &\gets S\_9 + S\_8 \\\\ +S\_{14} &\gets S\_{10} - S\_{11} \\\\ +S\_{15} &\gets S\_{10} + S\_{11} +\end{aligned} +$$ + +$$ +\begin{aligned} +X\_3 &\gets S\_{12} S\_{14} \\\\ +Y\_3 &\gets S\_{15} S\_{13} \\\\ +Z\_3 &\gets S\_{15} S\_{14} \\\\ +T\_3 &\gets S\_{12} S\_{13} +\end{aligned} +$$ + +to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). + +## Readdition + +If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we can precompute + +$$ +\begin{aligned} +S\_2 &\gets Y\_2 - X\_2 \\\\ +S\_3 &\gets Y\_2 + X\_2 +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_2' &\gets S\_2 \cdot 121666 \\\\ +S\_3' &\gets S\_3 \cdot 121666 \\\\ +Z\_2' &\gets Z\_2 \cdot 2 \cdot 121666 \\\\ +T\_2' &\gets T\_2 \cdot -2 \cdot 121665 \\\\ +\end{aligned} +$$ + +to obtain the `CachedPoint` \\( (S\_2', S\_3', Z\_2', T\_2') \\). +This precomputation is essentially the same as that suggested in +§3.1 of HWCD, with the difference that the multiplication by the curve +constant \\( -121665 / 121666 \\) is spread over all four +coordinates, to allow a vectorized computation of four +multiplications of small constants instead of a serial computation +of multiplication by a large constant. + +To perform readdition of \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and +\\(P_2 = (S\_2', S\_3', Z\_2', T\_2') \\), we compute + +$$ +\begin{aligned} +S\_0 &\gets Y\_1 - X\_1 \\\\ +S\_1 &\gets Y\_1 + X\_1 +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_8 &\gets S\_0 S\_2' \\\\ +S\_9 &\gets S\_1 S\_3' \\\\ +S\_{10} &\gets Z\_1 Z\_2' \\\\ +S\_{11} &\gets T\_1 T\_2' +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_{12} &\gets S\_9 - S\_8 \\\\ +S\_{13} &\gets S\_9 + S\_8 \\\\ +S\_{14} &\gets S\_{10} - S\_{11} \\\\ +S\_{15} &\gets S\_{10} + S\_{11} +\end{aligned} +$$ + +$$ +\begin{aligned} +X\_3 &\gets S\_{12} S\_{14} \\\\ +Y\_3 &\gets S\_{15} S\_{13} \\\\ +Z\_3 &\gets S\_{15} S\_{14} \\\\ +T\_3 &\gets S\_{12} S\_{13} +\end{aligned} +$$ + +to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). + +Compared to the addition formulas above, this saves \\( 1\mathbf D \\). + +## Doubling + +To double a point \\( P = (X\_1 : Y\_1 : Z\_1 : T\_1) \\), we compute + +$$ S\_0 \gets X\_1 + Y\_1 $$ + +$$ +\begin{aligned} +S\_1 &\gets X\_1\^2 \\\\ +S\_2 &\gets Y\_1\^2 \\\\ +S\_3 &\gets Z\_1\^2 \\\\ +S\_4 &\gets S\_0\^2 +\end{aligned} +$$ + +$$ +\begin{aligned} +S\_5 &\gets S\_1 + S\_2 \\\\ +S\_6 &\gets S\_1 - S\_2 \\\\ +S\_7 &\gets 2S\_3 \\\\ +S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\ +S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4 +\end{aligned} +$$ + +$$ +\begin{aligned} +X\_3 &\gets S\_8 S\_9 \\\\ +Y\_3 &\gets S\_5 S\_6 \\\\ +Z\_3 &\gets S\_8 S\_6 \\\\ +T\_3 &\gets S\_5 S\_9 +\end{aligned} +$$ + +to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = [2]P\_1 \\). + +Performing too many intermediate additions and subtractions grows +the bounds beyond what is allowed as input to multiplication, +forcing an extra carry pass. However, it is just possible to avoid +this by rearranging signs. + +Assume that the bounds on the limbs of each field element are +parameterized by \\( b \in \mathbb R \\) representing the excess +bits, so that each limb is bounded by either \\( 2\^{25} \\) or \\( +2\^{26} \\). + +The multiplication routine requires that its inputs are bounded by +\\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\) +into 32 bits. Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\) +when \\( x < 2\^{27.75} = 2\^{26 + 1.75} \\). However, this is only +required for one of the inputs; the other can grow up to \\( b < 2.5 +\\). + +Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as + +$$ +\begin{matrix} + & S\_1 & S\_1 & S\_1 & S\_1 \\\\ ++& S\_2 & & & S\_2 \\\\ ++& & & S\_3 & \\\\ ++& & & S\_3 & \\\\ ++& & 2p & 2p & 2p \\\\ +-& & S\_2 & S\_2 & \\\\ +-& & & & S\_4 \\\\ +=& S\_5 & S\_6 & S\_8 & S\_9 +\end{matrix} +$$ + +results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for +\\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute +are then + +$$ +\begin{aligned} +X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.00) \\\\ +Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ +Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ +T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 2.00) +\end{aligned} +$$ + +which are too large. However, if we flip the sign of \\( S\_4 = +S\_0\^2 \\) during squaring, so that we output \\(S\_4' = -S\_4 +\pmod p\\), then we can compute + +$$ +\begin{matrix} + & S\_1 & S\_1 & S\_1 & S\_1 \\\\ ++& S\_2 & & & S\_2 \\\\ ++& & & S\_3 & \\\\ ++& & & S\_3 & \\\\ ++& & & & S\_4' \\\\ ++& & 2p & 2p & \\\\ +-& & S\_2 & S\_2 & \\\\ +=& S\_5 & S\_6 & S\_8 & S\_9 +\end{matrix} +$$ + +resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for +\\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute +are then + +$$ +\begin{aligned} +X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.59) \\\\ +Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ +Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ +T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 1.59) +\end{aligned} +$$ + +whose right-hand sides are all bounded with \\( b < 1.75 \\) and +whose left-hand sides are all bounded with \\( b < 2.5 \\). + +# Field element representation + +The field element representation is oriented around the AVX2 +`vpmuluqdq` instruction, which multiplies the low 32 bits of each +64-bit lane of each operand to produce a 64-bit result. + +```text,no_run +(a1 ?? b1 ?? c1 ?? d1 ??) +(a2 ?? b2 ?? c2 ?? d2 ??) + +(a1*a2 b1*b2 c1*c2 d1*d2) +``` + +To unpack 32-bit values into 64-bit lanes for use in multiplication +it would be convenient to use the `vpunpck[lh]dq` instructions, +which unpack and interleave the low and high 32-bit lanes of two +source vectors. +However, the AVX2 versions of these instructions are designed to +operate only within 128-bit lanes of the 256-bit vectors, so that +interleaving the low lanes of `(a0 b0 c0 d0 a1 b1 c1 d1)` with zero +gives `(a0 00 b0 00 a1 00 b1 00)`. Instead, we pre-shuffle the data +layout as `(a0 b0 a1 b1 c0 d0 c1 d1)` so that we can unpack the +"low" and "high" parts as + +```text,no_run +(a0 00 b0 00 c0 00 d0 00) +(a1 00 b1 00 c1 00 d1 00) +``` + +The data layout for a vector of four field elements \\( (a,b,c,d) +\\) with limbs \\( a_0, a_1, \ldots, a_9 \\) is as `[u32x8; 5]` in +the form + +```text,no_run +(a0 b0 a1 b1 c0 d0 c1 d1) +(a2 b2 a3 b3 c2 d2 c3 d3) +(a4 b4 a5 b5 c4 d4 c5 d5) +(a6 b6 a7 b7 c6 d6 c7 d7) +(a8 b8 a9 b9 c8 d8 c9 d9) +``` + +Since this breaks cleanly into two 128-bit lanes, it may be possible +to adapt it to 128-bit vector instructions such as NEON without too +much difficulty. + +Going the other direction, to extend this to AVX512, we could either +run two point operations in parallel in lower and upper halves of +the registers, or use 2-way parallelism within a field operation. + +We don't attempt to use AVX2 for serial field element computations +such as inversion, since wherever we have AVX2 we also have `mulx`. +However, it might be useful for batched inverse square-root +computations, which can't be batched in the same way inversions can. + +# Implementation details + +The implementation uses the unstable `stdsimd` crate to provide AVX2 +intrinsics, and the code is not yet cleanly factored between the +field element parts and the point parts. + +When compiling with AVX512VL, LLVM is able to use the extra +`ymm16..ymm31` registers to reduce register pressure, and avoid +spills during field multiplication. This gives a small but +noticeable speedup. + +The addition and subtraction steps involve masking, to apply +operations to a single lane of the vector. AVX512VL extends the +predication features of AVX512 to AVX2 code and would probably be +beneficial. Unfortunately, LLVM is currently unable to lower `op + +blend` into an AVX512VL masked operation. However, the explicitly +masked versions of the intrinsics seem to produce the same LLVM IR +as an `op + blend`, so hopefully this will improve as the AVX512 +support in LLVM improves. + +When used for constant-time variable-base scalar multiplication, +this strategy (using AVX2) gives a significant speedup over the +serial implementation (using the \\(64 \times 64\\) multiplier) of +approximately 1.6x for Skylake-X with `target_cpu=skylake` (using AVX2), of +approximately 1.8x for Skylake-X with `target_cpu=skylake-avx512` (using the extra +`ymm16..ymm31` registers from AVX512VL), and of approximately 1.0x +for Ryzen (which implements AVX2 at half rate). + +When used for variable-time double-base scalar multiplication +\\( aA + bB \\) for fixed \\(B\\) (as in, e.g., signature verification), +this strategy provides a 1.4x speedup on Skylake-X over the same +operation as implemented in `ed25519-donna`, the fastest +production-quality Ed25519 implementation. + +[sandy2x]: https://eprint.iacr.org/2015/943.pdf +[avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28 +[hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf diff --git a/src/backend/avx2/mod.rs b/src/backend/avx2/mod.rs index bd3e860..94cf55a 100644 --- a/src/backend/avx2/mod.rs +++ b/src/backend/avx2/mod.rs @@ -8,471 +8,8 @@ // - Isis Agora Lovecruft // - Henry de Valence -//! An implementation of group operations on the twisted Edwards form of -//! Curve25519, using AVX2 to implement the 4-way parallel formulas of -//! Hisil, Wong, Carter, and Dawson (HWCD). -//! -//! Their 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08], which -//! introduced the extended coordinates used in other parts of `-dalek`, -//! also describes 4-way parallel formulas for point addition and -//! doubling: -//! -//! * a unified addition algorithm taking an effective \\(2\mathbf M + -//! 1\mathbf D\\); -//! -//! * a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf -//! S\\); -//! -//! * a dedicated (i.e., for distinct points) addition algorithm taking -//! an effective \\(2 \mathbf M \\). -//! -//! Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of -//! multiplication and squaring of generic field elements and \\(\mathbf -//! D\\) represents the cost of multiplication by a curve constant. -//! -//! Currently, this implementation uses only the first two algorithms. -//! -//! # Parallel formulas -//! -//! The doubling formula is presented in the HWCD paper as follows: -//! -//! | Cost | Processor 1 | Processor 2 | Processor 3 | Processor 4 | -//! |------------------|--------------------------------|--------------------------------|--------------------------------|--------------------------------| -//! | | idle | idle | idle | \\( R\_1 \gets X\_1 + Y\_1 \\) | -//! | \\(1\mathbf S\\) | \\( R\_2 \gets X\_1\^2 \\) | \\( R\_3 \gets Y\_1\^2 \\) | \\( R\_4 \gets Z\_1\^2 \\) | \\( R\_5 \gets R\_1\^2 \\) | -//! | | \\( R\_6 \gets R\_2 + R\_3 \\) | \\( R\_7 \gets R\_2 - R\_3 \\) | \\( R\_4 \gets 2 R\_4 \\) | idle | -//! | | idle | \\( R\_1 \gets R\_4 + R\_7 \\) | idle | \\( R\_2 \gets R\_6 - R\_5 \\) | -//! | \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_6 R\_7 \\) | \\( T\_3 \gets R\_2 R\_6 \\) | \\( Z\_3 \gets R\_1 R\_7 \\) | -//! -//! and the unified addition algorithm is presented as follows: -//! -//! | Cost | Processor 1 | Processor 2 | Processor 3 | Processor 4 | -//! |------------------|--------------------------------|--------------------------------|--------------------------------|--------------------------------| -//! | | \\( R\_1 \gets Y\_1 - X\_1 \\) | \\( R\_2 \gets Y\_2 - X\_2 \\) | \\( R\_3 \gets Y\_1 + X\_1 \\) | \\( R\_4 \gets Y\_2 + X\_2 \\) | -//! | \\(1\mathbf M\\) | \\( R\_5 \gets R\_1 R\_2 \\) | \\( R\_6 \gets R\_3 R\_4 \\) | \\( R\_7 \gets T\_1 T\_2 \\) | \\( R\_8 \gets Z\_1 Z\_2 \\) | -//! | \\(1\mathbf D\\) | idle | idle | \\( R\_7 \gets k R\_7 \\) | \\( R\_8 \gets 2 R\_8 \\) | -//! | | \\( R\_1 \gets R\_6 - R\_5 \\) | \\( R\_2 \gets R\_8 - R\_7 \\) | \\( R\_3 \gets R\_8 + R\_7 \\) | \\( R\_4 \gets R\_6 + R\_5 \\) | -//! | \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_3 R\_4 \\) | \\( T\_3 \gets R\_1 R\_4 \\) | \\( Z\_3 \gets R\_2 R\_3 \\) | -//! -//! Here \\( k = 2d \\) is a curve constant. -//! -//! # Implementation strategy -//! -//! For a software implementation, each "processor"'s operations are too -//! low-latency to parallelize across threads. However, the main cost -//! is in the multiplication and squaring steps, which share a single -//! instruction. -//! -//! Our strategy is to implement 4-wide multiplication and squaring -//! using one 64-bit AVX2 lane for each field element. Field elements -//! are represented in the usual way as 10 `u32` limbs in radix -//! \\(25.5\\) (i.e., alternating between \\(2\^{26}\\) for even limbs -//! and \\(2\^{25}\\) for odd limbs). This has the effect that passing -//! between the parallel 32-bit AVX2 representation and the serial -//! 64-bit representation amounts to regrouping digits. -//! -//! The addition and subtraction steps are done largely serially, using -//! masking to handle the instruction divergence. The remaining -//! obstacle to parallelism is the multiplication by the curve constant -//! \\(k = 2d\\). In the Curve25519 case, this is -//! -//! $$ k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. $$ -//! -//! HWCD suggest parallelising this step by breaking \\(k\\) into four -//! parts as \\(k = k_0 + 2\^n k_1 + 2\^{2n} k_2 + 2\^{3n} k_3 \\) and -//! computing \\(k_i R_7 \\) in parallel. However, this would be -//! somewhat awkward in our case, since we would normally represent -//! \\(k\\) as \\( 10 \\) 32-bit limbs, and \\(10 \\) is not divisible -//! by \\(4\\), so we would need a specialized routine to perform a -//! vectorized multiplication by 64-bit constants. -//! -//! Instead, since we are working projectively, we can multiply -//! \\(R_7\\) by \\( -2\cdot 121665 \\) and multiply the other three -//! variables by \\(121666\\). This trick was suggested by Mike -//! Hamburg. Ignoring the sign for the moment, since -//! \\(2 \cdot 121666 < 2\^{18}\\), all these constants fit in 32 bits, -//! so (up to sign) this can be done in parallel as four multiplications -//! by small constants \\( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \\). -//! -//! How do we handle the sign? -//! Since we're primarily interested in Ristretto performance, not -//! Curve25519 performance, we could alternately work on the -//! \\(4\\)-isogenous "IsoEd25519" curve, which has \\(d = 121665\\). -//! However, this would only save the negation step, since multiplying -//! one field element by a 32-bit constant is not much easier than -//! multiplying four field elements by 32-bit constants, and it would -//! prevent accelerating Curve25519, so we don't make this choice. -//! Instead, we just negate one lane, and move the \\(1 \mathbf D\\) -//! into precomputation (see below). -//! -//! The 4-wide formulas of the HWCD paper do not seem to have been -//! implemented using SIMD before. The HWCD paper also describes and -//! analyzes a 2-wide variant of the Montgomery ladder (for comparison -//! with parallel Edwards formulas); this strategy was used in 2015 by -//! Tung Chou's `sandy2x` implementation, which used a 2-wide field -//! implementation in 128-bit vector registers. -//! -//! Curiously, however, although the [`sandy2x` paper][sandy2x] also -//! implements Edwards arithmetic, and cites the HWCD paper, it doesn't -//! mention or discuss the parallel formulas from HWCD, or that the -//! 2-wide Montgomery formulas it uses were previously published there. -//! There is also a 2015 paper by Hernández and López on using AVX2 for -//! the X25519 Montgomery ladder, but neither the paper nor the code are -//! publicly available, and it apparently gives only a [slight -//! speedup][avx2trac], suggesting that it also overlooked the -//! HWCD formulas. -//! -//! HWCD also suggest using a mixed representation, passing between \\( -//! \mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\) -//! "projective" coordinates, where doubling is slightly cheaper (saving -//! about \\(\mathbf 1M\\). This approach is used for the -//! non-vectorized `u32` and `u64` backends, and more -//! details on the different coordinate systems can be found in the -//! `curve_models` module documentation. -//! -//! This optimization is not compatible with the parallel formulas, which are -//! therefore slightly less efficient when counting the total number of -//! field multiplications and squarings. In particular, vectorized doublings -//! are less efficient than serial doublings. -//! In addition, the parallel formulas can only use a \\( 32 \times 32 -//! \rightarrow 64 \\)-bit integer multiplier, so the speedup from -//! vectorization must overcome the disadvantage of losing the \\( 64 -//! \times 64 \rightarrow 128\\)-bit (serial) integer multiplier. -//! -//! # Tweaked formulas -//! -//! After tweaking the formulas as described above, we obtain the -//! following. To avoid confusion with the original HWCD formulas, -//! temporary variables are named \\(S\\) instead of \\(R\\) and are in -//! static single-assignment (SSA) form. -//! -//! ## Addition -//! -//! To add points \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and \\(P_2 = (X_2 -//! : Y_2 : Z_2 : T_2 ) \\), we compute -//! -//! $$ -//! \begin{aligned} -//! S\_0 &\gets Y\_1 - X\_1 \\\\ -//! S\_1 &\gets Y\_1 + X\_1 \\\\ -//! S\_2 &\gets Y\_2 - X\_2 \\\\ -//! S\_3 &\gets Y\_2 + X\_2 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_4 &\gets S\_0 S\_2 \\\\ -//! S\_5 &\gets S\_1 S\_3 \\\\ -//! S\_6 &\gets Z\_1 Z\_2 \\\\ -//! S\_7 &\gets T\_1 T\_2 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_8 &\gets S\_4 \cdot 121666 \\\\ -//! S\_9 &\gets S\_5 \cdot 121666 \\\\ -//! S\_{10} &\gets S\_6 \cdot 2 \cdot 121666 \\\\ -//! S\_{11} &\gets S\_7 \cdot -2 \cdot 121665 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_{12} &\gets S\_9 - S\_8 \\\\ -//! S\_{13} &\gets S\_9 + S\_8 \\\\ -//! S\_{14} &\gets S\_{10} - S\_{11} \\\\ -//! S\_{15} &\gets S\_{10} + S\_{11} -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! X\_3 &\gets S\_{12} S\_{14} \\\\ -//! Y\_3 &\gets S\_{15} S\_{13} \\\\ -//! Z\_3 &\gets S\_{15} S\_{14} \\\\ -//! T\_3 &\gets S\_{12} S\_{13} -//! \end{aligned} -//! $$ -//! -//! to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). -//! -//! ## Readdition -//! -//! If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we can precompute -//! -//! $$ -//! \begin{aligned} -//! S\_2 &\gets Y\_2 - X\_2 \\\\ -//! S\_3 &\gets Y\_2 + X\_2 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_2' &\gets S\_2 \cdot 121666 \\\\ -//! S\_3' &\gets S\_3 \cdot 121666 \\\\ -//! Z\_2' &\gets Z\_2 \cdot 2 \cdot 121666 \\\\ -//! T\_2' &\gets T\_2 \cdot -2 \cdot 121665 \\\\ -//! \end{aligned} -//! $$ -//! -//! to obtain the `CachedPoint` \\( (S\_2', S\_3', Z\_2', T\_2') \\). -//! This precomputation is essentially the same as that suggested in -//! §3.1 of HWCD, with the difference that the multiplication by the curve -//! constant \\( -121665 / 121666 \\) is spread over all four -//! coordinates, to allow a vectorized computation of four -//! multiplications of small constants instead of a serial computation -//! of multiplication by a large constant. -//! -//! To perform readdition of \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and -//! \\(P_2 = (S\_2', S\_3', Z\_2', T\_2') \\), we compute -//! -//! $$ -//! \begin{aligned} -//! S\_0 &\gets Y\_1 - X\_1 \\\\ -//! S\_1 &\gets Y\_1 + X\_1 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_8 &\gets S\_0 S\_2' \\\\ -//! S\_9 &\gets S\_1 S\_3' \\\\ -//! S\_{10} &\gets Z\_1 Z\_2' \\\\ -//! S\_{11} &\gets T\_1 T\_2' -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_{12} &\gets S\_9 - S\_8 \\\\ -//! S\_{13} &\gets S\_9 + S\_8 \\\\ -//! S\_{14} &\gets S\_{10} - S\_{11} \\\\ -//! S\_{15} &\gets S\_{10} + S\_{11} -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! X\_3 &\gets S\_{12} S\_{14} \\\\ -//! Y\_3 &\gets S\_{15} S\_{13} \\\\ -//! Z\_3 &\gets S\_{15} S\_{14} \\\\ -//! T\_3 &\gets S\_{12} S\_{13} -//! \end{aligned} -//! $$ -//! -//! to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). -//! -//! Compared to the addition formulas above, this saves \\( 1\mathbf D \\). -//! -//! ## Doubling -//! -//! To double a point \\( P = (X\_1 : Y\_1 : Z\_1 : T\_1) \\), we compute -//! -//! $$ S\_0 \gets X\_1 + Y\_1 $$ -//! -//! $$ -//! \begin{aligned} -//! S\_1 &\gets X\_1\^2 \\\\ -//! S\_2 &\gets Y\_1\^2 \\\\ -//! S\_3 &\gets Z\_1\^2 \\\\ -//! S\_4 &\gets S\_0\^2 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! S\_5 &\gets S\_1 + S\_2 \\\\ -//! S\_6 &\gets S\_1 - S\_2 \\\\ -//! S\_7 &\gets 2S\_3 \\\\ -//! S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\ -//! S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4 -//! \end{aligned} -//! $$ -//! -//! $$ -//! \begin{aligned} -//! X\_3 &\gets S\_8 S\_9 \\\\ -//! Y\_3 &\gets S\_5 S\_6 \\\\ -//! Z\_3 &\gets S\_8 S\_6 \\\\ -//! T\_3 &\gets S\_5 S\_9 -//! \end{aligned} -//! $$ -//! -//! to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = [2]P\_1 \\). -//! -//! Performing too many intermediate additions and subtractions grows -//! the bounds beyond what is allowed as input to multiplication, -//! forcing an extra carry pass. However, it is just possible to avoid -//! this by rearranging signs. -//! -//! Assume that the bounds on the limbs of each field element are -//! parameterized by \\( b \in \mathbb R \\) representing the excess -//! bits, so that each limb is bounded by either \\( 2\^{25} \\) or \\( -//! 2\^{26} \\). -//! -//! The multiplication routine requires that its inputs are bounded by -//! \\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\) -//! into 32 bits. Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\) -//! when \\( x < 2\^{27.75} = 2\^{26 + 1.75} \\). However, this is only -//! required for one of the inputs; the other can grow up to \\( b < 2.5 -//! \\). -//! -//! Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as -//! -//! $$ -//! \begin{matrix} -//! & S\_1 & S\_1 & S\_1 & S\_1 \\\\ -//! +& S\_2 & & & S\_2 \\\\ -//! +& & & S\_3 & \\\\ -//! +& & & S\_3 & \\\\ -//! +& & 2p & 2p & 2p \\\\ -//! -& & S\_2 & S\_2 & \\\\ -//! -& & & & S\_4 \\\\ -//! =& S\_5 & S\_6 & S\_8 & S\_9 -//! \end{matrix} -//! $$ -//! -//! results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for -//! \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute -//! are then -//! -//! $$ -//! \begin{aligned} -//! X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.00) \\\\ -//! Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ -//! Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ -//! T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 2.00) -//! \end{aligned} -//! $$ -//! -//! which are too large. However, if we flip the sign of \\( S\_4 = -//! S\_0\^2 \\) during squaring, so that we output \\(S\_4' = -S\_4 -//! \pmod p\\), then we can compute -//! -//! $$ -//! \begin{matrix} -//! & S\_1 & S\_1 & S\_1 & S\_1 \\\\ -//! +& S\_2 & & & S\_2 \\\\ -//! +& & & S\_3 & \\\\ -//! +& & & S\_3 & \\\\ -//! +& & & & S\_4' \\\\ -//! +& & 2p & 2p & \\\\ -//! -& & S\_2 & S\_2 & \\\\ -//! =& S\_5 & S\_6 & S\_8 & S\_9 -//! \end{matrix} -//! $$ -//! -//! resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for -//! \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute -//! are then -//! -//! $$ -//! \begin{aligned} -//! X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.59) \\\\ -//! Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ -//! Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ -//! T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 1.59) -//! \end{aligned} -//! $$ -//! -//! whose right-hand sides are all bounded with \\( b < 1.75 \\) and -//! whose left-hand sides are all bounded with \\( b < 2.5 \\). -//! -//! # Field element representation -//! -//! The field element representation is oriented around the AVX2 -//! `vpmuluqdq` instruction, which multiplies the low 32 bits of each -//! 64-bit lane of each operand to produce a 64-bit result. -//! -//! ```text,no_run -//! (a1 ?? b1 ?? c1 ?? d1 ??) -//! (a2 ?? b2 ?? c2 ?? d2 ??) -//! -//! (a1*a2 b1*b2 c1*c2 d1*d2) -//! ``` -//! -//! To unpack 32-bit values into 64-bit lanes for use in multiplication -//! it would be convenient to use the `vpunpck[lh]dq` instructions, -//! which unpack and interleave the low and high 32-bit lanes of two -//! source vectors. -//! However, the AVX2 versions of these instructions are designed to -//! operate only within 128-bit lanes of the 256-bit vectors, so that -//! interleaving the low lanes of `(a0 b0 c0 d0 a1 b1 c1 d1)` with zero -//! gives `(a0 00 b0 00 a1 00 b1 00)`. Instead, we pre-shuffle the data -//! layout as `(a0 b0 a1 b1 c0 d0 c1 d1)` so that we can unpack the -//! "low" and "high" parts as -//! -//! ```text,no_run -//! (a0 00 b0 00 c0 00 d0 00) -//! (a1 00 b1 00 c1 00 d1 00) -//! ``` -//! -//! The data layout for a vector of four field elements \\( (a,b,c,d) -//! \\) with limbs \\( a_0, a_1, \ldots, a_9 \\) is as `[u32x8; 5]` in -//! the form -//! -//! ```text,no_run -//! (a0 b0 a1 b1 c0 d0 c1 d1) -//! (a2 b2 a3 b3 c2 d2 c3 d3) -//! (a4 b4 a5 b5 c4 d4 c5 d5) -//! (a6 b6 a7 b7 c6 d6 c7 d7) -//! (a8 b8 a9 b9 c8 d8 c9 d9) -//! ``` -//! -//! Since this breaks cleanly into two 128-bit lanes, it may be possible -//! to adapt it to 128-bit vector instructions such as NEON without too -//! much difficulty. -//! -//! Going the other direction, to extend this to AVX512, we could either -//! run two point operations in parallel in lower and upper halves of -//! the registers, or use 2-way parallelism within a field operation. -//! -//! We don't attempt to use AVX2 for serial field element computations -//! such as inversion, since wherever we have AVX2 we also have `mulx`. -//! However, it might be useful for batched inverse square-root -//! computations, which can't be batched in the same way inversions can. -//! -//! # Implementation details -//! -//! The implementation uses the unstable `stdsimd` crate to provide AVX2 -//! intrinsics, and the code is not yet cleanly factored between the -//! field element parts and the point parts. -//! -//! When compiling with AVX512VL, LLVM is able to use the extra -//! `ymm16..ymm31` registers to reduce register pressure, and avoid -//! spills during field multiplication. This gives a small but -//! noticeable speedup. -//! -//! The addition and subtraction steps involve masking, to apply -//! operations to a single lane of the vector. AVX512VL extends the -//! predication features of AVX512 to AVX2 code and would probably be -//! beneficial. Unfortunately, LLVM is currently unable to lower `op + -//! blend` into an AVX512VL masked operation. However, the explicitly -//! masked versions of the intrinsics seem to produce the same LLVM IR -//! as an `op + blend`, so hopefully this will improve as the AVX512 -//! support in LLVM improves. -//! -//! When used for constant-time variable-base scalar multiplication, -//! this strategy (using AVX2) gives a significant speedup over the -//! serial implementation (using the \\(64 \times 64\\) multiplier) of -//! approximately 1.6x for Skylake-X with `target_cpu=skylake` (using AVX2), of -//! approximately 1.8x for Skylake-X with `target_cpu=skylake-avx512` (using the extra -//! `ymm16..ymm31` registers from AVX512VL), and of approximately 1.0x -//! for Ryzen (which implements AVX2 at half rate). -//! -//! When used for variable-time double-base scalar multiplication -//! \\( aA + bB \\) for fixed \\(B\\) (as in, e.g., signature verification), -//! this strategy provides a 1.4x speedup on Skylake-X over the same -//! operation as implemented in `ed25519-donna`, the fastest -//! production-quality Ed25519 implementation. -//! -//! [sandy2x]: https://eprint.iacr.org/2015/943.pdf -//! [avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28 -//! [hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf - +// See the comment above the ristretto::notes module. +#![cfg_attr(all(feature = "nightly", feature="precomputed_tables"), doc(include = "../docs/avx2-notes.md"))] pub(crate) mod field; diff --git a/src/ristretto.rs b/src/ristretto.rs index d30a3a2..2d89f84 100644 --- a/src/ristretto.rs +++ b/src/ristretto.rs @@ -159,6 +159,8 @@ // generating the lookup tables (in which case we're relative to the // location of build.rs, not lib.rs, so the markdown file appears // missing). +// +// This hack is also used in the avx2 notes. #[cfg_attr(all(feature = "nightly", feature="precomputed_tables"), doc(include = "../docs/ristretto-notes.md"))] mod notes { } From 1464c4101dfb5558d3385c16a81b48614b6d1a7e Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Sun, 8 Apr 2018 16:29:54 -0700 Subject: [PATCH 7/9] change LSB to least significant bit --- src/scalar.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/scalar.rs b/src/scalar.rs index c6b9917..d41d22b 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -534,7 +534,8 @@ impl Scalar { /// \\( \bar x \\) with \\( \bar x \equiv x \pmod{2^w} \\) and /// \\( -2^{w-1} \leq \bar x < 2^w \\). /// - /// We implement this by scanning across the bits of \\(k\\) from LSB to MSB. + /// We implement this by scanning across the bits of \\(k\\) from + /// least-significant bit to most-significant-bit. /// Write the bits of \\(k\\) as /// $$ /// k = \sum\_{i=0}\^m k\_i 2^i, From 67ba20183579a09de23b247e7dceef2095a366f0 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Sun, 8 Apr 2018 16:59:57 -0700 Subject: [PATCH 8/9] Update AVX2 documentation --- docs/avx2-notes.md | 270 ++++++++++++++++++++++----------------------- src/backend/mod.rs | 1 - 2 files changed, 135 insertions(+), 136 deletions(-) diff --git a/docs/avx2-notes.md b/docs/avx2-notes.md index 2e70d82..fb8085d 100644 --- a/docs/avx2-notes.md +++ b/docs/avx2-notes.md @@ -1,7 +1,6 @@ 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 @@ -20,9 +19,26 @@ Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of multiplication and squaring of generic field elements and \\(\mathbf D\\) represents the cost of multiplication by a curve constant. -Currently, this implementation uses only the first two algorithms. +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. -# Parallel formulas +The notes below describe a tweak to the \\( 2\mathbf M + 1\mathbf D \\) +unified addition formulas to give \\( 2\mathbf M \\) readdition with +\\(1\mathbf D\\) precomputation, and a tweak to the doubling formulas to +avoid an extra reduction. These tweaked formulas are the ones used by +the `avx2` backend of `curve25519-dalek`. + +# Parallel formulas in HWCD'08 The doubling formula is presented in the HWCD paper as follows: @@ -46,25 +62,17 @@ and the unified addition algorithm is presented as follows: Here \\( k = 2d \\) is a curve constant. -# Implementation strategy - -For a software implementation, each "processor"'s operations are too +For a software implementation, each processor's operations are too low-latency to parallelize across threads. However, the main cost -is in the multiplication and squaring steps, which share a single -instruction. +is in the multiplication and squaring steps, which are uniform, while +the divergent steps involve inexpensive additions and subtractions. -Our strategy is to implement 4-wide multiplication and squaring -using one 64-bit AVX2 lane for each field element. Field elements -are represented in the usual way as 10 `u32` limbs in radix -\\(25.5\\) (i.e., alternating between \\(2\^{26}\\) for even limbs -and \\(2\^{25}\\) for odd limbs). This has the effect that passing -between the parallel 32-bit AVX2 representation and the serial -64-bit representation amounts to regrouping digits. +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 addition and subtraction steps are done largely serially, using -masking to handle the instruction divergence. The remaining -obstacle to parallelism is the multiplication by the curve constant -\\(k = 2d\\). In the Curve25519 case, this is +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. $$ @@ -95,51 +103,18 @@ prevent accelerating Curve25519, so we don't make this choice. Instead, we just negate one lane, and move the \\(1 \mathbf D\\) into precomputation (see below). -The 4-wide formulas of the HWCD paper do not seem to have been -implemented using SIMD before. The HWCD paper also describes and -analyzes a 2-wide variant of the Montgomery ladder (for comparison -with parallel Edwards formulas); this strategy was used in 2015 by -Tung Chou's `sandy2x` implementation, which used a 2-wide field -implementation in 128-bit vector registers. - -Curiously, however, although the [`sandy2x` paper][sandy2x] also -implements Edwards arithmetic, and cites the HWCD paper, it doesn't -mention or discuss the parallel formulas from HWCD, or that the -2-wide Montgomery formulas it uses were previously published there. -There is also a 2015 paper by Hernández and López on using AVX2 for -the X25519 Montgomery ladder, but neither the paper nor the code are -publicly available, and it apparently gives only a [slight -speedup][avx2trac], suggesting that it also overlooked the -HWCD formulas. - -HWCD also suggest using a mixed representation, passing between \\( -\mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\) -"projective" coordinates, where doubling is slightly cheaper (saving -about \\(\mathbf 1M\\). This approach is used for the -non-vectorized `u32` and `u64` backends, and more -details on the different coordinate systems can be found in the -`curve_models` module documentation. - -This optimization is not compatible with the parallel formulas, which are -therefore slightly less efficient when counting the total number of -field multiplications and squarings. In particular, vectorized doublings -are less efficient than serial doublings. -In addition, the parallel formulas can only use a \\( 32 \times 32 -\rightarrow 64 \\)-bit integer multiplier, so the speedup from -vectorization must overcome the disadvantage of losing the \\( 64 -\times 64 \rightarrow 128\\)-bit (serial) integer multiplier. - # Tweaked formulas After tweaking the formulas as described above, we obtain the following. To avoid confusion with the original HWCD formulas, temporary variables are named \\(S\\) instead of \\(R\\) and are in -static single-assignment (SSA) form. +static single-assignment form. ## Addition -To add points \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and \\(P_2 = (X_2 -: Y_2 : Z_2 : T_2 ) \\), we compute +This implementation only implements readdition, but the tweaked addition +formulas are described first. To add points \\(P_1 = (X_1 : Y_1 : Z_1 : +T_1) \\) and \\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\), we compute $$ \begin{aligned} @@ -293,6 +268,81 @@ $$ to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = [2]P\_1 \\). +Unlike the (re)addition formulas, the divergent parts of these formulas +are less nice. However, with some careful bounds-juggling, it is +possible to implement them without inserting extra carry chains, as +described below. + +# Field element representation + +Our strategy is to implement 4-wide multiplication and squaring by +wordslicing, using one 64-bit AVX2 lane for each field element. Field +elements are represented in the usual way as 10 `u32` limbs in radix +\\(25.5\\) (i.e., alternating between \\(2\^{26}\\) for even limbs and +\\(2\^{25}\\) for odd limbs). This has the effect that passing between +the parallel 32-bit AVX2 representation and the serial 64-bit +representation (which uses radix \\(2^{51}\\)) amounts to regrouping +digits. + +The field element representation is oriented around the AVX2 +`vpmuluqdq` instruction, which multiplies the low 32 bits of each +64-bit lane of each operand to produce a 64-bit result. + +```text,no_run +(a1 ?? b1 ?? c1 ?? d1 ??) +(a2 ?? b2 ?? c2 ?? d2 ??) + +(a1*a2 b1*b2 c1*c2 d1*d2) +``` + +To unpack 32-bit values into 64-bit lanes for use in multiplication +it would be convenient to use the `vpunpck[lh]dq` instructions, +which unpack and interleave the low and high 32-bit lanes of two +source vectors. +However, the AVX2 versions of these instructions are designed to +operate only within 128-bit lanes of the 256-bit vectors, so that +interleaving the low lanes of `(a0 b0 c0 d0 a1 b1 c1 d1)` with zero +gives `(a0 00 b0 00 a1 00 b1 00)`. Instead, we pre-shuffle the data +layout as `(a0 b0 a1 b1 c0 d0 c1 d1)` so that we can unpack the +"low" and "high" parts as + +```text,no_run +(a0 00 b0 00 c0 00 d0 00) +(a1 00 b1 00 c1 00 d1 00) +``` + +The data layout for a vector of four field elements \\( (a,b,c,d) +\\) with limbs \\( a_0, a_1, \ldots, a_9 \\) is as `[u32x8; 5]` in +the form + +```text,no_run +(a0 b0 a1 b1 c0 d0 c1 d1) +(a2 b2 a3 b3 c2 d2 c3 d3) +(a4 b4 a5 b5 c4 d4 c5 d5) +(a6 b6 a7 b7 c6 d6 c7 d7) +(a8 b8 a9 b9 c8 d8 c9 d9) +``` + +Since this breaks cleanly into two 128-bit lanes, it may be possible +to adapt it to 128-bit vector instructions such as NEON without too +much difficulty. Going the other direction, to extend this to AVX512, +we could either run two point operations in parallel in lower and upper +halves of the registers, or use 2-way parallelism within a field operation. + +# Handling the Doubling Formulas + +The non-parallel portion of the doubling formulas is + +$$ +\begin{aligned} +S\_5 &\gets S\_1 + S\_2 \\\\ +S\_6 &\gets S\_1 - S\_2 \\\\ +S\_7 &\gets 2S\_3 \\\\ +S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\ +S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4 +\end{aligned} +$$ + Performing too many intermediate additions and subtractions grows the bounds beyond what is allowed as input to multiplication, forcing an extra carry pass. However, it is just possible to avoid @@ -300,8 +350,8 @@ this by rearranging signs. Assume that the bounds on the limbs of each field element are parameterized by \\( b \in \mathbb R \\) representing the excess -bits, so that each limb is bounded by either \\( 2\^{25} \\) or \\( -2\^{26} \\). +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 \\) @@ -371,93 +421,43 @@ $$ whose right-hand sides are all bounded with \\( b < 1.75 \\) and whose left-hand sides are all bounded with \\( b < 2.5 \\). -# Field element representation +# Comparison to non-vectorized formulas -The field element representation is oriented around the AVX2 -`vpmuluqdq` instruction, which multiplies the low 32 bits of each -64-bit lane of each operand to produce a 64-bit result. +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. -```text,no_run -(a1 ?? b1 ?? c1 ?? d1 ??) -(a2 ?? b2 ?? c2 ?? d2 ??) +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. -(a1*a2 b1*b2 c1*c2 d1*d2) -``` - -To unpack 32-bit values into 64-bit lanes for use in multiplication -it would be convenient to use the `vpunpck[lh]dq` instructions, -which unpack and interleave the low and high 32-bit lanes of two -source vectors. -However, the AVX2 versions of these instructions are designed to -operate only within 128-bit lanes of the 256-bit vectors, so that -interleaving the low lanes of `(a0 b0 c0 d0 a1 b1 c1 d1)` with zero -gives `(a0 00 b0 00 a1 00 b1 00)`. Instead, we pre-shuffle the data -layout as `(a0 b0 a1 b1 c0 d0 c1 d1)` so that we can unpack the -"low" and "high" parts as - -```text,no_run -(a0 00 b0 00 c0 00 d0 00) -(a1 00 b1 00 c1 00 d1 00) -``` - -The data layout for a vector of four field elements \\( (a,b,c,d) -\\) with limbs \\( a_0, a_1, \ldots, a_9 \\) is as `[u32x8; 5]` in -the form - -```text,no_run -(a0 b0 a1 b1 c0 d0 c1 d1) -(a2 b2 a3 b3 c2 d2 c3 d3) -(a4 b4 a5 b5 c4 d4 c5 d5) -(a6 b6 a7 b7 c6 d6 c7 d7) -(a8 b8 a9 b9 c8 d8 c9 d9) -``` - -Since this breaks cleanly into two 128-bit lanes, it may be possible -to adapt it to 128-bit vector instructions such as NEON without too -much difficulty. - -Going the other direction, to extend this to AVX512, we could either -run two point operations in parallel in lower and upper halves of -the registers, or use 2-way parallelism within a field operation. - -We don't attempt to use AVX2 for serial field element computations -such as inversion, since wherever we have AVX2 we also have `mulx`. -However, it might be useful for batched inverse square-root -computations, which can't be batched in the same way inversions can. - -# Implementation details - -The implementation uses the unstable `stdsimd` crate to provide AVX2 -intrinsics, and the code is not yet cleanly factored between the -field element parts and the point parts. +In addition, the parallel formulas can only use a \\( 32 \times 32 +\rightarrow 64 \\)-bit integer multiplier, so the speedup from +vectorization must overcome the disadvantage of losing the \\( 64 +\times 64 \rightarrow 128\\)-bit (serial) integer multiplier. When compiling with AVX512VL, LLVM is able to use the extra `ymm16..ymm31` registers to reduce register pressure, and avoid spills during field multiplication. This gives a small but noticeable speedup. -The addition and subtraction steps involve masking, to apply -operations to a single lane of the vector. AVX512VL extends the -predication features of AVX512 to AVX2 code and would probably be -beneficial. Unfortunately, LLVM is currently unable to lower `op + -blend` into an AVX512VL masked operation. However, the explicitly -masked versions of the intrinsics seem to produce the same LLVM IR -as an `op + blend`, so hopefully this will improve as the AVX512 -support in LLVM improves. +Another concern with AVX2 is that currently-available Intel processors +(particularly Skylake and Skylake-X microarchitectures) perform thermal +throttling when using wide vector instructions. For a mixed workload, +where point operations are interspersed with other tasks, this can +reduce overall performance. This probably means that this +implementation is not suitable for basic applications, like signatures, +but could still be worthwhile for complex applications, like +zero-knowledge proofs, which do enough work to make it worthwhile. -When used for constant-time variable-base scalar multiplication, -this strategy (using AVX2) gives a significant speedup over the -serial implementation (using the \\(64 \times 64\\) multiplier) of -approximately 1.6x for Skylake-X with `target_cpu=skylake` (using AVX2), of -approximately 1.8x for Skylake-X with `target_cpu=skylake-avx512` (using the extra -`ymm16..ymm31` registers from AVX512VL), and of approximately 1.0x -for Ryzen (which implements AVX2 at half rate). - -When used for variable-time double-base scalar multiplication -\\( aA + bB \\) for fixed \\(B\\) (as in, e.g., signature verification), -this strategy provides a 1.4x speedup on Skylake-X over the same -operation as implemented in `ed25519-donna`, the fastest -production-quality Ed25519 implementation. +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 diff --git a/src/backend/mod.rs b/src/backend/mod.rs index 847b1fe..bfa2a66 100644 --- a/src/backend/mod.rs +++ b/src/backend/mod.rs @@ -27,7 +27,6 @@ pub mod u32; #[cfg(feature="radix_51")] pub mod u64; -/// Code using AVX2. #[cfg(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2")))] pub mod avx2; From e492ac63d0eff26ffb36935aa3a7069b31b2b654 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Mon, 9 Apr 2018 09:50:03 -0700 Subject: [PATCH 9/9] Bump version to 0.16.3 --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index a20ee97..ee1d0b5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "curve25519-dalek" -version = "0.16.2" +version = "0.16.3" authors = ["Isis Lovecruft ", "Henry de Valence "] readme = "README.md"