From 42430f1a1f8bed1e22afeb6a7d0f7c0e1637c34e Mon Sep 17 00:00:00 2001 From: Tony Arcieri Date: Sat, 2 Jun 2018 11:17:10 -0700 Subject: [PATCH 01/27] Upgrade rand to (non-pre)release version 0.5 It was previously using pre-release version 0.5.0-pre.2 --- Cargo.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index dca1a7d..1103faf 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -41,7 +41,7 @@ harness = false # match exactly, since the build.rs uses the crate itself as a library. [dependencies] -rand = { version = "0.5.0-pre.2", default-features = false } +rand = { version = "0.5.0", default-features = false } byteorder = { version = "1", default-features = false } digest = "0.7" generic-array = "0.9" @@ -50,7 +50,7 @@ subtle = { version = "0.6", features = ["generic-impls"], default-features = fal serde = { version = "1.0", optional = true } [build-dependencies] -rand = { version = "0.5.0-pre.2", default-features = false } +rand = { version = "0.5.0", default-features = false } byteorder = { version = "1", default-features = false } digest = "0.7" generic-array = "0.9" From bd1e3c5f3e227b980768074fcbddb4d890a134ac Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:03:11 -0700 Subject: [PATCH 02/27] some rustfmt changes --- src/backend/avx2/edwards.rs | 84 +++++++++------- src/backend/avx2/field.rs | 186 ++++++++++++++++++++++++------------ 2 files changed, 178 insertions(+), 92 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index e317b65..702c64a 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -14,19 +14,19 @@ #![allow(bad_style)] use core::convert::From; -use core::ops::{Add, Sub, Neg}; +use core::ops::{Add, Neg, Sub}; -use core::simd::{IntoBits, u32x8}; +use core::simd::{u32x8, IntoBits}; -use subtle::ConditionallyAssignable; use subtle::Choice; +use subtle::ConditionallyAssignable; use edwards; use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use traits::Identity; -use backend::avx2::field::{D_LANES, Lanes, FieldElement32x4}; +use backend::avx2::field::{FieldElement32x4, Lanes, D_LANES}; use backend::avx2; @@ -43,7 +43,12 @@ impl From for ExtendedPoint { impl From for edwards::EdwardsPoint { fn from(P: ExtendedPoint) -> edwards::EdwardsPoint { let tmp = P.0.split(); - edwards::EdwardsPoint{X: tmp[0], Y: tmp[1], Z: tmp[2], T: tmp[3]} + edwards::EdwardsPoint { + X: tmp[0], + Y: tmp[1], + Z: tmp[2], + T: tmp[3], + } } } @@ -62,7 +67,7 @@ impl Default for ExtendedPoint { impl Identity for ExtendedPoint { fn identity() -> ExtendedPoint { ExtendedPoint(FieldElement32x4([ - u32x8::new(0,1,0,0,1,0,0,0), + u32x8::new(0, 1, 0, 0, 1, 0, 0, 0), u32x8::splat(0), u32x8::splat(0), u32x8::splat(0), @@ -74,9 +79,9 @@ impl Identity for ExtendedPoint { impl ExtendedPoint { pub fn double(&self) -> ExtendedPoint { unsafe { + use core::arch::x86_64::_mm256_blend_epi32; use core::arch::x86_64::_mm256_permute2x128_si256; use core::arch::x86_64::_mm256_permutevar8x32_epi32; - use core::arch::x86_64::_mm256_blend_epi32; use core::arch::x86_64::_mm256_shuffle_epi32; let P = &self.0; @@ -140,17 +145,22 @@ impl ExtendedPoint { // - | | S2 | S2 | | // ======================= // S5 S6 S8 S9 - // for i in 0..5 { let zero = u32x8::splat(0).into_bits(); let S1: u32x8 = _mm256_permutevar8x32_epi32(t1.0[i].into_bits(), c0).into_bits(); let S2: u32x8 = _mm256_permutevar8x32_epi32(t1.0[i].into_bits(), c1).into_bits(); - let S3_2: u32x8 = _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into_bits(), 0b01010000).into_bits(); + let S3_2: u32x8 = + _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into_bits(), 0b01010000) + .into_bits(); // tmp0 = (0 0 2*S3 -S4) - let tmp0: u32x8 = _mm256_blend_epi32(S3_2.into_bits(), t1.0[i].into_bits(), 0b10100000).into_bits(); + let tmp0: u32x8 = + _mm256_blend_epi32(S3_2.into_bits(), t1.0[i].into_bits(), 0b10100000) + .into_bits(); t0.0[i] = (avx2::constants::P_TIMES_2_MASKED.0[i] + tmp0) + S1; - let S2_pos: u32x8 = _mm256_blend_epi32(zero, S2.into_bits(), 0b10100101).into_bits(); - let S2_neg: u32x8 = _mm256_blend_epi32(S2.into_bits(), zero, 0b10100101).into_bits(); + let S2_pos: u32x8 = + _mm256_blend_epi32(zero, S2.into_bits(), 0b10100101).into_bits(); + let S2_neg: u32x8 = + _mm256_blend_epi32(S2.into_bits(), zero, 0b10100101).into_bits(); t0.0[i] = t0.0[i] + S2_pos; t0.0[i] = t0.0[i] - S2_neg; } @@ -290,7 +300,7 @@ impl<'a> From<&'a edwards::EdwardsPoint> for LookupTable { let P = ExtendedPoint::from(*point); let mut points = [CachedPoint::from(P); 8]; for i in 0..7 { - points[i+1] = (&P + &points[i]).into(); + points[i + 1] = (&P + &points[i]).into(); } LookupTable(points) } @@ -335,23 +345,23 @@ mod test { macro_rules! print_var { ($x:ident) => { println!("{} = {:?}", stringify!($x), $x.to_bytes()); - } + }; } - let S0 = &Y1 - &X1; // R1 - let S1 = &Y1 + &X1; // R3 - let S2 = &Y2 - &X2; // R2 - let S3 = &Y2 + &X2; // R4 + let S0 = &Y1 - &X1; // R1 + let S1 = &Y1 + &X1; // R3 + let S2 = &Y2 - &X2; // R2 + let S3 = &Y2 + &X2; // R4 print_var!(S0); print_var!(S1); print_var!(S2); print_var!(S3); println!(""); - let S4 = &S0 * &S2; // R5 = R1 * R2 - let S5 = &S1 * &S3; // R6 = R3 * R4 - let S6 = &Z1 * &Z2; // R8 - let S7 = &T1 * &T2; // R7 + let S4 = &S0 * &S2; // R5 = R1 * R2 + let S5 = &S1 * &S3; // R6 = R3 * R4 + let S6 = &Z1 * &Z2; // R8 + let S7 = &T1 * &T2; // R7 print_var!(S4); print_var!(S5); print_var!(S6); @@ -362,8 +372,8 @@ mod test { let S9 = &S5 * &FieldElement64([ 121666,0,0,0,0]); // R6 let S10 = &S6 * &FieldElement64([2*121666,0,0,0,0]); // R8 let S11 = &S7 * &(-&FieldElement64([2*121665,0,0,0,0])); // R7 - print_var!(S8 ); - print_var!(S9 ); + print_var!(S8); + print_var!(S9); print_var!(S10); print_var!(S11); println!(""); @@ -378,12 +388,17 @@ mod test { print_var!(S15); println!(""); - let X3 = &S12 * &S14; // R1 * R2 - let Y3 = &S15 * &S13; // R3 * R4 - let Z3 = &S15 * &S14; // R2 * R3 - let T3 = &S12 * &S13; // R1 * R4 + let X3 = &S12 * &S14; // R1 * R2 + let Y3 = &S15 * &S13; // R3 * R4 + let Z3 = &S15 * &S14; // R2 * R3 + let T3 = &S12 * &S13; // R1 * R4 - edwards::EdwardsPoint{X: X3, Y: Y3, Z: Z3, T: T3} + edwards::EdwardsPoint { + X: X3, + Y: Y3, + Z: Z3, + T: T3, + } } fn addition_test_helper(P: edwards::EdwardsPoint, Q: edwards::EdwardsPoint) { @@ -442,10 +457,10 @@ mod test { macro_rules! print_var { ($x:ident) => { println!("{} = {:?}", stringify!($x), $x.to_bytes()); - } + }; } - let S0 = &X1 + &Y1; // R1 + let S0 = &X1 + &Y1; // R1 print_var!(S0); println!(""); @@ -476,7 +491,12 @@ mod test { let Z3 = &S8 * &S6; let T3 = &S5 * &S9; - edwards::EdwardsPoint{X: X3, Y: Y3, Z: Z3, T: T3} + edwards::EdwardsPoint { + X: X3, + Y: Y3, + Z: Z3, + T: T3, + } } fn doubling_test_helper(P: edwards::EdwardsPoint) { diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 96b4baf..5dd8b3a 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -24,11 +24,11 @@ pub const D_LANES64: u8 = 0b11_00_00_00; pub const ALL_LANES: u8 = A_LANES | B_LANES | C_LANES | D_LANES; -use core::ops::Mul; -use core::simd::{IntoBits, u32x8, i32x8, u64x4}; +use core::ops::{Add, Mul}; +use core::simd::{i32x8, u32x8, u64x4, IntoBits}; +use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO}; use backend::u64::field::FieldElement64; -use backend::avx2::constants::{P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_16_LO, P_TIMES_16_HI}; #[derive(Copy, Clone)] pub enum Lanes { @@ -43,9 +43,23 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { use core::arch::x86_64::_mm256_blend_epi32; match control { - Lanes::AB => _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32).into_bits(), - Lanes::CD => _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32).into_bits(), - Lanes::ALL => _mm256_blend_epi32(x.into_bits(), y.into_bits(), ALL_LANES as i32).into_bits(), + Lanes::C => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), C_LANES as i32).into_bits() + } + Lanes::D => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), D_LANES as i32).into_bits() + } + Lanes::AB => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32) + .into_bits() + } + Lanes::CD => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32) + .into_bits() + } + Lanes::ALL => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), ALL_LANES as i32).into_bits() + } } } } @@ -54,8 +68,8 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { #[derive(Clone, Copy, Debug)] pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]); -use subtle::ConditionallyAssignable; use subtle::Choice; +use subtle::ConditionallyAssignable; impl ConditionallyAssignable for FieldElement32x4 { fn conditional_assign(&mut self, other: &FieldElement32x4, choice: Choice) { @@ -71,7 +85,6 @@ impl FieldElement32x4 { pub(crate) fn split(&self) -> [FieldElement64; 4] { let mut out = [FieldElement64::zero(); 4]; for i in 0..5 { - let a_2i = self.0[i].extract(0) as u64; // let b_2i = self.0[i].extract(1) as u64; // let a_2i_1 = self.0[i].extract(2) as u64; // `. @@ -91,11 +104,11 @@ impl FieldElement32x4 { } pub fn zero() -> FieldElement32x4 { - FieldElement32x4([u32x8::splat(0);5]) + FieldElement32x4([u32x8::splat(0); 5]) } pub fn splat(x: &FieldElement64) -> FieldElement32x4 { - FieldElement32x4::new(x,x,x,x) + FieldElement32x4::new(x, x, x, x) } pub fn new( @@ -130,11 +143,31 @@ impl FieldElement32x4 { pub fn negate_D_lazy(&mut self) { unsafe { use core::arch::x86_64::_mm256_blend_epi32; - self.0[0] = _mm256_blend_epi32(self.0[0].into_bits(), (P_TIMES_2_LO - self.0[0]).into_bits(), D_LANES as i32).into_bits(); - self.0[1] = _mm256_blend_epi32(self.0[1].into_bits(), (P_TIMES_2_HI - self.0[1]).into_bits(), D_LANES as i32).into_bits(); - self.0[2] = _mm256_blend_epi32(self.0[2].into_bits(), (P_TIMES_2_HI - self.0[2]).into_bits(), D_LANES as i32).into_bits(); - self.0[3] = _mm256_blend_epi32(self.0[3].into_bits(), (P_TIMES_2_HI - self.0[3]).into_bits(), D_LANES as i32).into_bits(); - self.0[4] = _mm256_blend_epi32(self.0[4].into_bits(), (P_TIMES_2_HI - self.0[4]).into_bits(), D_LANES as i32).into_bits(); + self.0[0] = _mm256_blend_epi32( + self.0[0].into_bits(), + (P_TIMES_2_LO - self.0[0]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[1] = _mm256_blend_epi32( + self.0[1].into_bits(), + (P_TIMES_2_HI - self.0[1]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[2] = _mm256_blend_epi32( + self.0[2].into_bits(), + (P_TIMES_2_HI - self.0[2]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[3] = _mm256_blend_epi32( + self.0[3].into_bits(), + (P_TIMES_2_HI - self.0[3]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[4] = _mm256_blend_epi32( + self.0[4].into_bits(), + (P_TIMES_2_HI - self.0[4]).into_bits(), + D_LANES as i32, + ).into_bits(); } } @@ -144,11 +177,31 @@ impl FieldElement32x4 { pub fn negate_D(&mut self) { unsafe { use core::arch::x86_64::_mm256_blend_epi32; - self.0[0] = _mm256_blend_epi32(self.0[0].into_bits(), (P_TIMES_16_LO - self.0[0]).into_bits(), D_LANES as i32).into_bits(); - self.0[1] = _mm256_blend_epi32(self.0[1].into_bits(), (P_TIMES_16_HI - self.0[1]).into_bits(), D_LANES as i32).into_bits(); - self.0[2] = _mm256_blend_epi32(self.0[2].into_bits(), (P_TIMES_16_HI - self.0[2]).into_bits(), D_LANES as i32).into_bits(); - self.0[3] = _mm256_blend_epi32(self.0[3].into_bits(), (P_TIMES_16_HI - self.0[3]).into_bits(), D_LANES as i32).into_bits(); - self.0[4] = _mm256_blend_epi32(self.0[4].into_bits(), (P_TIMES_16_HI - self.0[4]).into_bits(), D_LANES as i32).into_bits(); + self.0[0] = _mm256_blend_epi32( + self.0[0].into_bits(), + (P_TIMES_16_LO - self.0[0]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[1] = _mm256_blend_epi32( + self.0[1].into_bits(), + (P_TIMES_16_HI - self.0[1]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[2] = _mm256_blend_epi32( + self.0[2].into_bits(), + (P_TIMES_16_HI - self.0[2]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[3] = _mm256_blend_epi32( + self.0[3].into_bits(), + (P_TIMES_16_HI - self.0[3]).into_bits(), + D_LANES as i32, + ).into_bits(); + self.0[4] = _mm256_blend_epi32( + self.0[4].into_bits(), + (P_TIMES_16_HI - self.0[4]).into_bits(), + D_LANES as i32, + ).into_bits(); } self.reduce32(); } @@ -156,11 +209,12 @@ impl FieldElement32x4 { /// Given `self = (A,B,C,D)`, set `self = (B,A,C,D)` pub fn swap_AB(&mut self) { unsafe { - use core::arch::x86_64::_mm256_shuffle_epi32; use core::arch::x86_64::_mm256_blend_epi32; + use core::arch::x86_64::_mm256_shuffle_epi32; for i in 0..5 { let swapped = _mm256_shuffle_epi32(self.0[i].into_bits(), 0b10_11_00_01); - self.0[i] = _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b00001111).into_bits(); + self.0[i] = + _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b00001111).into_bits(); } } } @@ -168,11 +222,12 @@ impl FieldElement32x4 { /// Given `self = (A,B,C,D)`, set `self = (A,B,D,C)` pub fn swap_CD(&mut self) { unsafe { - use core::arch::x86_64::_mm256_shuffle_epi32; use core::arch::x86_64::_mm256_blend_epi32; + use core::arch::x86_64::_mm256_shuffle_epi32; for i in 0..5 { let swapped = _mm256_shuffle_epi32(self.0[i].into_bits(), 0b10_11_00_01); - self.0[i] = _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b11110000).into_bits(); + self.0[i] = + _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b11110000).into_bits(); } } } @@ -183,7 +238,7 @@ impl FieldElement32x4 { #[inline(always)] pub fn diff_sum(&mut self, control: Lanes) { unsafe { - use core::arch::x86_64::{_mm256_shuffle_epi32, _mm256_blend_epi32}; + use core::arch::x86_64::{_mm256_blend_epi32, _mm256_shuffle_epi32}; let shuffle = |v: u32x8| -> u32x8 { _mm256_shuffle_epi32(v.into_bits(), 0b10_11_00_01).into_bits() @@ -192,36 +247,41 @@ impl FieldElement32x4 { let x01 = self.0[0]; let x01_shuf = shuffle(x01); let v1 = (x01_shuf + P_TIMES_2_LO) - x01; - let v2 = x01_shuf + x01; - let diffsum01 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); + let v2 = x01_shuf + x01; + let diffsum01 = + _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); self.0[0] = blend_lanes(x01, diffsum01, control); let x23 = self.0[1]; let x23_shuf = shuffle(x23); let v1 = (x23_shuf + P_TIMES_2_HI) - x23; - let v2 = x23_shuf + x23; - let diffsum23 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); + let v2 = x23_shuf + x23; + let diffsum23 = + _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); self.0[1] = blend_lanes(x23, diffsum23, control); let x45 = self.0[2]; let x45_shuf = shuffle(x45); let v1 = (x45_shuf + P_TIMES_2_HI) - x45; - let v2 = x45_shuf + x45; - let diffsum45 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); + let v2 = x45_shuf + x45; + let diffsum45 = + _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); self.0[2] = blend_lanes(x45, diffsum45, control); let x67 = self.0[3]; let x67_shuf = shuffle(x67); let v1 = (x67_shuf + P_TIMES_2_HI) - x67; - let v2 = x67_shuf + x67; - let diffsum67 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); + let v2 = x67_shuf + x67; + let diffsum67 = + _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); self.0[3] = blend_lanes(x67, diffsum67, control); let x89 = self.0[4]; let x89_shuf = shuffle(x89); let v1 = (x89_shuf + P_TIMES_2_HI) - x89; - let v2 = x89_shuf + x89; - let diffsum89 = _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); + let v2 = x89_shuf + x89; + let diffsum89 = + _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); self.0[4] = blend_lanes(x89, diffsum89, control); } } @@ -233,7 +293,7 @@ impl FieldElement32x4 { pub fn scale_by_curve_constants(&mut self) { let mut b = [u64x4::splat(0); 10]; - let consts = u32x8::new(121666, 0, 121666, 0, 2*121666, 0, 2*121665, 0); + let consts = u32x8::new(121666, 0, 121666, 0, 2 * 121666, 0, 2 * 121665, 0); unsafe { use core::arch::x86_64::_mm256_mul_epu32; @@ -263,10 +323,17 @@ impl FieldElement32x4 { } pub fn reduce32(&mut self) { - - let shifts = i32x8::new(26,26,25,25,26,26,25,25); - let masks = u32x8::new((1<<26)-1, (1<<26)-1, (1<<25)-1, (1<<25)-1, - (1<<26)-1, (1<<26)-1, (1<<25)-1, (1<<25)-1); + let shifts = i32x8::new(26, 26, 25, 25, 26, 26, 25, 25); + let masks = u32x8::new( + (1 << 26) - 1, + (1 << 26) - 1, + (1 << 25) - 1, + (1 << 25) - 1, + (1 << 26) - 1, + (1 << 26) - 1, + (1 << 25) - 1, + (1 << 25) - 1, + ); let carry = |v: u32x8| -> u32x8 { unsafe { @@ -319,19 +386,19 @@ impl FieldElement32x4 { pub fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { // These aren't const because splat isn't a const fn - let LOW_25_BITS: u64x4 = u64x4::splat((1<<25)-1); - let LOW_26_BITS: u64x4 = u64x4::splat((1<<26)-1); + let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1); + let LOW_26_BITS: u64x4 = u64x4::splat((1 << 26) - 1); // Carry the value from limb i = 0..8 to limb i+1 let carry = |z: &mut [u64x4; 10], i: usize| { debug_assert!(i < 9); if i % 2 == 0 { // Even limbs have 26 bits - z[i+1] = z[i+1] + (z[i] >> 26); + z[i + 1] = z[i + 1] + (z[i] >> 26); z[i] = z[i] & LOW_26_BITS; } else { // Odd limbs have 25 bits - z[i+1] = z[i+1] + (z[i] >> 25); + z[i + 1] = z[i + 1] + (z[i] >> 25); z[i] = z[i] & LOW_25_BITS; } }; @@ -386,21 +453,21 @@ impl FieldElement32x4 { pub fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { let a: u32x8; let b: u32x8; - let zero = i32x8::new(0,0,0,0,0,0,0,0); + let zero = i32x8::new(0, 0, 0, 0, 0, 0, 0, 0); unsafe { use core::arch::x86_64::_mm256_unpackhi_epi32; use core::arch::x86_64::_mm256_unpacklo_epi32; a = _mm256_unpacklo_epi32(src.into_bits(), zero.into_bits()).into_bits(); b = _mm256_unpackhi_epi32(src.into_bits(), zero.into_bits()).into_bits(); } - (a,b) + (a, b) } #[inline(always)] pub fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { unsafe { - use core::arch::x86_64::_mm256_shuffle_epi32; use core::arch::x86_64::_mm256_blend_epi32; + use core::arch::x86_64::_mm256_shuffle_epi32; // Input: x = (a0, 0, b0, 0, c0, 0, d0) // Input: y = (a1, 0, b1, 0, c1, 0, d1) @@ -429,16 +496,16 @@ impl FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } #[inline(always)] fn m_lo(x: u32x8, y: u32x8) -> u32x8 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } - let v19 = u32x8::new(19,0,19,0,19,0,19,0); + let v19 = u32x8::new(19, 0, 19, 0, 19, 0, 19, 0); let (x0, x1) = unpack_pair(self.0[0]); let (x2, x3) = unpack_pair(self.0[1]); @@ -520,13 +587,13 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } #[inline(always)] fn m_lo(x: u32x8, y: u32x8) -> u32x8 { use core::arch::x86_64::_mm256_mul_epu32; - unsafe { _mm256_mul_epu32(x.into_bits(),y.into_bits()).into_bits() } + unsafe { _mm256_mul_epu32(x.into_bits(), y.into_bits()).into_bits() } } let (x0, x1) = unpack_pair(self.0[0]); @@ -541,7 +608,7 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let (y6, y7) = unpack_pair(_rhs.0[3]); let (y8, y9) = unpack_pair(_rhs.0[4]); - let v19 = u32x8::new(19,0,19,0,19,0,19,0); + let v19 = u32x8::new(19, 0, 19, 0, 19, 0, 19, 0); let y1_19 = m_lo(v19, y1); // This fits in a u32 let y2_19 = m_lo(v19, y2); // iff 26 + b + lg(19) < 32 @@ -585,10 +652,10 @@ mod test { x.scale_by_curve_constants(); let xs = x.split(); - assert_eq!(xs[0], FieldElement64([ 121666,0,0,0,0])); - assert_eq!(xs[1], FieldElement64([ 121666,0,0,0,0])); - assert_eq!(xs[2], FieldElement64([2*121666,0,0,0,0])); - assert_eq!(xs[3], FieldElement64([2*121665,0,0,0,0])); + assert_eq!(xs[0], FieldElement64([121666, 0, 0, 0, 0])); + assert_eq!(xs[1], FieldElement64([121666, 0, 0, 0, 0])); + assert_eq!(xs[2], FieldElement64([2 * 121666, 0, 0, 0, 0])); + assert_eq!(xs[3], FieldElement64([2 * 121665, 0, 0, 0, 0])); } #[test] @@ -636,7 +703,6 @@ mod test { assert_eq!(result[3], -&(&x3 * &x3)); } - #[test] fn multiply_vs_serial() { let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]); @@ -666,7 +732,7 @@ mod test { let src = vec.0[0]; - let (a,b) = unpack_pair(src); + let (a, b) = unpack_pair(src); let expected_a = u32x8::new(10000, 0, 10100, 0, 10200, 0, 10300, 0); let expected_b = u32x8::new(10001, 0, 10101, 0, 10201, 0, 10301, 0); @@ -674,7 +740,7 @@ mod test { assert_eq!(a, expected_a); assert_eq!(b, expected_b); - let expected_src = repack_pair(a,b); + let expected_src = repack_pair(a, b); assert_eq!(src, expected_src); } From c8dc2a6418bffaa93c07be82ea008a8f2f0120dd Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:07:18 -0700 Subject: [PATCH 03/27] Implement addition for FieldElement32x4 --- src/backend/avx2/field.rs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 5dd8b3a..6751b81 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -581,6 +581,20 @@ impl FieldElement32x4 { } } +impl Add for FieldElement32x4 { + type Output = FieldElement32x4; + #[inline] + fn add(self, rhs: FieldElement32x4) -> FieldElement32x4 { + FieldElement32x4([ + self.0[0] + rhs.0[0], + self.0[1] + rhs.0[1], + self.0[2] + rhs.0[2], + self.0[3] + rhs.0[3], + self.0[4] + rhs.0[4], + ]) + } +} + impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { type Output = FieldElement32x4; fn mul(self, _rhs: &'b FieldElement32x4) -> FieldElement32x4 { From 14ce6d3da640114628c49b54b27ac07ae9280a4c Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:07:33 -0700 Subject: [PATCH 04/27] Add a shuffling abstraction for FieldElement32x4 --- src/backend/avx2/edwards.rs | 83 ++++++++++--------------------------- src/backend/avx2/field.rs | 58 ++++++++++++++++++++++++++ 2 files changed, 81 insertions(+), 60 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 702c64a..be154dd 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -26,7 +26,7 @@ use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use traits::Identity; -use backend::avx2::field::{FieldElement32x4, Lanes, D_LANES}; +use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle, D_LANES}; use backend::avx2; @@ -96,41 +96,21 @@ impl ExtendedPoint { // and then adding. // Set t0 = (X1 Y1 X1 Y1) - t0.0[0] = _mm256_permute2x128_si256(P.0[0].into_bits(), P.0[0].into_bits(), 0b0000_0000).into_bits(); - t0.0[1] = _mm256_permute2x128_si256(P.0[1].into_bits(), P.0[1].into_bits(), 0b0000_0000).into_bits(); - t0.0[2] = _mm256_permute2x128_si256(P.0[2].into_bits(), P.0[2].into_bits(), 0b0000_0000).into_bits(); - t0.0[3] = _mm256_permute2x128_si256(P.0[3].into_bits(), P.0[3].into_bits(), 0b0000_0000).into_bits(); - t0.0[4] = _mm256_permute2x128_si256(P.0[4].into_bits(), P.0[4].into_bits(), 0b0000_0000).into_bits(); + let mut t0 = P.shuffle(Shuffle::ABAB); // Set t1 = (Y1 X1 Y1 X1) - t1.0[0] = _mm256_shuffle_epi32(t0.0[0].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[1] = _mm256_shuffle_epi32(t0.0[1].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[2] = _mm256_shuffle_epi32(t0.0[2].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[3] = _mm256_shuffle_epi32(t0.0[3].into_bits(), 0b10_11_00_01).into_bits(); - t1.0[4] = _mm256_shuffle_epi32(t0.0[4].into_bits(), 0b10_11_00_01).into_bits(); - - // Set t0 = (X1+Y1 X1+Y1 X1+Y1 X1+Y1) - t0.0[0] = t0.0[0] + t1.0[0]; - t0.0[1] = t0.0[1] + t1.0[1]; - t0.0[2] = t0.0[2] + t1.0[2]; - t0.0[3] = t0.0[3] + t1.0[3]; - t0.0[4] = t0.0[4] + t1.0[4]; + let mut t1 = t0.shuffle(Shuffle::BADC); // Set t0 = (X1 Y1 Z1 X1+Y1) - // why does this intrinsic take an i32 for the imm8 ??? - t0.0[0] = _mm256_blend_epi32(P.0[0].into_bits(), t0.0[0].into_bits(), D_LANES as i32).into_bits(); - t0.0[1] = _mm256_blend_epi32(P.0[1].into_bits(), t0.0[1].into_bits(), D_LANES as i32).into_bits(); - t0.0[2] = _mm256_blend_epi32(P.0[2].into_bits(), t0.0[2].into_bits(), D_LANES as i32).into_bits(); - t0.0[3] = _mm256_blend_epi32(P.0[3].into_bits(), t0.0[3].into_bits(), D_LANES as i32).into_bits(); - t0.0[4] = _mm256_blend_epi32(P.0[4].into_bits(), t0.0[4].into_bits(), D_LANES as i32).into_bits(); + t0 = P.blend(&t0 + &t1, Lanes::D); // Set t1 = t0^2, negating the D values t1 = t0.square_and_negate_D(); // Now t1 = (S1 S2 S3 -S4) - let c0 = u32x8::new(0,0,2,2,0,0,2,2).into_bits(); // (ABCD) -> (AAAA) - let c1 = u32x8::new(1,1,3,3,1,1,3,3).into_bits(); // (ABCD) -> (BBBB) + let c0 = u32x8::new(0, 0, 2, 2, 0, 0, 2, 2).into_bits(); // (ABCD) -> (AAAA) + let c1 = u32x8::new(1, 1, 3, 3, 1, 1, 3, 3).into_bits(); // (ABCD) -> (BBBB) // See discussion of bounds in the module-level documentation. // @@ -165,14 +145,8 @@ impl ExtendedPoint { t0.0[i] = t0.0[i] - S2_neg; } - let c0 = u32x8::new(4,0,6,2,4,0,6,2).into_bits(); // (ABCD) -> (CACA) - let c1 = u32x8::new(5,1,7,3,1,5,3,7).into_bits(); // (ABCD) -> (DBBD) - - for i in 0..5 { - let tmp = t0.0[i]; - t0.0[i] = _mm256_permutevar8x32_epi32(tmp.into_bits(), c0).into_bits(); - t1.0[i] = _mm256_permutevar8x32_epi32(tmp.into_bits(), c1).into_bits(); - } + t1 = t0.shuffle(Shuffle::DBBD); + t0 = t0.shuffle(Shuffle::CACA); ExtendedPoint(&t0 * &t1) } @@ -248,38 +222,27 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { /// Uses a slight tweak of the parallel unified formulas of HWCD'08 fn add(self, other: &'b CachedPoint) -> ExtendedPoint { - unsafe { - use core::arch::x86_64::_mm256_permutevar8x32_epi32; + let mut tmp = self.0; - let mut tmp = self.0; + // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) + tmp.diff_sum(Lanes::AB); - // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) - tmp.diff_sum(Lanes::AB); + // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) + tmp = &tmp * &other.0; - // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) - tmp = &tmp * &other.0; + // tmp = (S8 S9 S11 S10) + tmp.swap_CD(); - // tmp = (S8 S9 S11 S10) - tmp.swap_CD(); + // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) + tmp.diff_sum(Lanes::ALL); - // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - tmp.diff_sum(Lanes::ALL); + // set t0 = (S12 S15 S15 S12) + let t0 = tmp.shuffle(Shuffle::ADDA); + // set t1 = (S14 S13 S14 S13) + let t1 = tmp.shuffle(Shuffle::CBCB); - let c0 = u32x8::new(0,5,2,7,5,0,7,2); // (ABCD) -> (ADDA) - let c1 = u32x8::new(4,1,6,3,4,1,6,3); // (ABCD) -> (CBCB) - - // set t0 = (S12 S15 S15 S12) - // set t1 = (S14 S13 S14 S13) - let mut t0 = FieldElement32x4::zero(); - let mut t1 = FieldElement32x4::zero(); - for i in 0..5 { - t0.0[i] = _mm256_permutevar8x32_epi32(tmp.0[i].into_bits(), c0.into_bits()).into_bits(); - t1.0[i] = _mm256_permutevar8x32_epi32(tmp.0[i].into_bits(), c1.into_bits()).into_bits(); - } - - // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) - ExtendedPoint(&t0 * &t1) - } + // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + ExtendedPoint(&t0 * &t1) } } diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 6751b81..d5ac6a8 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -32,6 +32,8 @@ use backend::u64::field::FieldElement64; #[derive(Copy, Clone)] pub enum Lanes { + C, + D, AB, CD, ALL, @@ -64,6 +66,18 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { } } +#[derive(Copy, Clone)] +pub enum Shuffle { + AAAA, + BBBB, + CACA, + DBBD, + ADDA, + CBCB, + ABAB, + BADC, +} + /// A vector of four `FieldElements`, implemented using AVX2. #[derive(Clone, Copy, Debug)] pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]); @@ -103,6 +117,50 @@ impl FieldElement32x4 { out } + #[inline(always)] + pub fn blend(&self, b: FieldElement32x4, control: Lanes) -> FieldElement32x4 { + FieldElement32x4([ + blend_lanes(self.0[0], b.0[0], control), + blend_lanes(self.0[1], b.0[1], control), + blend_lanes(self.0[2], b.0[2], control), + blend_lanes(self.0[3], b.0[3], control), + blend_lanes(self.0[4], b.0[4], control), + ]) + } + + #[inline(always)] + pub fn shuffle(&self, control: Shuffle) -> FieldElement32x4 { + #[inline(always)] + fn shuffle_lanes(x: u32x8, control: Shuffle) -> u32x8 { + unsafe { + use core::arch::x86_64::_mm256_permutevar8x32_epi32; + + let c: u32x8 = match control { + Shuffle::AAAA => u32x8::new(0, 0, 2, 2, 0, 0, 2, 2), + Shuffle::BBBB => u32x8::new(1, 1, 3, 3, 1, 1, 3, 3), + Shuffle::CACA => u32x8::new(4, 0, 6, 2, 4, 0, 6, 2), + Shuffle::DBBD => u32x8::new(5, 1, 7, 3, 1, 5, 3, 7), + Shuffle::ADDA => u32x8::new(0, 5, 2, 7, 5, 0, 7, 2), + Shuffle::CBCB => u32x8::new(4, 1, 6, 3, 4, 1, 6, 3), + Shuffle::ABAB => u32x8::new(0, 1, 2, 3, 0, 1, 2, 3), + Shuffle::BADC => u32x8::new(1, 0, 3, 2, 5, 4, 7, 6), + }; + // Note that this gets turned into a generic LLVM + // shuffle-by-constants, which can be lowered to a simpler + // instruction than a generic permute. + _mm256_permutevar8x32_epi32(x.into_bits(), c.into_bits()).into_bits() + } + } + + FieldElement32x4([ + shuffle_lanes(self.0[0], control), + shuffle_lanes(self.0[1], control), + shuffle_lanes(self.0[2], control), + shuffle_lanes(self.0[3], control), + shuffle_lanes(self.0[4], control), + ]) + } + pub fn zero() -> FieldElement32x4 { FieldElement32x4([u32x8::splat(0); 5]) } From 00d8b6ea4f177a7dcf9ba06c8c08cc0e303e6fdf Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 13 Jun 2018 12:32:42 -0700 Subject: [PATCH 05/27] Change `negate_D`, `negate_D_lazy` to `impl Neg`, `negate_lazy` --- src/backend/avx2/edwards.rs | 9 ++-- src/backend/avx2/field.rs | 101 +++++++++++++----------------------- 2 files changed, 40 insertions(+), 70 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index be154dd..73db3db 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -176,7 +176,7 @@ impl From for CachedPoint { x.scale_by_curve_constants(); // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) - x.negate_D(); + x = x.blend(-x, Lanes::D); CachedPoint(x) } @@ -210,10 +210,9 @@ impl<'a> Neg for &'a CachedPoint { type Output = CachedPoint; fn neg(self) -> CachedPoint { - let mut neg = *self; - neg.0.swap_AB(); - neg.0.negate_D_lazy(); - neg + let mut coords = self.0; + coords.swap_AB(); + CachedPoint(coords.blend(coords.negate_lazy(), Lanes::D)) } } diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index d5ac6a8..c7ee80b 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -24,7 +24,7 @@ pub const D_LANES64: u8 = 0b11_00_00_00; pub const ALL_LANES: u8 = A_LANES | B_LANES | C_LANES | D_LANES; -use core::ops::{Add, Mul}; +use core::ops::{Add, Mul, Neg}; use core::simd::{i32x8, u32x8, u64x4, IntoBits}; use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO}; @@ -195,73 +195,20 @@ impl FieldElement32x4 { return out; } - /// Negate the \\(D\\) variable of \\((A,B,C,D)\\). + /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), without + /// performing a reduction. /// /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. - pub fn negate_D_lazy(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - self.0[0] = _mm256_blend_epi32( - self.0[0].into_bits(), - (P_TIMES_2_LO - self.0[0]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[1] = _mm256_blend_epi32( - self.0[1].into_bits(), - (P_TIMES_2_HI - self.0[1]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[2] = _mm256_blend_epi32( - self.0[2].into_bits(), - (P_TIMES_2_HI - self.0[2]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[3] = _mm256_blend_epi32( - self.0[3].into_bits(), - (P_TIMES_2_HI - self.0[3]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[4] = _mm256_blend_epi32( - self.0[4].into_bits(), - (P_TIMES_2_HI - self.0[4]).into_bits(), - D_LANES as i32, - ).into_bits(); - } - } - - /// Negate the \\(D\\) variable of \\((A,B,C,D)\\). /// - /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. - pub fn negate_D(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - self.0[0] = _mm256_blend_epi32( - self.0[0].into_bits(), - (P_TIMES_16_LO - self.0[0]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[1] = _mm256_blend_epi32( - self.0[1].into_bits(), - (P_TIMES_16_HI - self.0[1]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[2] = _mm256_blend_epi32( - self.0[2].into_bits(), - (P_TIMES_16_HI - self.0[2]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[3] = _mm256_blend_epi32( - self.0[3].into_bits(), - (P_TIMES_16_HI - self.0[3]).into_bits(), - D_LANES as i32, - ).into_bits(); - self.0[4] = _mm256_blend_epi32( - self.0[4].into_bits(), - (P_TIMES_16_HI - self.0[4]).into_bits(), - D_LANES as i32, - ).into_bits(); - } - self.reduce32(); + /// The output limbs are bounded by \\(2p\\). + pub fn negate_lazy(&self) -> FieldElement32x4 { + FieldElement32x4([ + P_TIMES_2_LO - self.0[0], + P_TIMES_2_HI - self.0[1], + P_TIMES_2_HI - self.0[2], + P_TIMES_2_HI - self.0[3], + P_TIMES_2_HI - self.0[4], + ]) } /// Given `self = (A,B,C,D)`, set `self = (B,A,C,D)` @@ -639,6 +586,30 @@ impl FieldElement32x4 { } } +impl Neg for FieldElement32x4 { + type Output = FieldElement32x4; + + /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), and + /// perform a reduction. + /// + /// The input limbs can be any size. + /// + /// The output limbs are freshly reduced. + #[inline] + fn neg(self) -> FieldElement32x4 { + let mut neg = FieldElement32x4([ + P_TIMES_16_LO - self.0[0], + P_TIMES_16_HI - self.0[1], + P_TIMES_16_HI - self.0[2], + P_TIMES_16_HI - self.0[3], + P_TIMES_16_HI - self.0[4], + ]); + neg.reduce32(); + + neg + } +} + impl Add for FieldElement32x4 { type Output = FieldElement32x4; #[inline] From 64cb9998661976b7fe3aeefd8603735555643cc5 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:19:26 -0700 Subject: [PATCH 06/27] Make platform-vector lanes constants private --- src/backend/avx2/edwards.rs | 2 +- src/backend/avx2/field.rs | 18 ++++++++---------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 73db3db..7dd4353 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -26,7 +26,7 @@ use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use traits::Identity; -use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle, D_LANES}; +use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle}; use backend::avx2; diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index c7ee80b..331baf3 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -12,17 +12,15 @@ #![allow(bad_style)] -pub const A_LANES: u8 = 0b0000_0101; -pub const B_LANES: u8 = 0b0000_1010; -pub const C_LANES: u8 = 0b0101_0000; -pub const D_LANES: u8 = 0b1010_0000; +const A_LANES: u8 = 0b0000_0101; +const B_LANES: u8 = 0b0000_1010; +const C_LANES: u8 = 0b0101_0000; +const D_LANES: u8 = 0b1010_0000; -pub const A_LANES64: u8 = 0b00_00_00_11; -pub const B_LANES64: u8 = 0b00_00_11_00; -pub const C_LANES64: u8 = 0b00_11_00_00; -pub const D_LANES64: u8 = 0b11_00_00_00; - -pub const ALL_LANES: u8 = A_LANES | B_LANES | C_LANES | D_LANES; +const A_LANES64: u8 = 0b00_00_00_11; +const B_LANES64: u8 = 0b00_00_11_00; +const C_LANES64: u8 = 0b00_11_00_00; +const D_LANES64: u8 = 0b11_00_00_00; use core::ops::{Add, Mul, Neg}; use core::simd::{i32x8, u32x8, u64x4, IntoBits}; From 30a2b01c0527f3a47005d889d77ac33f0f2c2808 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:19:46 -0700 Subject: [PATCH 07/27] Add more selectors to the Lanes enum --- src/backend/avx2/field.rs | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 331baf3..13a6156 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -34,6 +34,8 @@ pub enum Lanes { D, AB, CD, + AD, + BC, ALL, } @@ -49,6 +51,10 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { Lanes::D => { _mm256_blend_epi32(x.into_bits(), y.into_bits(), D_LANES as i32).into_bits() } + Lanes::AD => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | D_LANES) as i32) + .into_bits() + } Lanes::AB => { _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32) .into_bits() @@ -57,9 +63,15 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32) .into_bits() } - Lanes::ALL => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), ALL_LANES as i32).into_bits() + Lanes::BC => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (B_LANES | C_LANES) as i32) + .into_bits() } + Lanes::ALL => _mm256_blend_epi32( + x.into_bits(), + y.into_bits(), + (A_LANES | B_LANES | C_LANES | D_LANES) as i32, + ).into_bits(), } } } From 794ed5c8e3c2ebed34f3b2f6039a81b3e77c793b Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:20:00 -0700 Subject: [PATCH 08/27] Rewrite the doubling horrorshow --- src/backend/avx2/edwards.rs | 108 +++++++++++++++--------------------- 1 file changed, 45 insertions(+), 63 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 7dd4353..563d3d6 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -78,78 +78,60 @@ impl Identity for ExtendedPoint { impl ExtendedPoint { pub fn double(&self) -> ExtendedPoint { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - use core::arch::x86_64::_mm256_permute2x128_si256; - use core::arch::x86_64::_mm256_permutevar8x32_epi32; - use core::arch::x86_64::_mm256_shuffle_epi32; + // Want to compute (X1 Y1 Z1 X1+Y1). + // Not sure how to do this less expensively than computing + // (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1) + // (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1) + // and then adding. - let P = &self.0; + // Set tmp0 = (X1 Y1 X1 Y1) + let mut tmp0 = self.0.shuffle(Shuffle::ABAB); - let mut t0 = FieldElement32x4::zero(); - let mut t1 = FieldElement32x4::zero(); + // Set tmp1 = (Y1 X1 Y1 X1) + let mut tmp1 = tmp0.shuffle(Shuffle::BADC); - // Want to compute (X1 Y1 Z1 X1+Y1). - // Not sure how to do this less expensively than computing - // (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1) - // (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1) - // and then adding. + // Set tmp0 = (X1 Y1 Z1 X1+Y1) + tmp0 = self.0.blend(tmp0 + tmp1, Lanes::D); - // Set t0 = (X1 Y1 X1 Y1) - let mut t0 = P.shuffle(Shuffle::ABAB); + // Set tmp1 = tmp0^2, negating the D values + tmp1 = tmp0.square_and_negate_D(); + // Now tmp1 = (S1 S2 S3 -S4) - // Set t1 = (Y1 X1 Y1 X1) - let mut t1 = t0.shuffle(Shuffle::BADC); + // See discussion of bounds in the module-level documentation. + // We want to compute + // + // + | S1 | S1 | S1 | S1 | + // + | S2 | | | S2 | + // + | | | S3 | | + // + | | | S3 | | + // + | | | |-S4 | + // + | | 2p | 2p | | + // - | | S2 | S2 | | + // ======================= + // S5 S6 S8 S9 - // Set t0 = (X1 Y1 Z1 X1+Y1) - t0 = P.blend(&t0 + &t1, Lanes::D); + let zero = FieldElement32x4::zero(); + let S_1 = tmp1.shuffle(Shuffle::AAAA); + let S_2 = tmp1.shuffle(Shuffle::BBBB); - // Set t1 = t0^2, negating the D values - t1 = t0.square_and_negate_D(); + // tmp0 = (0, 0, 2S_3, 0) + tmp0 = zero.blend(tmp1 + tmp1, Lanes::C); + // tmp0 = (0, 0, 2S_3, -S_4) + tmp0 = tmp0.blend(tmp1, Lanes::D); + // tmp0 = (S_1, S_1, S_1 + 2S_3, S_1 - S_4) + tmp0 = tmp0 + S_1; + // tmp0 = (S_1 + S_2, S_1, S_1 + 2S_3, S_1 + S_2 - S_4) + tmp0 = tmp0 + zero.blend(S_2, Lanes::AD); + // tmp0 = (S_1 + S_2, S_1 - S_2, S_1 - S_2 + 2S_3, S_1 + S_2 - S_4) + tmp0 = tmp0 + zero.blend(S_2.negate_lazy(), Lanes::BC); + // Now tmp0 = (S_5, S_6, S_8, S_9) - // Now t1 = (S1 S2 S3 -S4) + // Set tmp1 = (S_9, S_6, S_6, S_9) + tmp1 = tmp0.shuffle(Shuffle::DBBD); + // Set tmp1 = (S_8, S_5, S_8, S_5) + tmp0 = tmp0.shuffle(Shuffle::CACA); - let c0 = u32x8::new(0, 0, 2, 2, 0, 0, 2, 2).into_bits(); // (ABCD) -> (AAAA) - let c1 = u32x8::new(1, 1, 3, 3, 1, 1, 3, 3).into_bits(); // (ABCD) -> (BBBB) - - // See discussion of bounds in the module-level documentation. - // - // We want to compute - // - // + | S1 | S1 | S1 | S1 | - // + | S2 | | | S2 | - // + | | | S3 | | - // + | | | S3 | | - // + | | | |-S4 | - // + | | 2p | 2p | | - // - | | S2 | S2 | | - // ======================= - // S5 S6 S8 S9 - for i in 0..5 { - let zero = u32x8::splat(0).into_bits(); - let S1: u32x8 = _mm256_permutevar8x32_epi32(t1.0[i].into_bits(), c0).into_bits(); - let S2: u32x8 = _mm256_permutevar8x32_epi32(t1.0[i].into_bits(), c1).into_bits(); - let S3_2: u32x8 = - _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into_bits(), 0b01010000) - .into_bits(); - // tmp0 = (0 0 2*S3 -S4) - let tmp0: u32x8 = - _mm256_blend_epi32(S3_2.into_bits(), t1.0[i].into_bits(), 0b10100000) - .into_bits(); - t0.0[i] = (avx2::constants::P_TIMES_2_MASKED.0[i] + tmp0) + S1; - let S2_pos: u32x8 = - _mm256_blend_epi32(zero, S2.into_bits(), 0b10100101).into_bits(); - let S2_neg: u32x8 = - _mm256_blend_epi32(S2.into_bits(), zero, 0b10100101).into_bits(); - t0.0[i] = t0.0[i] + S2_pos; - t0.0[i] = t0.0[i] - S2_neg; - } - - t1 = t0.shuffle(Shuffle::DBBD); - t0 = t0.shuffle(Shuffle::CACA); - - ExtendedPoint(&t0 * &t1) - } + ExtendedPoint(&tmp0 * &tmp1) } pub fn mul_by_pow_2(&self, k: u32) -> ExtendedPoint { From 28f10bc183b50ff41bc42fb6877921a65f55a47f Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:21:52 -0700 Subject: [PATCH 09/27] Change Lanes::ALL to Lanes::ABCD for consistency --- src/backend/avx2/edwards.rs | 2 +- src/backend/avx2/field.rs | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 563d3d6..1c7c616 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -215,7 +215,7 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { tmp.swap_CD(); // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - tmp.diff_sum(Lanes::ALL); + tmp.diff_sum(Lanes::ABCD); // set t0 = (S12 S15 S15 S12) let t0 = tmp.shuffle(Shuffle::ADDA); diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 13a6156..5caabed 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -36,7 +36,7 @@ pub enum Lanes { CD, AD, BC, - ALL, + ABCD, } #[inline(always)] @@ -67,7 +67,7 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { _mm256_blend_epi32(x.into_bits(), y.into_bits(), (B_LANES | C_LANES) as i32) .into_bits() } - Lanes::ALL => _mm256_blend_epi32( + Lanes::ABCD => _mm256_blend_epi32( x.into_bits(), y.into_bits(), (A_LANES | B_LANES | C_LANES | D_LANES) as i32, @@ -719,7 +719,7 @@ mod test { let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); - vec.diff_sum(Lanes::ALL); + vec.diff_sum(Lanes::ABCD); let result = vec.split(); From cc8728b2a6e604270012db81334aaeb349209188 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:22:28 -0700 Subject: [PATCH 10/27] Add comment about rustc-constant-info to blend function --- src/backend/avx2/field.rs | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 5caabed..d25d997 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -44,6 +44,26 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { unsafe { use core::arch::x86_64::_mm256_blend_epi32; + // This would be much cleaner if we could factor out the match + // statement on the control. Unfortunately, rustc forgets + // constant-info very quickly, so we can't even write + // ``` + // match control { + // Lanes::C => { + // let imm = C_LANES as i32; + // _mm256_blend_epi32(..., imm) + // ``` + // let alone + // ``` + // let imm = match control { + // Lanes::C => C_LANES as i32, + // } + // _mm256_blend_epi32(..., imm) + // ``` + // even though both of these would be constant-folded by LLVM + // at a lower level (as happens in the shuffle implementation, + // which does not require a shuffle immediate but *is* lowered + // to immediate shuffles anyways). match control { Lanes::C => { _mm256_blend_epi32(x.into_bits(), y.into_bits(), C_LANES as i32).into_bits() From 25d9f3f6cab38a5df86fcfb059a75249f371aa15 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 13 Jun 2018 16:41:57 -0700 Subject: [PATCH 11/27] Eliminate vector constants from `edwards` module --- src/backend/avx2/constants.rs | 18 ++++++++++++++++++ src/backend/avx2/edwards.rs | 21 +++------------------ src/backend/avx2/field.rs | 2 -- 3 files changed, 21 insertions(+), 20 deletions(-) diff --git a/src/backend/avx2/constants.rs b/src/backend/avx2/constants.rs index e99e312..f08daa7 100644 --- a/src/backend/avx2/constants.rs +++ b/src/backend/avx2/constants.rs @@ -16,6 +16,24 @@ use backend::avx2::edwards::{CachedPoint, ExtendedPoint}; use backend::avx2::field::FieldElement32x4; use scalar_mul::window::NafLookupTable8; +/// The identity element as an `ExtendedPoint`. +pub(crate) static EXTENDEDPOINT_IDENTITY: ExtendedPoint = ExtendedPoint(FieldElement32x4([ + u32x8::new(0, 1, 0, 0, 1, 0, 0, 0), + u32x8::splat(0), + u32x8::splat(0), + u32x8::splat(0), + u32x8::splat(0), +])); + +/// The identity element as a `CachedPoint`. +pub(crate) static CACHEDPOINT_IDENTITY: CachedPoint = CachedPoint(FieldElement32x4([ + u32x8::new(121647, 121666, 0, 0, 243332, 67108845, 0, 33554431), + u32x8::new(67108864, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), +])); + /// The low limbs of (2p, 2p, 2p, 2p), so that /// ```no_run /// (2p, 2p, 2p, 2p) = [P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI] diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 1c7c616..a04c785 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -16,8 +16,6 @@ use core::convert::From; use core::ops::{Add, Neg, Sub}; -use core::simd::{u32x8, IntoBits}; - use subtle::Choice; use subtle::ConditionallyAssignable; @@ -27,8 +25,7 @@ use scalar_mul::window::{LookupTable, NafLookupTable5, NafLookupTable8}; use traits::Identity; use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle}; - -use backend::avx2; +use backend::avx2::constants; /// A point on Curve25519, represented in an AVX2-friendly format. #[derive(Copy, Clone, Debug)] @@ -66,13 +63,7 @@ impl Default for ExtendedPoint { impl Identity for ExtendedPoint { fn identity() -> ExtendedPoint { - ExtendedPoint(FieldElement32x4([ - u32x8::new(0, 1, 0, 0, 1, 0, 0, 0), - u32x8::splat(0), - u32x8::splat(0), - u32x8::splat(0), - u32x8::splat(0), - ])) + constants::EXTENDEDPOINT_IDENTITY } } @@ -172,13 +163,7 @@ impl Default for CachedPoint { impl Identity for CachedPoint { fn identity() -> CachedPoint { - CachedPoint(FieldElement32x4([ - u32x8::new(121647, 121666, 0, 0, 243332, 67108845, 0, 33554431), - u32x8::new(67108864, 0, 33554431, 0, 0, 67108863, 0, 33554431), - u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), - u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), - u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), - ])) + constants::CACHEDPOINT_IDENTITY } } diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index d25d997..082287d 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -526,8 +526,6 @@ impl FieldElement32x4 { /// /// Limbs must be bounded by bit-excess \\( b < 2.0 \\). pub fn square_and_negate_D(&self) -> FieldElement32x4 { - let neg_mask = D_LANES64; - #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; From 02296fafb5ff50749c274ceb04a15e9fb9439711 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 13 Jun 2018 16:42:22 -0700 Subject: [PATCH 12/27] Delete unused constant --- src/backend/avx2/constants.rs | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/backend/avx2/constants.rs b/src/backend/avx2/constants.rs index f08daa7..304a1be 100644 --- a/src/backend/avx2/constants.rs +++ b/src/backend/avx2/constants.rs @@ -94,14 +94,6 @@ pub(crate) static P_TIMES_16_HI: u32x8 = u32x8::new( 33554431 << 4, ); -pub(crate) static P_TIMES_2_MASKED: FieldElement32x4 = FieldElement32x4([ - u32x8::new(0, 134217690, 0, 67108862, 134217690, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), - u32x8::new(0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), -]); - /// Odd multiples of the Ed25519 basepoint: pub(crate) static BASEPOINT_ODD_LOOKUP_TABLE: NafLookupTable8 = NafLookupTable8([ CachedPoint(FieldElement32x4([ From eea3eadf5b5b10249c77cdc0e351f5803bd6b68b Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 11:13:46 -0700 Subject: [PATCH 13/27] Replace special-case `swap_{AB,CD}` methods with general shuffles --- src/backend/avx2/edwards.rs | 7 +++---- src/backend/avx2/field.rs | 30 ++++-------------------------- 2 files changed, 7 insertions(+), 30 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index a04c785..df446b7 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -177,9 +177,8 @@ impl<'a> Neg for &'a CachedPoint { type Output = CachedPoint; fn neg(self) -> CachedPoint { - let mut coords = self.0; - coords.swap_AB(); - CachedPoint(coords.blend(coords.negate_lazy(), Lanes::D)) + let swapped = self.0.shuffle(Shuffle::BACD); + CachedPoint(swapped.blend(swapped.negate_lazy(), Lanes::D)) } } @@ -197,7 +196,7 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { tmp = &tmp * &other.0; // tmp = (S8 S9 S11 S10) - tmp.swap_CD(); + tmp = tmp.shuffle(Shuffle::ABDC); // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) tmp.diff_sum(Lanes::ABCD); diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 082287d..ea6e6cb 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -106,6 +106,8 @@ pub enum Shuffle { CBCB, ABAB, BADC, + BACD, + ABDC, } /// A vector of four `FieldElements`, implemented using AVX2. @@ -174,6 +176,8 @@ impl FieldElement32x4 { Shuffle::CBCB => u32x8::new(4, 1, 6, 3, 4, 1, 6, 3), Shuffle::ABAB => u32x8::new(0, 1, 2, 3, 0, 1, 2, 3), Shuffle::BADC => u32x8::new(1, 0, 3, 2, 5, 4, 7, 6), + Shuffle::BACD => u32x8::new(1, 0, 3, 2, 4, 5, 6, 7), + Shuffle::ABDC => u32x8::new(0, 1, 2, 3, 5, 4, 7, 6), }; // Note that this gets turned into a generic LLVM // shuffle-by-constants, which can be lowered to a simpler @@ -241,32 +245,6 @@ impl FieldElement32x4 { ]) } - /// Given `self = (A,B,C,D)`, set `self = (B,A,C,D)` - pub fn swap_AB(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - use core::arch::x86_64::_mm256_shuffle_epi32; - for i in 0..5 { - let swapped = _mm256_shuffle_epi32(self.0[i].into_bits(), 0b10_11_00_01); - self.0[i] = - _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b00001111).into_bits(); - } - } - } - - /// Given `self = (A,B,C,D)`, set `self = (A,B,D,C)` - pub fn swap_CD(&mut self) { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - use core::arch::x86_64::_mm256_shuffle_epi32; - for i in 0..5 { - let swapped = _mm256_shuffle_epi32(self.0[i].into_bits(), 0b10_11_00_01); - self.0[i] = - _mm256_blend_epi32(self.0[i].into_bits(), swapped, 0b11110000).into_bits(); - } - } - } - /// Given `self = (A,B,C,D)`, set `self = (B - A, B + A, D - C, D + C)` according to `mask`. /// /// This is `#[inline(always)]` because the `mask` parameter should be an immediate. From c46ec9638ceadd7c0f8493bf39b9e33c8e95917e Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 11:14:12 -0700 Subject: [PATCH 14/27] Add documentation --- src/backend/avx2/field.rs | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index ea6e6cb..1f64775 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -28,7 +28,13 @@ use core::simd::{i32x8, u32x8, u64x4, IntoBits}; use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO}; use backend::u64::field::FieldElement64; -#[derive(Copy, Clone)] +/// The `Lanes` enum represents a subset of the lanes `A,B,C,D` of a +/// `FieldElement32x4`. +/// +/// It's used to specify blend operations without +/// having to know details about the data layout of the +/// `FieldElement32x4`. +#[derive(Copy, Clone, Debug)] pub enum Lanes { C, D, @@ -96,7 +102,8 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { } } -#[derive(Copy, Clone)] +/// The `Shuffle` enum represents a shuffle of a `FieldElement32x4`. +#[derive(Copy, Clone, Debug)] pub enum Shuffle { AAAA, BBBB, From fe51adad319f2c7dbea15fbbcf872d2b32ee01c4 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:32:13 -0700 Subject: [PATCH 15/27] Don't expose u32x8 unpacking functions --- src/backend/avx2/field.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 1f64775..63a1c3f 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -470,7 +470,7 @@ impl FieldElement32x4 { } #[inline(always)] -pub fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { +fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { let a: u32x8; let b: u32x8; let zero = i32x8::new(0, 0, 0, 0, 0, 0, 0, 0); @@ -484,7 +484,7 @@ pub fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { } #[inline(always)] -pub fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { +fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { unsafe { use core::arch::x86_64::_mm256_blend_epi32; use core::arch::x86_64::_mm256_shuffle_epi32; From 64b1b481ba7467dea5ae3ef7fa4e0f6f8aa8363e Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 14:59:11 -0700 Subject: [PATCH 16/27] Rewrite diff_sum in terms of shuffle, blend, negate --- src/backend/avx2/edwards.rs | 6 +-- src/backend/avx2/field.rs | 80 ++++++++----------------------------- 2 files changed, 19 insertions(+), 67 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index df446b7..298317a 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -143,7 +143,7 @@ impl From for CachedPoint { let mut x = P.0; // x = (S2 S3 Z2 T2) - x.diff_sum(Lanes::AB); + x = x.blend(x.diff_sum(), Lanes::AB); // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) x.scale_by_curve_constants(); @@ -190,7 +190,7 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { let mut tmp = self.0; // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) - tmp.diff_sum(Lanes::AB); + tmp = tmp.blend(tmp.diff_sum(), Lanes::AB); // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) tmp = &tmp * &other.0; @@ -199,7 +199,7 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { tmp = tmp.shuffle(Shuffle::ABDC); // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - tmp.diff_sum(Lanes::ABCD); + tmp = tmp.diff_sum(); // set t0 = (S12 S15 S15 S12) let t0 = tmp.shuffle(Shuffle::ADDA); diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 63a1c3f..ea2a647 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -39,6 +39,7 @@ pub enum Lanes { C, D, AB, + AC, CD, AD, BC, @@ -85,6 +86,10 @@ fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32) .into_bits() } + Lanes::AC => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | C_LANES) as i32) + .into_bits() + } Lanes::CD => { _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32) .into_bits() @@ -242,6 +247,7 @@ impl FieldElement32x4 { /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. /// /// The output limbs are bounded by \\(2p\\). + #[inline] pub fn negate_lazy(&self) -> FieldElement32x4 { FieldElement32x4([ P_TIMES_2_LO - self.0[0], @@ -252,58 +258,15 @@ impl FieldElement32x4 { ]) } - /// Given `self = (A,B,C,D)`, set `self = (B - A, B + A, D - C, D + C)` according to `mask`. - /// - /// This is `#[inline(always)]` because the `mask` parameter should be an immediate. - #[inline(always)] - pub fn diff_sum(&mut self, control: Lanes) { - unsafe { - use core::arch::x86_64::{_mm256_blend_epi32, _mm256_shuffle_epi32}; - - let shuffle = |v: u32x8| -> u32x8 { - _mm256_shuffle_epi32(v.into_bits(), 0b10_11_00_01).into_bits() - }; - - let x01 = self.0[0]; - let x01_shuf = shuffle(x01); - let v1 = (x01_shuf + P_TIMES_2_LO) - x01; - let v2 = x01_shuf + x01; - let diffsum01 = - _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[0] = blend_lanes(x01, diffsum01, control); - - let x23 = self.0[1]; - let x23_shuf = shuffle(x23); - let v1 = (x23_shuf + P_TIMES_2_HI) - x23; - let v2 = x23_shuf + x23; - let diffsum23 = - _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[1] = blend_lanes(x23, diffsum23, control); - - let x45 = self.0[2]; - let x45_shuf = shuffle(x45); - let v1 = (x45_shuf + P_TIMES_2_HI) - x45; - let v2 = x45_shuf + x45; - let diffsum45 = - _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[2] = blend_lanes(x45, diffsum45, control); - - let x67 = self.0[3]; - let x67_shuf = shuffle(x67); - let v1 = (x67_shuf + P_TIMES_2_HI) - x67; - let v2 = x67_shuf + x67; - let diffsum67 = - _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[3] = blend_lanes(x67, diffsum67, control); - - let x89 = self.0[4]; - let x89_shuf = shuffle(x89); - let v1 = (x89_shuf + P_TIMES_2_HI) - x89; - let v2 = x89_shuf + x89; - let diffsum89 = - _mm256_blend_epi32(v1.into_bits(), v2.into_bits(), 0b10101010).into_bits(); - self.0[4] = blend_lanes(x89, diffsum89, control); - } + /// Given `self = (A,B,C,D)`, compute `(B - A, B + A, D - C, D + C)`. + #[inline] + pub fn diff_sum(&self) -> FieldElement32x4 { + // tmp1 = (B, A, D, C) + let tmp1 = self.shuffle(Shuffle::BADC); + // tmp2 = (-A, B, -C, D) + let tmp2 = self.blend(self.negate_lazy(), Lanes::AC); + // (B - A, B + A, D - C, D + C) + tmp1 + tmp2 } /// Let `self` \\(= (A, B, C, D) \\). @@ -721,8 +684,7 @@ mod test { let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]); let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); - let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); - vec.diff_sum(Lanes::ABCD); + let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3).diff_sum(); let result = vec.split(); @@ -730,16 +692,6 @@ mod test { assert_eq!(result[1], &x1 + &x0); assert_eq!(result[2], &x3 - &x2); assert_eq!(result[3], &x3 + &x2); - - let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); - vec.diff_sum(Lanes::AB); // leave C,D unchanged - - let result = vec.split(); - - assert_eq!(result[0], &x1 - &x0); - assert_eq!(result[1], &x1 + &x0); - assert_eq!(result[2], x2); - assert_eq!(result[3], x3); } #[test] From 4dc219910aa377c6a1b8da728e058491e60e78b9 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 15:13:53 -0700 Subject: [PATCH 17/27] Move blend_lanes into the blend function --- src/backend/avx2/field.rs | 146 +++++++++++++++++++------------------- 1 file changed, 73 insertions(+), 73 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index ea2a647..06a22fd 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -46,67 +46,6 @@ pub enum Lanes { ABCD, } -#[inline(always)] -fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - - // This would be much cleaner if we could factor out the match - // statement on the control. Unfortunately, rustc forgets - // constant-info very quickly, so we can't even write - // ``` - // match control { - // Lanes::C => { - // let imm = C_LANES as i32; - // _mm256_blend_epi32(..., imm) - // ``` - // let alone - // ``` - // let imm = match control { - // Lanes::C => C_LANES as i32, - // } - // _mm256_blend_epi32(..., imm) - // ``` - // even though both of these would be constant-folded by LLVM - // at a lower level (as happens in the shuffle implementation, - // which does not require a shuffle immediate but *is* lowered - // to immediate shuffles anyways). - match control { - Lanes::C => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), C_LANES as i32).into_bits() - } - Lanes::D => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), D_LANES as i32).into_bits() - } - Lanes::AD => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | D_LANES) as i32) - .into_bits() - } - Lanes::AB => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32) - .into_bits() - } - Lanes::AC => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | C_LANES) as i32) - .into_bits() - } - Lanes::CD => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32) - .into_bits() - } - Lanes::BC => { - _mm256_blend_epi32(x.into_bits(), y.into_bits(), (B_LANES | C_LANES) as i32) - .into_bits() - } - Lanes::ABCD => _mm256_blend_epi32( - x.into_bits(), - y.into_bits(), - (A_LANES | B_LANES | C_LANES | D_LANES) as i32, - ).into_bits(), - } - } -} - /// The `Shuffle` enum represents a shuffle of a `FieldElement32x4`. #[derive(Copy, Clone, Debug)] pub enum Shuffle { @@ -161,18 +100,7 @@ impl FieldElement32x4 { out } - #[inline(always)] - pub fn blend(&self, b: FieldElement32x4, control: Lanes) -> FieldElement32x4 { - FieldElement32x4([ - blend_lanes(self.0[0], b.0[0], control), - blend_lanes(self.0[1], b.0[1], control), - blend_lanes(self.0[2], b.0[2], control), - blend_lanes(self.0[3], b.0[3], control), - blend_lanes(self.0[4], b.0[4], control), - ]) - } - - #[inline(always)] + #[inline] pub fn shuffle(&self, control: Shuffle) -> FieldElement32x4 { #[inline(always)] fn shuffle_lanes(x: u32x8, control: Shuffle) -> u32x8 { @@ -207,6 +135,78 @@ impl FieldElement32x4 { ]) } + #[inline] + pub fn blend(&self, b: FieldElement32x4, control: Lanes) -> FieldElement32x4 { + #[inline(always)] + fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { + unsafe { + use core::arch::x86_64::_mm256_blend_epi32; + + // This would be much cleaner if we could factor out the match + // statement on the control. Unfortunately, rustc forgets + // constant-info very quickly, so we can't even write + // ``` + // match control { + // Lanes::C => { + // let imm = C_LANES as i32; + // _mm256_blend_epi32(..., imm) + // ``` + // let alone + // ``` + // let imm = match control { + // Lanes::C => C_LANES as i32, + // } + // _mm256_blend_epi32(..., imm) + // ``` + // even though both of these would be constant-folded by LLVM + // at a lower level (as happens in the shuffle implementation, + // which does not require a shuffle immediate but *is* lowered + // to immediate shuffles anyways). + match control { + Lanes::C => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), C_LANES as i32).into_bits() + } + Lanes::D => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), D_LANES as i32).into_bits() + } + Lanes::AD => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | D_LANES) as i32) + .into_bits() + } + Lanes::AB => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | B_LANES) as i32) + .into_bits() + } + Lanes::AC => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (A_LANES | C_LANES) as i32) + .into_bits() + } + Lanes::CD => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (C_LANES | D_LANES) as i32) + .into_bits() + } + Lanes::BC => { + _mm256_blend_epi32(x.into_bits(), y.into_bits(), (B_LANES | C_LANES) as i32) + .into_bits() + } + Lanes::ABCD => _mm256_blend_epi32( + x.into_bits(), + y.into_bits(), + (A_LANES | B_LANES | C_LANES | D_LANES) as i32, + ).into_bits(), + } + } + } + + FieldElement32x4([ + blend_lanes(self.0[0], b.0[0], control), + blend_lanes(self.0[1], b.0[1], control), + blend_lanes(self.0[2], b.0[2], control), + blend_lanes(self.0[3], b.0[3], control), + blend_lanes(self.0[4], b.0[4], control), + ]) + } + pub fn zero() -> FieldElement32x4 { FieldElement32x4([u32x8::splat(0); 5]) } From 97292fef91c789682a01e78d351f44c02afaab44 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 15:21:51 -0700 Subject: [PATCH 18/27] Move packing functions to top of the module --- src/backend/avx2/field.rs | 86 +++++++++++++++++++++++---------------- 1 file changed, 51 insertions(+), 35 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 06a22fd..45db01b 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -28,6 +28,57 @@ use core::simd::{i32x8, u32x8, u64x4, IntoBits}; use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO}; use backend::u64::field::FieldElement64; +/// Unpack +/// ``` +/// (a0, b0, a1, b1, c0, d0, c1, d1) +/// ``` +/// into +/// ``` +/// (a0, 0, b0, 0, c0, 0, d0, 0) +/// (a1, 0, b1, 0, c1, 0, d1, 0) +/// ``` +#[inline(always)] +fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { + let a: u32x8; + let b: u32x8; + let zero = i32x8::new(0, 0, 0, 0, 0, 0, 0, 0); + unsafe { + use core::arch::x86_64::_mm256_unpackhi_epi32; + use core::arch::x86_64::_mm256_unpacklo_epi32; + a = _mm256_unpacklo_epi32(src.into_bits(), zero.into_bits()).into_bits(); + b = _mm256_unpackhi_epi32(src.into_bits(), zero.into_bits()).into_bits(); + } + (a, b) +} + +/// Repack +/// ``` +/// (a0, 0, b0, 0, c0, 0, d0, 0) +/// (a1, 0, b1, 0, c1, 0, d1, 0) +/// ``` +/// into +/// ``` +/// (a0, b0, a1, b1, c0, d0, c1, d1) +/// ``` +#[inline(always)] +fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { + unsafe { + use core::arch::x86_64::_mm256_blend_epi32; + use core::arch::x86_64::_mm256_shuffle_epi32; + + // Input: x = (a0, 0, b0, 0, c0, 0, d0, 0) + // Input: y = (a1, 0, b1, 0, c1, 0, d1, 0) + + let x_shuffled = _mm256_shuffle_epi32(x.into_bits(), 0b11_01_10_00); + let y_shuffled = _mm256_shuffle_epi32(y.into_bits(), 0b10_00_11_01); + + // x' = (a0, b0, 0, 0, c0, d0, 0, 0) + // y' = ( 0, 0, a1, b1, 0, 0, c1, d1) + + return _mm256_blend_epi32(x_shuffled, y_shuffled, 0b11001100).into_bits(); + } +} + /// The `Lanes` enum represents a subset of the lanes `A,B,C,D` of a /// `FieldElement32x4`. /// @@ -430,42 +481,7 @@ impl FieldElement32x4 { repack_pair(z[8].into_bits(), z[9].into_bits()), ]) } -} -#[inline(always)] -fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { - let a: u32x8; - let b: u32x8; - let zero = i32x8::new(0, 0, 0, 0, 0, 0, 0, 0); - unsafe { - use core::arch::x86_64::_mm256_unpackhi_epi32; - use core::arch::x86_64::_mm256_unpacklo_epi32; - a = _mm256_unpacklo_epi32(src.into_bits(), zero.into_bits()).into_bits(); - b = _mm256_unpackhi_epi32(src.into_bits(), zero.into_bits()).into_bits(); - } - (a, b) -} - -#[inline(always)] -fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { - unsafe { - use core::arch::x86_64::_mm256_blend_epi32; - use core::arch::x86_64::_mm256_shuffle_epi32; - - // Input: x = (a0, 0, b0, 0, c0, 0, d0) - // Input: y = (a1, 0, b1, 0, c1, 0, d1) - - let x_shuffled = _mm256_shuffle_epi32(x.into_bits(), 0b11_01_10_00); - let y_shuffled = _mm256_shuffle_epi32(y.into_bits(), 0b10_00_11_01); - - // x' = (a0, b0, 0, 0, c0, d0, 0, 0) - // y' = ( 0, 0, a1, b1, 0, 0, c1, d1) - - return _mm256_blend_epi32(x_shuffled, y_shuffled, 0b11001100).into_bits(); - } -} - -impl FieldElement32x4 { /// Square this field element, then conditionally negate according /// to `neg_mask`. This parameter is hardcoded as `neg_mask = /// D_LANES64` to negate the \\( D \\) value. From 6ef9e9dcfbdef851402e27918fd58965fa7dd10d Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 15:22:16 -0700 Subject: [PATCH 19/27] Make publicity a little more consistent --- src/backend/avx2/field.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 45db01b..a3906b4 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -114,7 +114,7 @@ pub enum Shuffle { /// A vector of four `FieldElements`, implemented using AVX2. #[derive(Clone, Copy, Debug)] -pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]); +pub struct FieldElement32x4(pub(crate) [u32x8; 5]); use subtle::Choice; use subtle::ConditionallyAssignable; @@ -130,7 +130,7 @@ impl ConditionallyAssignable for FieldElement32x4 { } impl FieldElement32x4 { - pub(crate) fn split(&self) -> [FieldElement64; 4] { + pub fn split(&self) -> [FieldElement64; 4] { let mut out = [FieldElement64::zero(); 4]; for i in 0..5 { let a_2i = self.0[i].extract(0) as u64; // @@ -418,7 +418,7 @@ impl FieldElement32x4 { v[0] = v[0] + c9_19; } - pub fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { + fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { // These aren't const because splat isn't a const fn let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1); let LOW_26_BITS: u64x4 = u64x4::splat((1 << 26) - 1); From bab1ebbeb8f9d5208cc2343383b4d7ee3eef0648 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 15:39:32 -0700 Subject: [PATCH 20/27] Replace scale_by_curve_constants by a Mul<(u32,u32,u32,u32)> impl --- src/backend/avx2/edwards.rs | 2 +- src/backend/avx2/field.rs | 70 +++++++++++++++++-------------------- 2 files changed, 34 insertions(+), 38 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 298317a..d688fca 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -146,7 +146,7 @@ impl From for CachedPoint { x = x.blend(x.diff_sum(), Lanes::AB); // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) - x.scale_by_curve_constants(); + x = x * (121666, 121666, 2*121666, 2*121665); // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) x = x.blend(-x, Lanes::D); diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index a3906b4..0c3aa75 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -320,42 +320,6 @@ impl FieldElement32x4 { tmp1 + tmp2 } - /// Let `self` \\(= (A, B, C, D) \\). - /// - /// Compute - /// $$( 121666A, 121666B, 2\cdot 121666C, 2\cdot 121665 D).$$ - pub fn scale_by_curve_constants(&mut self) { - let mut b = [u64x4::splat(0); 10]; - - let consts = u32x8::new(121666, 0, 121666, 0, 2 * 121666, 0, 2 * 121665, 0); - - unsafe { - use core::arch::x86_64::_mm256_mul_epu32; - - let (b0, b1) = unpack_pair(self.0[0]); - b[0] = _mm256_mul_epu32(b0.into_bits(), consts.into_bits()).into_bits(); - b[1] = _mm256_mul_epu32(b1.into_bits(), consts.into_bits()).into_bits(); - - let (b2, b3) = unpack_pair(self.0[1]); - b[2] = _mm256_mul_epu32(b2.into_bits(), consts.into_bits()).into_bits(); - b[3] = _mm256_mul_epu32(b3.into_bits(), consts.into_bits()).into_bits(); - - let (b4, b5) = unpack_pair(self.0[2]); - b[4] = _mm256_mul_epu32(b4.into_bits(), consts.into_bits()).into_bits(); - b[5] = _mm256_mul_epu32(b5.into_bits(), consts.into_bits()).into_bits(); - - let (b6, b7) = unpack_pair(self.0[3]); - b[6] = _mm256_mul_epu32(b6.into_bits(), consts.into_bits()).into_bits(); - b[7] = _mm256_mul_epu32(b7.into_bits(), consts.into_bits()).into_bits(); - - let (b8, b9) = unpack_pair(self.0[4]); - b[8] = _mm256_mul_epu32(b8.into_bits(), consts.into_bits()).into_bits(); - b[9] = _mm256_mul_epu32(b9.into_bits(), consts.into_bits()).into_bits(); - } - - *self = FieldElement32x4::reduce64(b); - } - pub fn reduce32(&mut self) { let shifts = i32x8::new(26, 26, 25, 25, 26, 26, 25, 25); let masks = u32x8::new( @@ -616,6 +580,37 @@ impl Add for FieldElement32x4 { } } +impl Mul<(u32, u32, u32, u32)> for FieldElement32x4 { + type Output = FieldElement32x4; + #[inline] + fn mul(self, scalars: (u32, u32, u32, u32)) -> FieldElement32x4 { + unsafe { + use core::arch::x86_64::_mm256_mul_epu32; + + let consts = u32x8::new(scalars.0, 0, scalars.1, 0, scalars.2, 0, scalars.3, 0); + + let (b0, b1) = unpack_pair(self.0[0]); + let (b2, b3) = unpack_pair(self.0[1]); + let (b4, b5) = unpack_pair(self.0[2]); + let (b6, b7) = unpack_pair(self.0[3]); + let (b8, b9) = unpack_pair(self.0[4]); + + FieldElement32x4::reduce64([ + _mm256_mul_epu32(b0.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b1.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b2.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b3.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b4.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b5.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b6.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b7.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b8.into_bits(), consts.into_bits()).into_bits(), + _mm256_mul_epu32(b9.into_bits(), consts.into_bits()).into_bits(), + ]) + } + } +} + impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { type Output = FieldElement32x4; fn mul(self, _rhs: &'b FieldElement32x4) -> FieldElement32x4 { @@ -684,7 +679,8 @@ mod test { #[test] fn scale_by_curve_constants() { let mut x = FieldElement32x4::splat(&FieldElement64::one()); - x.scale_by_curve_constants(); + + x = x * (121666, 121666, 2*121666, 2*121665); let xs = x.split(); assert_eq!(xs[0], FieldElement64([121666, 0, 0, 0, 0])); From 9f5bd8c4c0794fee6d633c89d56b257ab86fc7a5 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 15:54:36 -0700 Subject: [PATCH 21/27] Rename reduce32 to reduce and have it return its result. This means that all FieldElement32x4 operations return values, vs mutating interior state. --- src/backend/avx2/field.rs | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 0c3aa75..297ba34 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -287,9 +287,10 @@ impl FieldElement32x4 { buf[i] = u32x8::new(a_2i, b_2i, a_2i_1, b_2i_1, c_2i, d_2i, c_2i_1, d_2i_1); } - let mut out = FieldElement32x4(buf); - out.reduce32(); - return out; + // We don't know that the original `FieldElement64`s were + // fully reduced, so the odd limbs may exceed 2^25. + // Reduce them to be sure. + FieldElement32x4(buf).reduce() } /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), without @@ -320,7 +321,9 @@ impl FieldElement32x4 { tmp1 + tmp2 } - pub fn reduce32(&mut self) { + /// Compute the reduced form of the field elements. + #[inline] + pub fn reduce(&self) -> FieldElement32x4 { let shifts = i32x8::new(26, 26, 25, 25, 26, 26, 25, 25); let masks = u32x8::new( (1 << 26) - 1, @@ -333,6 +336,7 @@ impl FieldElement32x4 { (1 << 25) - 1, ); + // Compute the carryout of v. let carry = |v: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_srlv_epi32; @@ -340,6 +344,7 @@ impl FieldElement32x4 { } }; + // Swap adjacent 32-bit lanes. let swap_lanes = |v: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_shuffle_epi32; @@ -354,7 +359,7 @@ impl FieldElement32x4 { } }; - let v = &mut self.0; + let mut v = self.0; let c10 = swap_lanes(carry(v[0])); v[0] = (v[0] & masks) + combine(u32x8::splat(0), c10); @@ -380,6 +385,8 @@ impl FieldElement32x4 { } v[0] = v[0] + c9_19; + + FieldElement32x4(v) } fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { @@ -553,16 +560,13 @@ impl Neg for FieldElement32x4 { /// The output limbs are freshly reduced. #[inline] fn neg(self) -> FieldElement32x4 { - let mut neg = FieldElement32x4([ + FieldElement32x4([ P_TIMES_16_LO - self.0[0], P_TIMES_16_HI - self.0[1], P_TIMES_16_HI - self.0[2], P_TIMES_16_HI - self.0[3], P_TIMES_16_HI - self.0[4], - ]); - neg.reduce32(); - - neg + ]).reduce() } } From 7198719419b39479850548d0b5e98cbea5645120 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 14 Jun 2018 17:15:16 -0700 Subject: [PATCH 22/27] Document bounds on FieldElement32x4 functions --- src/backend/avx2/edwards.rs | 57 +++++---- src/backend/avx2/field.rs | 235 ++++++++++++++++++++++++++++-------- 2 files changed, 221 insertions(+), 71 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index d688fca..7c5ff70 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -86,7 +86,7 @@ impl ExtendedPoint { // Set tmp1 = tmp0^2, negating the D values tmp1 = tmp0.square_and_negate_D(); - // Now tmp1 = (S1 S2 S3 -S4) + // Now tmp1 = (S1 S2 S3 -S4) with b < 0.007 // See discussion of bounds in the module-level documentation. // We want to compute @@ -105,23 +105,27 @@ impl ExtendedPoint { let S_1 = tmp1.shuffle(Shuffle::AAAA); let S_2 = tmp1.shuffle(Shuffle::BBBB); - // tmp0 = (0, 0, 2S_3, 0) tmp0 = zero.blend(tmp1 + tmp1, Lanes::C); - // tmp0 = (0, 0, 2S_3, -S_4) + // tmp0 = (0, 0, 2S_3, 0) tmp0 = tmp0.blend(tmp1, Lanes::D); - // tmp0 = (S_1, S_1, S_1 + 2S_3, S_1 - S_4) + // tmp0 = (0, 0, 2S_3, -S_4) tmp0 = tmp0 + S_1; - // tmp0 = (S_1 + S_2, S_1, S_1 + 2S_3, S_1 + S_2 - S_4) + // tmp0 = ( S_1, S_1, S_1 + 2S_3, S_1 - S_4) tmp0 = tmp0 + zero.blend(S_2, Lanes::AD); - // tmp0 = (S_1 + S_2, S_1 - S_2, S_1 - S_2 + 2S_3, S_1 + S_2 - S_4) + // tmp0 = (S_1 + S_2, S_1, S_1 + 2S_3, S_1 + S_2 - S_4) tmp0 = tmp0 + zero.blend(S_2.negate_lazy(), Lanes::BC); + // tmp0 = (S_1 + S_2, S_1 - S_2, S_1 - S_2 + 2S_3, S_1 + S_2 - S_4) + // b < ( 1.01, 1.6, 2.33, 1.6) // Now tmp0 = (S_5, S_6, S_8, S_9) - // Set tmp1 = (S_9, S_6, S_6, S_9) + // Set tmp1 = ( S_9, S_6, S_6, S_9) + // b < ( 1.6, 1.6, 1.6, 1.6) tmp1 = tmp0.shuffle(Shuffle::DBBD); - // Set tmp1 = (S_8, S_5, S_8, S_5) + // Set tmp1 = ( S_8, S_5, S_8, S_5) + // b < (2.33, 1.01, 2.33, 1.01) tmp0 = tmp0.shuffle(Shuffle::CACA); + // Bounds on (tmp0, tmp1) are (2.33, 1.6) < (2.5, 1.75). ExtendedPoint(&tmp0 * &tmp1) } @@ -142,15 +146,16 @@ impl From for CachedPoint { fn from(P: ExtendedPoint) -> CachedPoint { let mut x = P.0; - // x = (S2 S3 Z2 T2) x = x.blend(x.diff_sum(), Lanes::AB); + // x = (X1 - Y1, X2 + Y2, Z2, T2) = (S2 S3 Z2 T2) - // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) x = x * (121666, 121666, 2*121666, 2*121665); + // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) - // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) x = x.blend(-x, Lanes::D); + // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) + // The coefficients of the output are bounded with b < 0.007. CachedPoint(x) } } @@ -175,7 +180,12 @@ impl ConditionallyAssignable for CachedPoint { impl<'a> Neg for &'a CachedPoint { type Output = CachedPoint; - + /// Lazily negate the point. + /// + /// # Warning + /// + /// Because this method does not perform a reduction, it is not + /// safe to repeatedly negate a point. fn neg(self) -> CachedPoint { let swapped = self.0.shuffle(Shuffle::BACD); CachedPoint(swapped.blend(swapped.negate_lazy(), Lanes::D)) @@ -185,28 +195,30 @@ impl<'a> Neg for &'a CachedPoint { impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { type Output = ExtendedPoint; - /// Uses a slight tweak of the parallel unified formulas of HWCD'08 + /// Add an `ExtendedPoint` and a `CachedPoint`. fn add(self, other: &'b CachedPoint) -> ExtendedPoint { let mut tmp = self.0; - // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) tmp = tmp.blend(tmp.diff_sum(), Lanes::AB); + // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) with b < 1.6 - // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) + // (tmp, other) bounded with b < (1.6, 1.0) < (2.5, 1.75). tmp = &tmp * &other.0; + // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) - // tmp = (S8 S9 S11 S10) tmp = tmp.shuffle(Shuffle::ABDC); + // tmp = (S8 S9 S11 S10) - // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) tmp = tmp.diff_sum(); + // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - // set t0 = (S12 S15 S15 S12) let t0 = tmp.shuffle(Shuffle::ADDA); - // set t1 = (S14 S13 S14 S13) + // t0 = (S12 S15 S15 S12) let t1 = tmp.shuffle(Shuffle::CBCB); + // t1 = (S14 S13 S14 S13) - // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + // All coefficients of t0, t1 are bounded with b < 1.6. + // Return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) ExtendedPoint(&t0 * &t1) } } @@ -216,8 +228,9 @@ impl<'a, 'b> Sub<&'b CachedPoint> for &'a ExtendedPoint { /// Implement subtraction by negating the point and adding. /// - /// Empirically, this seems about the same cost as a custom subtraction impl (maybe because the - /// benefit is cancelled by increased code size?) + /// Empirically, this seems about the same cost as a custom + /// subtraction impl (maybe because the benefit is cancelled by + /// increased code size?) fn sub(self, other: &'b CachedPoint) -> ExtendedPoint { self + &(-other) } diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 297ba34..1484d75 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -112,7 +112,7 @@ pub enum Shuffle { ABDC, } -/// A vector of four `FieldElements`, implemented using AVX2. +/// A vector of four field elements, in an AVX2-friendly format. #[derive(Clone, Copy, Debug)] pub struct FieldElement32x4(pub(crate) [u32x8; 5]); @@ -258,14 +258,21 @@ impl FieldElement32x4 { ]) } + /// Construct a vector of zeros. pub fn zero() -> FieldElement32x4 { FieldElement32x4([u32x8::splat(0); 5]) } + /// Convenience wrapper around `new(x,x,x,x)`. pub fn splat(x: &FieldElement64) -> FieldElement32x4 { FieldElement32x4::new(x, x, x, x) } + /// Create a `FieldElement32x4` from four `FieldElement64`s. + /// + /// # Postconditions + /// + /// The resulting `FieldElement32x4` is bounded with \\( b < 0.0002 \\). pub fn new( x0: &FieldElement64, x1: &FieldElement64, @@ -296,11 +303,18 @@ impl FieldElement32x4 { /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), without /// performing a reduction. /// - /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. + /// # Preconditions /// - /// The output limbs are bounded by \\(2p\\). + /// The coefficients of `self` must be bounded with \\( b < 0.999 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 1 \\). #[inline] pub fn negate_lazy(&self) -> FieldElement32x4 { + // The limbs of self are bounded with b < 0.999, while the + // smallest limb of 2*p is 67108845 > 2^{26+0.9999}, so + // underflows are not possible. FieldElement32x4([ P_TIMES_2_LO - self.0[0], P_TIMES_2_HI - self.0[1], @@ -311,17 +325,29 @@ impl FieldElement32x4 { } /// Given `self = (A,B,C,D)`, compute `(B - A, B + A, D - C, D + C)`. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 0.01 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 1.6 \\). #[inline] pub fn diff_sum(&self) -> FieldElement32x4 { // tmp1 = (B, A, D, C) let tmp1 = self.shuffle(Shuffle::BADC); // tmp2 = (-A, B, -C, D) let tmp2 = self.blend(self.negate_lazy(), Lanes::AC); - // (B - A, B + A, D - C, D + C) + // (B - A, B + A, D - C, D + C) bounded with b < 1.6 tmp1 + tmp2 } - /// Compute the reduced form of the field elements. + /// Reduce this vector of field elements \\(\mathrm{mod} p\\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.0002 \\). #[inline] pub fn reduce(&self) -> FieldElement32x4 { let shifts = i32x8::new(26, 26, 25, 25, 26, 26, 25, 25); @@ -336,22 +362,36 @@ impl FieldElement32x4 { (1 << 25) - 1, ); - // Compute the carryout of v. - let carry = |v: u32x8| -> u32x8 { + // Let c(x) denote the carryout of the coefficient x. + // + // Given ( x0, y0, x1, y1, z0, w0, z1, w1), + // compute (c(x1), c(y1), c(x0), c(y0), c(z1), c(w1), c(z0), c(w0)). + // + // The carryouts are bounded by 2^(32 - 25) = 2^7. + let rotated_carryout = |v: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_srlv_epi32; - _mm256_srlv_epi32(v.into_bits(), shifts.into_bits()).into_bits() - } - }; - - // Swap adjacent 32-bit lanes. - let swap_lanes = |v: u32x8| -> u32x8 { - unsafe { use core::arch::x86_64::_mm256_shuffle_epi32; - _mm256_shuffle_epi32(v.into_bits(), 0b01_00_11_10).into_bits() + + let c = _mm256_srlv_epi32(v.into_bits(), shifts.into_bits()); + _mm256_shuffle_epi32(c, 0b01_00_11_10).into_bits() } }; + // Combine (lo, lo, lo, lo, lo, lo, lo, lo) + // with (hi, hi, hi, hi, hi, hi, hi, hi) + // to (lo, lo, hi, hi, lo, lo, hi, hi) + // + // This allows combining carryouts, e.g., + // + // lo (c(x1), c(y1), c(x0), c(y0), c(z1), c(w1), c(z0), c(w0)) + // hi (c(x3), c(y3), c(x2), c(y2), c(z3), c(w3), c(z2), c(w2)) + // -> (c(x1), c(y1), c(x2), c(y2), c(z1), c(w1), c(z2), c(w2)) + // + // which is exactly the vector of carryins for + // + // ( x2, y2, x3, y3, z2, w2, z3, w3). + // let combine = |v_lo: u32x8, v_hi: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_blend_epi32; @@ -361,34 +401,62 @@ impl FieldElement32x4 { let mut v = self.0; - let c10 = swap_lanes(carry(v[0])); + let c10 = rotated_carryout(v[0]); v[0] = (v[0] & masks) + combine(u32x8::splat(0), c10); - let c32 = swap_lanes(carry(v[1])); + + let c32 = rotated_carryout(v[1]); v[1] = (v[1] & masks) + combine(c10, c32); - let c54 = swap_lanes(carry(v[2])); + + let c54 = rotated_carryout(v[2]); v[2] = (v[2] & masks) + combine(c32, c54); - let c76 = swap_lanes(carry(v[3])); + + let c76 = rotated_carryout(v[3]); v[3] = (v[3] & masks) + combine(c54, c76); - let c98 = swap_lanes(carry(v[4])); + + let c98 = rotated_carryout(v[4]); v[4] = (v[4] & masks) + combine(c76, c98); - // Still need to account for c9 - // c98 = (c9, c9, c8, c8, c9, c9, c8, c8) - // - let c9_19: u32x8; - unsafe { + let c9_19: u32x8 = unsafe { use core::arch::x86_64::_mm256_mul_epu32; use core::arch::x86_64::_mm256_shuffle_epi32; - let c9_spread = _mm256_shuffle_epi32(c98.into_bits(), 0b11_01_10_00); - let c9_19_spread = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into_bits()); - c9_19 = _mm256_shuffle_epi32(c9_19_spread, 0b11_01_10_00).into_bits(); - } + // Need to rearrange c98, since vpmuludq uses the low + // 32-bits of each 64-bit lane to compute the product: + // + // c98 = (c(x9), c(y9), c(x8), c(y8), c(z9), c(w9), c(z8), c(w8)); + // c9_spread = (c(x9), c(x8), c(y9), c(y8), c(z9), c(z8), c(w9), c(w8)). + let c9_spread = _mm256_shuffle_epi32(c98.into_bits(), 0b11_01_10_00); + + // Since the carryouts are bounded by 2^7, their products with 19 + // are bounded by 2^11.25. This means that + // + // c9_19_spread = (19*c(x9), 0, 19*c(y9), 0, 19*c(z9), 0, 19*c(w9), 0). + let c9_19_spread = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into_bits()); + + // Unshuffle: + // c9_19 = (19*c(x9), 19*c(y9), 0, 0, 19*c(z9), 19*c(w9), 0, 0). + _mm256_shuffle_epi32(c9_19_spread, 0b11_01_10_00).into_bits() + }; + + // Add the final carryin. v[0] = v[0] + c9_19; + // Each output coefficient has exactly one carryin, which is + // bounded by 2^11.25, so they are bounded as + // + // c_even < 2^26 + 2^11.25 < 26.00006 < 2^{26+b} + // c_odd < 2^25 + 2^11.25 < 25.0001 < 2^{25+b} + // + // where b = 0.0002. FieldElement32x4(v) } + /// Given an array of wide coefficients, reduce them to a `FieldElement32x4`. + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + #[inline] fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { // These aren't const because splat isn't a const fn let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1); @@ -442,8 +510,13 @@ impl FieldElement32x4 { z[1] = z[1] + c1; // z1 < 2^25 + 2^17.25 < 2^25.0067 carry(&mut z, 0); // z0 < 2^26, z1 < 2^25.0067 + 2^4.33 = 2^25.007 - // Now repack the [u64x4; 10] into a FieldElement32x4 - + // The output coefficients are bounded with + // + // b = 0.007 for z[1] + // b = 0.0004 for z[5] + // b = 0 for other z[i]. + // + // So the packed result is bounded with b = 0.007. FieldElement32x4([ repack_pair(z[0].into_bits(), z[1].into_bits()), repack_pair(z[2].into_bits(), z[3].into_bits()), @@ -453,13 +526,15 @@ impl FieldElement32x4 { ]) } - /// Square this field element, then conditionally negate according - /// to `neg_mask`. This parameter is hardcoded as `neg_mask = - /// D_LANES64` to negate the \\( D \\) value. + /// Square this field element, and negate the result's \\(D\\) value. /// - /// # Precondition + /// # Preconditions /// - /// Limbs must be bounded by bit-excess \\( b < 2.0 \\). + /// The coefficients of `self` must be bounded with \\( b < 1.5 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). pub fn square_and_negate_D(&self) -> FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { @@ -552,12 +627,18 @@ impl FieldElement32x4 { impl Neg for FieldElement32x4 { type Output = FieldElement32x4; - /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), and - /// perform a reduction. + /// Negate this field element, performing a reduction. /// - /// The input limbs can be any size. + /// If the coefficients are known to be small, use `negate_lazy` + /// to avoid performing a reduction. /// - /// The output limbs are freshly reduced. + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 4.0 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.0002 \\). #[inline] fn neg(self) -> FieldElement32x4 { FieldElement32x4([ @@ -572,6 +653,7 @@ impl Neg for FieldElement32x4 { impl Add for FieldElement32x4 { type Output = FieldElement32x4; + /// Add two `FieldElement32x4`s, without performing a reduction. #[inline] fn add(self, rhs: FieldElement32x4) -> FieldElement32x4 { FieldElement32x4([ @@ -586,6 +668,11 @@ impl Add for FieldElement32x4 { impl Mul<(u32, u32, u32, u32)> for FieldElement32x4 { type Output = FieldElement32x4; + /// Perform a multiplication by a vector of small constants. + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). #[inline] fn mul(self, scalars: (u32, u32, u32, u32)) -> FieldElement32x4 { unsafe { @@ -617,7 +704,19 @@ impl Mul<(u32, u32, u32, u32)> for FieldElement32x4 { impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { type Output = FieldElement32x4; - fn mul(self, _rhs: &'b FieldElement32x4) -> FieldElement32x4 { + /// Multiply `self` by `rhs`. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 2.5 \\). + /// + /// The coefficients of `rhs` must be bounded with \\( b < 1.75 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + /// + fn mul(self, rhs: &'b FieldElement32x4) -> FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; @@ -636,11 +735,11 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let (x6, x7) = unpack_pair(self.0[3]); let (x8, x9) = unpack_pair(self.0[4]); - let (y0, y1) = unpack_pair(_rhs.0[0]); - let (y2, y3) = unpack_pair(_rhs.0[1]); - let (y4, y5) = unpack_pair(_rhs.0[2]); - let (y6, y7) = unpack_pair(_rhs.0[3]); - let (y8, y9) = unpack_pair(_rhs.0[4]); + let (y0, y1) = unpack_pair(rhs.0[0]); + let (y2, y3) = unpack_pair(rhs.0[1]); + let (y4, y5) = unpack_pair(rhs.0[2]); + let (y6, y7) = unpack_pair(rhs.0[3]); + let (y8, y9) = unpack_pair(rhs.0[4]); let v19 = u32x8::new(19, 0, 19, 0, 19, 0, 19, 0); @@ -648,9 +747,9 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let y2_19 = m_lo(v19, y2); // iff 26 + b + lg(19) < 32 let y3_19 = m_lo(v19, y3); // if b < 32 - 26 - 4.248 = 1.752 let y4_19 = m_lo(v19, y4); - let y5_19 = m_lo(v19, y5); // below, b<2.5: this is a bottleneck, - let y6_19 = m_lo(v19, y6); // could be avoided by promoting to - let y7_19 = m_lo(v19, y7); // u64 here instead of in m() + let y5_19 = m_lo(v19, y5); + let y6_19 = m_lo(v19, y6); + let y7_19 = m_lo(v19, y7); let y8_19 = m_lo(v19, y8); let y9_19 = m_lo(v19, y9); @@ -671,6 +770,44 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let z8 = m(x0,y8) + m(x1_2,y7) + m(x2,y6) + m(x3_2,y5) + m(x4,y4) + m(x5_2,y3) + m(x6,y2) + m(x7_2,y1) + m(x8,y0) + m(x9_2,y9_19); let z9 = m(x0,y9) + m(x1,y8) + m(x2,y7) + m(x3,y6) + m(x4,y5) + m(x5,y4) + m(x6,y3) + m(x7,y2) + m(x8,y1) + m(x9,y0); + // The bounds on z[i] are the same as in the serial 32-bit code + // and the comment below is copied from there: + + // How big is the contribution to z[i+j] from x[i], y[j]? + // + // Using the bounds above, we get: + // + // i even, j even: x[i]*y[j] < 2^(26+b)*2^(26+b) = 2*2^(51+2*b) + // i odd, j even: x[i]*y[j] < 2^(25+b)*2^(26+b) = 1*2^(51+2*b) + // i even, j odd: x[i]*y[j] < 2^(26+b)*2^(25+b) = 1*2^(51+2*b) + // i odd, j odd: 2*x[i]*y[j] < 2*2^(25+b)*2^(25+b) = 1*2^(51+2*b) + // + // We perform inline reduction mod p by replacing 2^255 by 19 + // (since 2^255 - 19 = 0 mod p). This adds a factor of 19, so + // we get the bounds (z0 is the biggest one, but calculated for + // posterity here in case finer estimation is needed later): + // + // z0 < ( 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 249*2^(51 + 2*b) + // z1 < ( 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 154*2^(51 + 2*b) + // z2 < ( 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 195*2^(51 + 2*b) + // z3 < ( 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 118*2^(51 + 2*b) + // z4 < ( 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 141*2^(51 + 2*b) + // z5 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 82*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 87*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 )*2^(51 + 2b) = 46*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 )*2^(51 + 2b) = 33*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 )*2^(51 + 2b) = 10*2^(51 + 2*b) + // + // So z[0] fits into a u64 if 51 + 2*b + lg(249) < 64 + // if b < 2.5. + + // In fact this bound is slightly sloppy, since it treats both + // inputs x and y as being bounded by the same parameter b, + // while they are in fact bounded by b_x and b_y, and we + // already require that b_y < 1.75 in order to fit the + // multiplications by 19 into a u32. The tighter bound on b_y + // means we could get a tighter bound on the outputs, or a + // looser bound on b_x. FieldElement32x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9]) } } From 15f97221ba72fa800d463a7bbb75d6d4f5cf7e55 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 15 Jun 2018 13:44:31 -0700 Subject: [PATCH 23/27] Suppress extraneous warnings --- src/backend/avx2/edwards.rs | 5 ++--- src/backend/avx2/field.rs | 6 +++++- src/backend/avx2/scalar_mul/straus.rs | 2 ++ 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index 7c5ff70..a053414 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -10,8 +10,7 @@ //! Extended Twisted Edwards for Curve25519, using AVX2. -// just going to own it -#![allow(bad_style)] +#![allow(non_snake_case)] use core::convert::From; use core::ops::{Add, Neg, Sub}; @@ -393,7 +392,7 @@ mod test { } fn serial_double(P: edwards::EdwardsPoint) -> edwards::EdwardsPoint { - let (X1, Y1, Z1, T1) = (P.X, P.Y, P.Z, P.T); + let (X1, Y1, Z1, _T1) = (P.X, P.Y, P.Z, P.T); macro_rules! print_var { ($x:ident) => { diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 1484d75..b5bf4b3 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -10,16 +10,20 @@ //! 4-way vectorized 32bit field arithmetic using AVX2. -#![allow(bad_style)] +#![allow(non_snake_case)] const A_LANES: u8 = 0b0000_0101; const B_LANES: u8 = 0b0000_1010; const C_LANES: u8 = 0b0101_0000; const D_LANES: u8 = 0b1010_0000; +#[allow(unused)] const A_LANES64: u8 = 0b00_00_00_11; +#[allow(unused)] const B_LANES64: u8 = 0b00_00_11_00; +#[allow(unused)] const C_LANES64: u8 = 0b00_11_00_00; +#[allow(unused)] const D_LANES64: u8 = 0b11_00_00_00; use core::ops::{Add, Mul, Neg}; diff --git a/src/backend/avx2/scalar_mul/straus.rs b/src/backend/avx2/scalar_mul/straus.rs index 52c3079..053713f 100644 --- a/src/backend/avx2/scalar_mul/straus.rs +++ b/src/backend/avx2/scalar_mul/straus.rs @@ -8,6 +8,8 @@ // - Isis Agora Lovecruft // - Henry de Valence +#![allow(non_snake_case)] + use core::borrow::Borrow; use clear_on_drop::ClearOnDrop; From d791047aac8046daa3c8dc8e8684574523a22f6c Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 15 Jun 2018 16:35:35 -0700 Subject: [PATCH 24/27] Rewrite notes and documentation. --- docs/avx2-notes.md | 563 ++++++++++++++++++------------------ src/backend/avx2/edwards.rs | 48 ++- src/backend/avx2/field.rs | 55 +++- 3 files changed, 370 insertions(+), 296 deletions(-) diff --git a/docs/avx2-notes.md b/docs/avx2-notes.md index fb8085d..d902650 100644 --- a/docs/avx2-notes.md +++ b/docs/avx2-notes.md @@ -1,42 +1,48 @@ -An implementation of group operations on the twisted Edwards form of -Curve25519, using AVX2 to implement the 4-way parallel formulas of -Hisil, Wong, Carter, and Dawson (HWCD). -Their 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08], which -introduced the extended coordinates used in other parts of `-dalek`, -also describes 4-way parallel formulas for point addition and -doubling: +A vectorized implementation of group operations on the twisted Edwards +form of Curve25519, using a modification of the 4-way parallel +formulas of Hisil, Wong, Carter, and Dawson. -* a unified addition algorithm taking an effective \\(2\mathbf M + -1\mathbf D\\); +# Overview -* a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf -S\\); +The 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08] by Hisil, +Wong, Carter, and Dawson (HWCD) introduced the “extended coordinates” +and mixed-model representations which are used by most Edwards curve +implementations. -* a dedicated (i.e., for distinct points) addition algorithm taking -an effective \\(2 \mathbf M \\). +However, they also describe 4-way parallel formulas for point addition +and doubling: a unified addition algorithm taking an effective +\\(2\mathbf M + 1\mathbf D\\), a doubling algorithm taking an +effective \\(1\mathbf M + 1\mathbf S\\), and a dedicated (i.e., for +distinct points) addition algorithm taking an effective \\(2 \mathbf M +\\). They compare these formulas with a 2-way parallel variant of the +Montgomery ladder. -Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of -multiplication and squaring of generic field elements and \\(\mathbf -D\\) represents the cost of multiplication by a curve constant. +Unlike their serial formulas, which are used widely, their parallel +formulas do not seem to have been implemented in software before. The +2-way parallel Montgomery ladder was used in 2015 by Tung Chou's +`sandy2x` implementation. Curiously, however, although the [`sandy2x` +paper][sandy2x] also implements Edwards arithmetic, and cites HWCD08, +it doesn't mention their parallel Edwards formulas. +A 2015 paper by Hernández and López describes an AVX2 implementation +of X25519. Neither the paper nor the code are publicly available, but +it apparently gives only a [slight speedup][avx2trac], suggesting that +it uses a 4-way parallel Montgomery ladder rather than parallel +Edwards formulas. -These formulas do not seem to have been implemented using SIMD before. -A 2015 paper by Hernández and López mentions using AVX2 for the X25519 -Montgomery ladder, but neither the paper nor the code are publicly -available, and it apparently gives only a [slight speedup][avx2trac]. -The 2008 HWCD paper also describes and analyzes a 2-wide variant of the -Montgomery ladder (for comparison with parallel Edwards formulas); this -strategy was used in 2015 by Tung Chou's `sandy2x` implementation, which -used a 2-wide field implementation in 128-bit vector registers. -Curiously, however, although the [`sandy2x` paper][sandy2x] also -implements Edwards arithmetic, and cites the HWCD paper, it doesn't -mention the parallel formulas from HWCD, suggesting that they have been -overlooked for software implementations. +The reason may be that HWCD08 describe their formulas as operating on +four independent processors, which would make a software +implementation impractical: all of the operations are too low-latency +to effectively synchronize. But a closer inspection reveals that the +(more expensive) multiplication and squaring steps are uniform, while +the instruction divergence occurs in the (much cheaper) addition and +subtraction steps. This means that a SIMD implementation can perform +the expensive steps uniformly, and handle divergence in the +inexpensive steps using masking. -The notes below describe a tweak to the \\( 2\mathbf M + 1\mathbf D \\) -unified addition formulas to give \\( 2\mathbf M \\) readdition with -\\(1\mathbf D\\) precomputation, and a tweak to the doubling formulas to -avoid an extra reduction. These tweaked formulas are the ones used by -the `avx2` backend of `curve25519-dalek`. +These notes describe modifications to the original parallel formulas +to allow a SIMD implementation, and this module contains an +implementation of the modified formulas using 256-bit AVX2 vector +operations. # Parallel formulas in HWCD'08 @@ -60,218 +66,153 @@ and the unified addition algorithm is presented as follows: | | \\( R\_1 \gets R\_6 - R\_5 \\) | \\( R\_2 \gets R\_8 - R\_7 \\) | \\( R\_3 \gets R\_8 + R\_7 \\) | \\( R\_4 \gets R\_6 + R\_5 \\) | | \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_3 R\_4 \\) | \\( T\_3 \gets R\_1 R\_4 \\) | \\( Z\_3 \gets R\_2 R\_3 \\) | -Here \\( k = 2d \\) is a curve constant. +Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of +multiplication and squaring of generic field elements, \\(\mathbf D\\) +represents the cost of multiplication by a curve constant (in this +case \\( k = 2d \\)). -For a software implementation, each processor's operations are too -low-latency to parallelize across threads. However, the main cost -is in the multiplication and squaring steps, which are uniform, while -the divergent steps involve inexpensive additions and subtractions. +Notice that the \\(1\mathbf M\\) and \\(1\mathbf S\\) steps are +uniform. The non-uniform steps are all inexpensive additions or +subtractions, with the exception of the multiplication by the curve +constant \\(k = 2d\\): +$$ +R\_7 \gets 2 d R\_7. +$$ -This means we can use SIMD to implement the expensive portions in -parallel, and handle the instruction divergence on the inexpensive parts -using masking. - -The remaining obstacle to parallelism is the multiplication by the curve -constant \\(k = 2d\\). In the Curve25519 case, this is - -$$ k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. $$ - -HWCD suggest parallelising this step by breaking \\(k\\) into four +HWCD suggest parallelising this step by breaking \\(k = 2d\\) into four parts as \\(k = k_0 + 2\^n k_1 + 2\^{2n} k_2 + 2\^{3n} k_3 \\) and -computing \\(k_i R_7 \\) in parallel. However, this would be -somewhat awkward in our case, since we would normally represent -\\(k\\) as \\( 10 \\) 32-bit limbs, and \\(10 \\) is not divisible -by \\(4\\), so we would need a specialized routine to perform a -vectorized multiplication by 64-bit constants. +computing \\(k_i R_7 \\) in parallel. This is quite awkward, but if +the curve constant is a ratio \\( d = d\_1/d\_2 \\), then projective +coordinates allow us to instead compute +$$ +(R\_5, R\_6, R\_7, R\_8) \gets (d\_2 R\_5, d\_2 R\_6, 2d\_1 R\_7, d\_2 R\_8). +$$ +This can be performed as a uniform multiplication by a vector of +constants, and if \\(d\_1, d\_2\\) are small, it is relatively +inexpensive. (This trick was suggested by Mike Hamburg). +In the Curve25519 case, we have +$$ +d = \frac{d\_1}{d\_2} = \frac{-121665}{121666}; +$$ +Since \\(2 \cdot 121666 < 2\^{18}\\), all the constants above fit (up +to sign) in 32 bits, so this can be done in parallel as four +multiplications by small constants \\( (121666, 121666, 2\cdot 121665, +2\cdot 121666) \\), followed by a negation to compute \\( - 2\cdot 121665\\). -Instead, since we are working projectively, we can multiply -\\(R_7\\) by \\( -2\cdot 121665 \\) and multiply the other three -variables by \\(121666\\). This trick was suggested by Mike -Hamburg. Ignoring the sign for the moment, since -\\(2 \cdot 121666 < 2\^{18}\\), all these constants fit in 32 bits, -so (up to sign) this can be done in parallel as four multiplications -by small constants \\( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \\). +# Modified parallel formulas -How do we handle the sign? -Since we're primarily interested in Ristretto performance, not -Curve25519 performance, we could alternately work on the -\\(4\\)-isogenous "IsoEd25519" curve, which has \\(d = 121665\\). -However, this would only save the negation step, since multiplying -one field element by a 32-bit constant is not much easier than -multiplying four field elements by 32-bit constants, and it would -prevent accelerating Curve25519, so we don't make this choice. -Instead, we just negate one lane, and move the \\(1 \mathbf D\\) -into precomputation (see below). - -# Tweaked formulas - -After tweaking the formulas as described above, we obtain the -following. To avoid confusion with the original HWCD formulas, -temporary variables are named \\(S\\) instead of \\(R\\) and are in -static single-assignment form. +Using the modifications sketched above, we can write SIMD-friendly +versions of the parallel formulas as follows. To avoid confusion with +the original formulas, temporary variables are named \\(S\\) instead +of \\(R\\) and are in static single-assignment form. ## Addition -This implementation only implements readdition, but the tweaked addition -formulas are described first. To add points \\(P_1 = (X_1 : Y_1 : Z_1 : -T_1) \\) and \\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\), we compute - +To add points +\\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) +and +\\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\), +we compute $$ \begin{aligned} -S\_0 &\gets Y\_1 - X\_1 \\\\ -S\_1 &\gets Y\_1 + X\_1 \\\\ -S\_2 &\gets Y\_2 - X\_2 \\\\ -S\_3 &\gets Y\_2 + X\_2 +(S\_0 &&,&& S\_1 &&,&& S\_2 &&,&& S\_3 ) +&\gets +(Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Y\_2 - X\_2&&,&& Y\_2 + X\_2) +\\\\ +(S\_4 &&,&& S\_5 &&,&& S\_6 &&,&& S\_7 ) +&\gets +(S\_0 \cdot S\_2&&,&& S\_1 \cdot S\_3&&,&& Z\_1 \cdot Z\_2&&,&& T\_1 \cdot T\_2) +\\\\ +(S\_8 &&,&& S\_9 &&,&& S\_{10} &&,&& S\_{11} ) +&\gets +(d\_2 \cdot S\_4 &&,&& d\_2 \cdot S\_5 &&,&& 2 d\_2 \cdot S\_6 &&,&& 2 d\_1 \cdot S\_7 ) +\\\\ +(S\_{12} &&,&& S\_{13} &&,&& S\_{14} &&,&& S\_{15}) +&\gets +(S\_9 - S\_8&&,&& S\_9 + S\_8&&,&& S\_{10} - S\_{11}&&,&& S\_{10} + S\_{11}) +\\\\ +(X\_3&&,&& Y\_3&&,&& Z\_3&&,&& T\_3) +&\gets +(S\_{12} \cdot S\_{14}&&,&& S\_{15} \cdot S\_{13}&&,&& S\_{15} \cdot S\_{14}&&,&& S\_{12} \cdot S\_{13}) \end{aligned} $$ - -$$ -\begin{aligned} -S\_4 &\gets S\_0 S\_2 \\\\ -S\_5 &\gets S\_1 S\_3 \\\\ -S\_6 &\gets Z\_1 Z\_2 \\\\ -S\_7 &\gets T\_1 T\_2 -\end{aligned} -$$ - -$$ -\begin{aligned} -S\_8 &\gets S\_4 \cdot 121666 \\\\ -S\_9 &\gets S\_5 \cdot 121666 \\\\ -S\_{10} &\gets S\_6 \cdot 2 \cdot 121666 \\\\ -S\_{11} &\gets S\_7 \cdot -2 \cdot 121665 -\end{aligned} -$$ - -$$ -\begin{aligned} -S\_{12} &\gets S\_9 - S\_8 \\\\ -S\_{13} &\gets S\_9 + S\_8 \\\\ -S\_{14} &\gets S\_{10} - S\_{11} \\\\ -S\_{15} &\gets S\_{10} + S\_{11} -\end{aligned} -$$ - -$$ -\begin{aligned} -X\_3 &\gets S\_{12} S\_{14} \\\\ -Y\_3 &\gets S\_{15} S\_{13} \\\\ -Z\_3 &\gets S\_{15} S\_{14} \\\\ -T\_3 &\gets S\_{12} S\_{13} -\end{aligned} -$$ - to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). +This costs \\( 2\mathbf M + 1 \mathbf D\\). ## Readdition -If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we can precompute - +If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we +can cache the multiplication of the curve constants by computing $$ \begin{aligned} -S\_2 &\gets Y\_2 - X\_2 \\\\ -S\_3 &\gets Y\_2 + X\_2 +(S\_2' &&,&& S\_3' &&,&& Z\_2' &&,&& T\_2' ) +&\gets +(d\_2 \cdot (Y\_2 - X\_2)&&,&& d\_2 \cdot (Y\_1 + X\_1)&&,&& 2d\_2 \cdot Z\_2 &&,&& 2d\_1 \cdot T\_2). \end{aligned} $$ - +This costs \\( 1\mathbf D\\); with \\( (S\_2', S\_3', Z\_2', T\_2')\\) +in hand, the addition formulas above become $$ \begin{aligned} -S\_2' &\gets S\_2 \cdot 121666 \\\\ -S\_3' &\gets S\_3 \cdot 121666 \\\\ -Z\_2' &\gets Z\_2 \cdot 2 \cdot 121666 \\\\ -T\_2' &\gets T\_2 \cdot -2 \cdot 121665 \\\\ +(S\_0 &&,&& S\_1 &&,&& Z\_1 &&,&& T\_1 ) +&\gets +(Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Z\_1 &&,&& T\_1) +\\\\ +(S\_8 &&,&& S\_9 &&,&& S\_{10} &&,&& S\_{11} ) +&\gets +(S\_0 \cdot S\_2' &&,&& S\_1 \cdot S\_3'&&,&& Z\_1 \cdot Z\_2' &&,&& T\_1 \cdot T\_2') +\\\\ +(S\_{12} &&,&& S\_{13} &&,&& S\_{14} &&,&& S\_{15}) +&\gets +(S\_9 - S\_8&&,&& S\_9 + S\_8&&,&& S\_{10} - S\_{11}&&,&& S\_{10} + S\_{11}) +\\\\ +(X\_3&&,&& Y\_3&&,&& Z\_3&&,&& T\_3) +&\gets +(S\_{12} \cdot S\_{14}&&,&& S\_{15} \cdot S\_{13}&&,&& S\_{15} \cdot S\_{14}&&,&& S\_{12} \cdot S\_{13}) \end{aligned} $$ - -to obtain the `CachedPoint` \\( (S\_2', S\_3', Z\_2', T\_2') \\). -This precomputation is essentially the same as that suggested in -§3.1 of HWCD, with the difference that the multiplication by the curve -constant \\( -121665 / 121666 \\) is spread over all four -coordinates, to allow a vectorized computation of four -multiplications of small constants instead of a serial computation -of multiplication by a large constant. - -To perform readdition of \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and -\\(P_2 = (S\_2', S\_3', Z\_2', T\_2') \\), we compute - -$$ -\begin{aligned} -S\_0 &\gets Y\_1 - X\_1 \\\\ -S\_1 &\gets Y\_1 + X\_1 -\end{aligned} -$$ - -$$ -\begin{aligned} -S\_8 &\gets S\_0 S\_2' \\\\ -S\_9 &\gets S\_1 S\_3' \\\\ -S\_{10} &\gets Z\_1 Z\_2' \\\\ -S\_{11} &\gets T\_1 T\_2' -\end{aligned} -$$ - -$$ -\begin{aligned} -S\_{12} &\gets S\_9 - S\_8 \\\\ -S\_{13} &\gets S\_9 + S\_8 \\\\ -S\_{14} &\gets S\_{10} - S\_{11} \\\\ -S\_{15} &\gets S\_{10} + S\_{11} -\end{aligned} -$$ - -$$ -\begin{aligned} -X\_3 &\gets S\_{12} S\_{14} \\\\ -Y\_3 &\gets S\_{15} S\_{13} \\\\ -Z\_3 &\gets S\_{15} S\_{14} \\\\ -T\_3 &\gets S\_{12} S\_{13} -\end{aligned} -$$ - -to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). - -Compared to the addition formulas above, this saves \\( 1\mathbf D \\). +which costs only \\( 2\mathbf M \\). This precomputation is +essentially similar to the precomputation that HWCD suggest for their +serial formulas. Because the cost of precomputation and then +readdition is the same as addition, it's sufficient to only +implement caching and readdition. ## Doubling +The non-uniform portions of the (re)addition formulas have a fairly +regular structure. Unfortunately, this is not the case for the +doubling formulas, which are much less nice. + To double a point \\( P = (X\_1 : Y\_1 : Z\_1 : T\_1) \\), we compute - -$$ S\_0 \gets X\_1 + Y\_1 $$ - $$ \begin{aligned} -S\_1 &\gets X\_1\^2 \\\\ -S\_2 &\gets Y\_1\^2 \\\\ -S\_3 &\gets Z\_1\^2 \\\\ -S\_4 &\gets S\_0\^2 +(X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& S\_0) +&\gets +(X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& X\_1 + Y\_1) +\\\\ +(S\_1 &&,&& S\_2 &&,&& S\_3 &&,&& S\_4 ) +&\gets +(X\_1\^2 &&,&& Y\_1\^2&&,&& Z\_1\^2 &&,&& S\_0\^2) +\\\\ +(S\_5 &&,&& S\_6 &&,&& S\_8 &&,&& S\_9 ) +&\gets +(S\_1 + S\_2 &&,&& S\_1 - S\_2 &&,&& S\_1 + 2S\_3 - S\_2 &&,&& S\_1 + S\_2 - S\_4) +\\\\ +(X\_3 &&,&& Y\_3 &&,&& Z\_3 &&,&& T\_3 ) +&\gets +(S\_8 \cdot S\_9 &&,&& S\_5 \cdot S\_6 &&,&& S\_8 \cdot S\_6 &&,&& S\_5 \cdot S\_9) \end{aligned} $$ - -$$ -\begin{aligned} -S\_5 &\gets S\_1 + S\_2 \\\\ -S\_6 &\gets S\_1 - S\_2 \\\\ -S\_7 &\gets 2S\_3 \\\\ -S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\ -S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4 -\end{aligned} -$$ - -$$ -\begin{aligned} -X\_3 &\gets S\_8 S\_9 \\\\ -Y\_3 &\gets S\_5 S\_6 \\\\ -Z\_3 &\gets S\_8 S\_6 \\\\ -T\_3 &\gets S\_5 S\_9 -\end{aligned} -$$ - to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = [2]P\_1 \\). -Unlike the (re)addition formulas, the divergent parts of these formulas -are less nice. However, with some careful bounds-juggling, it is -possible to implement them without inserting extra carry chains, as -described below. +The intermediate step between the squaring and multiplication requires +a long chain of additions, but with some care and finesse, +described below, it is possible (in our case) to arrange this +computation without requiring an intermediate reduction. + +However, it does mean that the doubling formulas have proportionately +more vectorization overhead than the (re)addition formulas. The +effects of this are discussed in the comparison section below. # Field element representation @@ -329,39 +270,36 @@ much difficulty. Going the other direction, to extend this to AVX512, we could either run two point operations in parallel in lower and upper halves of the registers, or use 2-way parallelism within a field operation. -# Handling the Doubling Formulas +# Avoiding Overflow in Doubling -The non-parallel portion of the doubling formulas is +To analyze the size of the field element coefficients during the +computations, we can parameterize the bounds on the limbs of each +field element by \\( b \in \mathbb R \\) representing the excess bits +above that limb's radix, so that each limb is bounded by either +\\(2\^{25+b} \\) or \\( 2\^{26+b} \\), as appropriate. -$$ -\begin{aligned} -S\_5 &\gets S\_1 + S\_2 \\\\ -S\_6 &\gets S\_1 - S\_2 \\\\ -S\_7 &\gets 2S\_3 \\\\ -S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\ -S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4 -\end{aligned} -$$ - -Performing too many intermediate additions and subtractions grows -the bounds beyond what is allowed as input to multiplication, -forcing an extra carry pass. However, it is just possible to avoid -this by rearranging signs. - -Assume that the bounds on the limbs of each field element are -parameterized by \\( b \in \mathbb R \\) representing the excess -bits, so that each limb is bounded by either -\\( 2\^{25+b} \\) or \\( 2\^{26+b} \\). - -The multiplication routine requires that its inputs are bounded by +The multiplication routine requires that its inputs are bounded with \\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\) into 32 bits. Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\) when \\( x < 2\^{27.75} = 2\^{26 + 1.75} \\). However, this is only required for one of the inputs; the other can grow up to \\( b < 2.5 \\). -Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as +In addition, the multiplication and squaring routines do not +canonically reduce their outputs, but can leave some small uncarried +excesses, so that their reduced outputs are bounded with +\\( b < 0.007 \\). +The non-parallel portion of the doubling formulas is +$$ +\begin{aligned} +(S\_5 &&,&& S\_6 &&,&& S\_8 &&,&& S\_9 ) +&\gets +(S\_1 + S\_2 &&,&& S\_1 - S\_2 &&,&& S\_1 + 2S\_3 - S\_2 &&,&& S\_1 + S\_2 - S\_4) +\end{aligned} +$$ + +Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as $$ \begin{matrix} & S\_1 & S\_1 & S\_1 & S\_1 \\\\ @@ -374,24 +312,22 @@ $$ =& S\_5 & S\_6 & S\_8 & S\_9 \end{matrix} $$ - -results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for +results in bit-excesses \\( < (1.01, 1.60, 2.33, 2.01)\\) for \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute are then - $$ \begin{aligned} -X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.00) \\\\ -Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ -Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ -T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 2.00) +X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.01) \\\\ +Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\ +Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\ +T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 2.01) \end{aligned} $$ - -which are too large. However, if we flip the sign of \\( S\_4 = -S\_0\^2 \\) during squaring, so that we output \\(S\_4' = -S\_4 -\pmod p\\), then we can compute - +which are too large: it's not possible to arrange the multiplicands so +that one vector has \\(b < 2.5\\) and the other has \\( b < 1.75 \\). +However, if we flip the sign of \\( S\_4 = S\_0\^2 \\) during +squaring, so that we output \\(S\_4' = -S\_4 \pmod p\\), then we can +compute $$ \begin{matrix} & S\_1 & S\_1 & S\_1 & S\_1 \\\\ @@ -404,61 +340,120 @@ $$ =& S\_5 & S\_6 & S\_8 & S\_9 \end{matrix} $$ - -resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for +resulting in bit-excesses \\( < (1.01, 1.60, 2.33, 1.60)\\) for \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute are then - $$ \begin{aligned} -X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.59) \\\\ -Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ -Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ -T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 1.59) +X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.60) \\\\ +Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\ +Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\ +T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 1.60) \end{aligned} $$ - whose right-hand sides are all bounded with \\( b < 1.75 \\) and -whose left-hand sides are all bounded with \\( b < 2.5 \\). +whose left-hand sides are all bounded with \\( b < 2.5 \\), +so that we can avoid any intermediate reductions. # Comparison to non-vectorized formulas -HWCD also suggest using a mixed representation, passing between \\( -\mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\) -"projective" coordinates, where doubling is slightly cheaper (saving -about \\(\mathbf 1M\\). This approach is used for the -non-vectorized `u32` and `u64` backends, and more -details on the different coordinate systems can be found in the -`curve_models` module documentation. +In theory, the parallel Edwards formulas seem to allow a \\(4\\)-way +speedup from parallelism. However, an actual vectorized +implementation has several slowdowns that cut into this speedup. -This optimization is not compatible with the parallel formulas, which are -therefore slightly less efficient when counting the total number of -field multiplications and squarings. In particular, vectorized doublings -are less efficient than serial doublings. - -In addition, the parallel formulas can only use a \\( 32 \times 32 +First, the parallel formulas can only use a \\( 32 \times 32 \rightarrow 64 \\)-bit integer multiplier, so the speedup from vectorization must overcome the disadvantage of losing the \\( 64 -\times 64 \rightarrow 128\\)-bit (serial) integer multiplier. +\times 64 \rightarrow 128\\)-bit (serial) integer multiplier. The +effect of this slowdown is microarchitecture-dependent, since it +requires accounting for the total number of multiplications and +additions and their relative costs. In the future, it will probably +be possible to avoid this slowdown by using the `IFMA52` instructions, +whose parallelism is perfectly suited to these formulas. -When compiling with AVX512VL, LLVM is able to use the extra -`ymm16..ymm31` registers to reduce register pressure, and avoid -spills during field multiplication. This gives a small but -noticeable speedup. +Second, the parallel doubling formulas incur both a theoretical and +practical slowdown. The parallel formulas described above work on the +\\( \mathbb P\^3 \\) “extended” coordinates. The \\( \mathbb P\^2 \\) +model introduced earlier by [Bernstein, Birkner, Joye, Lange, and +Peters][bbjlp08] allows slightly faster doublings, so HWCD suggest +mixing coordinate systems while performing scalar multiplication +(attributing the idea to [a 1998 paper][cmo98] by Cohen, Miyagi, and +Ono). The \\( T \\) coordinate is not required for doublings, so when +doublings are followed by doublings, its computation can be skipped. +More details on this approach and the different coordinate systems can +be found in the [`curve_models` module documentation][curve_models]. -Another concern with AVX2 is that currently-available Intel processors -(particularly Skylake and Skylake-X microarchitectures) perform thermal -throttling when using wide vector instructions. For a mixed workload, +Unfortunately, this optimization is not compatible with the parallel +formulas, which cannot save time by skipping a single variable, so the +parallel doubling formulas do slightly more work when counting the +total number of field multiplications and squarings. + +In addition, the parallel doubling formulas have a less regular +pattern of additions and subtractions than the parallel addition +formulas, so the vectorization overhead is proportionately greater. +Both the parallel addition and parallel doubling formulas also require +some shuffling to rearrange data within the vectors, which places more +pressure on the shuffle unit than is desirable. + +This means that the speedup from using a vectorized implementation of +parallel Edwards formulas is likely to be greatest in applications +that do fewer doublings and more additions (like a large multiscalar +multiplication) rather than applications that do fewer additions and +more doublings (like a double-base scalar multiplication). + +Third, current Intel CPUs perform thermal throttling when using wide +vector instructions. A detailed description can be found in §15.26 of +[the Intel Optimization Manual][intel], but using wide vector +instructions prevents the core from operating at higher frequencies. +The core can return to the higher-frequency state after 2 +milliseconds, but this timer is reset every time high-power +instructions are used. + +Any speedup from vectorization therefore has to be weighed against a +slowdown for the next few million instructions. For a mixed workload, where point operations are interspersed with other tasks, this can -reduce overall performance. This probably means that this -implementation is not suitable for basic applications, like signatures, -but could still be worthwhile for complex applications, like -zero-knowledge proofs, which do enough work to make it worthwhile. +reduce overall performance. This implementation is therefore probably +not suitable for basic applications, like signatures, but is +worthwhile for complex applications, like zero-knowledge proofs, which +do sustained work. + +For this reason, the AVX2 backend is not enabled by default, but can +be selected using the `avx2_backend` feature. + +# Future work + +There are several directions for future improvement: + +* Using the vectorized field arithmetic code to parallelize across + point operations rather than within a single point operation. This + is less flexible, but would give a speedup both from allowing use of + the faster mixed-model arithmetic and from reducing shuffle + pressure. One approach in this direction would be to implement + batched scalar-point operations using vectors of points (AoSoA + layout). This less generally useful but would give a speedup for + Bulletproofs. + +* Extending the implementation to use the full width of AVX512, either + handling the extra parallelism internally to a single point + operation (by using a 2-way parallel implementation of field + arithmetic instead of a wordsliced one), or externally, + parallelizing across point operations. Internal parallelism would + be preferable but might require too much shuffle pressure. + +* Generalizing the implementation to non-AVX2 instructions, + particularly NEON. The current point arithmetic code is written in + terms of field element vectors, which are in turn implemented using + platform SIMD vectors. It should be possible to write an alternate + implementation of the `FieldElement32x4` using NEON without changing + the point arithmetic. NEON has 128-bit vectors rather than 256-bit + vectors, but this may still be worthwhile compared to a serial + implementation. -On AMD's Zen microarchitecture, thermal throttling is not a concern, -since AVX2 is implemented at half rate, so there is no penalty for mixed -workloads (but also no speedup). [sandy2x]: https://eprint.iacr.org/2015/943.pdf [avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28 [hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf +[curve_models]: https://doc-internal.dalek.rs/curve25519_dalek/curve_models/index.html +[bbjlp08]: https://eprint.iacr.org/2008/013 +[cmo98]: https://link.springer.com/content/pdf/10.1007%2F3-540-49649-1_6.pdf +[intel]: https://software.intel.com/sites/default/files/managed/9e/bc/64-ia-32-architectures-optimization-manual.pdf \ No newline at end of file diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index a053414..43f02f7 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -8,7 +8,29 @@ // - Isis Agora Lovecruft // - Henry de Valence -//! Extended Twisted Edwards for Curve25519, using AVX2. +//! Parallel Edwards Arithmetic for Curve25519. +//! +//! This module currently has two point types: +//! +//! * `ExtendedPoint`: a point stored in vector-friendly format, with +//! vectorized doubling and addition; +//! +//! * `CachedPoint`: used for readdition. +//! +//! Details on the formulas can be found in the documentation for the +//! parent `avx2` module. +//! +//! This API is designed to be safe: vectorized points can only be +//! created from serial points (which do validation on decompression), +//! and operations on valid points return valid points, so invalid +//! point states should be unrepresentable. +//! +//! This design goal is met, with one exception: the `Neg` +//! implementation for the `CachedPoint` performs a lazy negation, so +//! that subtraction can be efficiently implemented as a negation and +//! an addition. Repeatedly negating a `CachedPoint` will cause its +//! coefficients to grow and eventually overflow. Repeatedly negating +//! a point should not be necessary anyways. #![allow(non_snake_case)] @@ -26,7 +48,13 @@ use traits::Identity; use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle}; use backend::avx2::constants; -/// A point on Curve25519, represented in an AVX2-friendly format. +/// A point on Curve25519, using parallel Edwards formulas for curve +/// operations. +/// +/// # Invariant +/// +/// The coefficients of an `ExtendedPoint` are bounded with +/// \\( b < 0.007 \\). #[derive(Copy, Clone, Debug)] pub struct ExtendedPoint(pub(super) FieldElement32x4); @@ -67,6 +95,7 @@ impl Identity for ExtendedPoint { } impl ExtendedPoint { + /// Compute the double of this point. pub fn double(&self) -> ExtendedPoint { // Want to compute (X1 Y1 Z1 X1+Y1). // Not sure how to do this less expensively than computing @@ -138,6 +167,15 @@ impl ExtendedPoint { } /// A cached point with some precomputed variables used for readdition. +/// +/// # Warning +/// +/// It is not safe to negate this point more than once. +/// +/// # Invariant +/// +/// As long as the `CachedPoint` is not repeatedly negated, its +/// coefficients will be bounded with \\( b < 1.0 \\). #[derive(Copy, Clone, Debug)] pub struct CachedPoint(pub(super) FieldElement32x4); @@ -196,6 +234,12 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { /// Add an `ExtendedPoint` and a `CachedPoint`. fn add(self, other: &'b CachedPoint) -> ExtendedPoint { + // The coefficients of an `ExtendedPoint` are reduced after + // every operation. If the `CachedPoint` was negated, its + // coefficients grow by one bit. So on input, `self` is + // bounded with `b < 0.007` and `other` is bounded with + // `b < 1.0`. + let mut tmp = self.0; tmp = tmp.blend(tmp.diff_sum(), Lanes::AB); diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index b5bf4b3..175ff14 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -8,7 +8,20 @@ // - Isis Agora Lovecruft // - Henry de Valence -//! 4-way vectorized 32bit field arithmetic using AVX2. +//! An implementation of 4-way vectorized 32bit field arithmetic using +//! AVX2. +//! +//! The `FieldElement32x4` struct provides a vector of four field +//! elements, implemented using AVX2 operations. Its API is designed +//! to abstract away the platform-dependent details, so that point +//! arithmetic can be implemented only in terms of a vector of field +//! elements. +//! +//! At this level, the API is optimized for speed and not safety. The +//! `FieldElement32x4` does not always perform reductions. The pre- +//! and post-conditions on the bounds of the coefficients are +//! documented for each method, but it is the caller's responsibility +//! to ensure that there are no overflows. #![allow(non_snake_case)] @@ -32,7 +45,7 @@ use core::simd::{i32x8, u32x8, u64x4, IntoBits}; use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO}; use backend::u64::field::FieldElement64; -/// Unpack +/// Unpack 32-bit lanes into 64-bit lanes: /// ``` /// (a0, b0, a1, b1, c0, d0, c1, d1) /// ``` @@ -55,7 +68,7 @@ fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { (a, b) } -/// Repack +/// Repack 64-bit lanes into 32-bit lanes: /// ``` /// (a0, 0, b0, 0, c0, 0, d0, 0) /// (a1, 0, b1, 0, c1, 0, d1, 0) @@ -102,6 +115,10 @@ pub enum Lanes { } /// The `Shuffle` enum represents a shuffle of a `FieldElement32x4`. +/// +/// The enum variants are named by what they do to a vector \\( +/// (A,B,C,D) \\); for instance, `Shuffle::BADC` turns \\( (A, B, C, +/// D) \\) into \\( (B, A, D, C) \\). #[derive(Copy, Clone, Debug)] pub enum Shuffle { AAAA, @@ -116,7 +133,13 @@ pub enum Shuffle { ABDC, } -/// A vector of four field elements, in an AVX2-friendly format. +/// A vector of four field elements. +/// +/// Each operation on a `FieldElement32x4` has documented effects on +/// the bounds of the coefficients. This API is designed for speed +/// and not safety; it is the caller's responsibility to ensure that +/// the post-conditions of one operation are compatible with the +/// pre-conditions of the next. #[derive(Clone, Copy, Debug)] pub struct FieldElement32x4(pub(crate) [u32x8; 5]); @@ -134,6 +157,8 @@ impl ConditionallyAssignable for FieldElement32x4 { } impl FieldElement32x4 { + /// Split this vector into an array of four (serial) field + /// elements. pub fn split(&self) -> [FieldElement64; 4] { let mut out = [FieldElement64::zero(); 4]; for i in 0..5 { @@ -155,6 +180,11 @@ impl FieldElement32x4 { out } + /// Rearrange the elements of this vector according to `control`. + /// + /// The `control` parameter should be a compile-time constant, so + /// that when this function is inlined, LLVM is able to lower the + /// shuffle using an immediate. #[inline] pub fn shuffle(&self, control: Shuffle) -> FieldElement32x4 { #[inline(always)] @@ -190,8 +220,13 @@ impl FieldElement32x4 { ]) } + /// Blend `self` with `other`, taking lanes specified in `control` from `other`. + /// + /// The `control` parameter should be a compile-time constant, so + /// that this function can be inlined and LLVM can lower it to a + /// blend instruction using an immediate. #[inline] - pub fn blend(&self, b: FieldElement32x4, control: Lanes) -> FieldElement32x4 { + pub fn blend(&self, other: FieldElement32x4, control: Lanes) -> FieldElement32x4 { #[inline(always)] fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 { unsafe { @@ -254,11 +289,11 @@ impl FieldElement32x4 { } FieldElement32x4([ - blend_lanes(self.0[0], b.0[0], control), - blend_lanes(self.0[1], b.0[1], control), - blend_lanes(self.0[2], b.0[2], control), - blend_lanes(self.0[3], b.0[3], control), - blend_lanes(self.0[4], b.0[4], control), + blend_lanes(self.0[0], other.0[0], control), + blend_lanes(self.0[1], other.0[1], control), + blend_lanes(self.0[2], other.0[2], control), + blend_lanes(self.0[3], other.0[3], control), + blend_lanes(self.0[4], other.0[4], control), ]) } From 16f39c82e557ad0924cd839c622ec08d489c574e Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Mon, 18 Jun 2018 13:30:45 -0700 Subject: [PATCH 25/27] Remove `yolocrypto` from `avx2_backend` --- .travis.yml | 4 ++-- README.md | 13 ++++++++----- build.rs | 2 +- src/backend/mod.rs | 2 +- src/lib.rs | 2 +- 5 files changed, 13 insertions(+), 10 deletions(-) diff --git a/.travis.yml b/.travis.yml index cf6db51..50e9a3c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -10,7 +10,7 @@ env: # Tests the u64 backend - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std u64_backend' # Tests the avx2 backend - - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend yolocrypto' + - TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend' # Tests serde support and default feature selection - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='serde' # Tests building without std. We have to select a backend, so we select the one @@ -21,7 +21,7 @@ matrix: exclude: # Test the avx2 backend only on nightly - rust: stable - env: TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend yolocrypto' + env: TEST_COMMAND=test EXTRA_FLAGS='--no-default-features' FEATURES='std avx2_backend' # Test no_std only on nightly. - rust: stable env: TEST_COMMAND=build EXTRA_FLAGS=--no-default-features FEATURES='u32_backend' diff --git a/README.md b/README.md index 7dba401..0943fd7 100644 --- a/README.md +++ b/README.md @@ -71,15 +71,17 @@ Curve arithmetic is implemented using one of the following backends: * a `u32` backend using `u64` products; * a `u64` backend using `u128` products; -* an experimental AVX2 backend, available using the `yolocrypto` feature when - compiling for a target with `target_feature=+avx2`. +* an `avx2` backend using parallel formulas, available when compiling for a + target with `target_feature=+avx2`. By default the `u64` backend is selected. To select a specific backend, use: ```sh cargo build --no-default-features --features "std u32_backend" cargo build --no-default-features --features "std u64_backend" -cargo build --no-default-features --features "std avx2_backend yolocrypto" +cargo build --no-default-features --features "std avx2_backend" ``` +Crates using `curve25519-dalek` can either select a backend on behalf of their +users, or expose feature flags that control the `curve25519-dalek` backend. Benchmarks are run using [`criterion.rs`][criterion]: @@ -88,7 +90,7 @@ Benchmarks are run using [`criterion.rs`][criterion]: export RUSTFLAGS="-C target_cpu=native" cargo bench --no-default-features --features "std u32_backend" cargo bench --no-default-features --features "std u64_backend" -cargo bench --no-default-features --features "std avx2_backend yolocrypto" +cargo bench --no-default-features --features "std avx2_backend" ``` # Contributing @@ -117,7 +119,8 @@ to the Dalek race.* Portions of this library were originally a port of [Adam Langley's Golang ed25519 library](https://github.com/agl/ed25519), which was in -turn a port of the reference `ref10` implementation. +turn a port of the reference `ref10` implementation. Most of this code, +including the 32-bit field arithmetic, has since been rewritten. The fast `u32` and `u64` scalar arithmetic was implemented by Andrew Moon, and the addition chain for scalar inversion was provided by Brian Smith. diff --git a/build.rs b/build.rs index 64aec5d..e0a2da7 100644 --- a/build.rs +++ b/build.rs @@ -1,5 +1,5 @@ #![cfg_attr(feature = "nightly", feature(cfg_target_feature))] -#![cfg_attr(all(feature = "nightly", feature = "yolocrypto"), feature(stdsimd))] +#![cfg_attr(all(feature = "nightly", feature = "avx2_backend"), feature(stdsimd))] #![allow(unused_variables)] #![allow(non_snake_case)] #![allow(dead_code)] diff --git a/src/backend/mod.rs b/src/backend/mod.rs index b095eef..aa44d52 100644 --- a/src/backend/mod.rs +++ b/src/backend/mod.rs @@ -27,6 +27,6 @@ pub mod u32; #[cfg(feature = "u64_backend")] pub mod u64; -#[cfg(all(feature = "avx2_backend", feature = "yolocrypto", target_feature = "avx2"))] +#[cfg(all(feature = "avx2_backend", target_feature = "avx2"))] pub mod avx2; diff --git a/src/lib.rs b/src/lib.rs index 4ebffe7..70f80f5 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -14,7 +14,7 @@ #![cfg_attr(feature = "nightly", feature(cfg_target_feature))] #![cfg_attr(feature = "nightly", feature(external_doc))] -#![cfg_attr(all(feature = "nightly", feature = "yolocrypto"), feature(stdsimd))] +#![cfg_attr(all(feature = "nightly", feature = "avx2_backend"), feature(stdsimd))] // Refuse to compile if documentation is missing, but only on nightly. // From 791506a29fb46bd779065ce5e643062f2ee22598 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 21 Jun 2018 12:01:09 -0700 Subject: [PATCH 26/27] Add avx2_backend to make doc FEATURES --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 8b376cb..e447175 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -FEATURES := nightly yolocrypto +FEATURES := nightly yolocrypto avx2_backend doc: cargo rustdoc --features "$(FEATURES)" -- --html-in-header docs/assets/rustdoc-include-katex-header.html From 28b7ed5709a4b127f02ea08d21a7cc06f563ee8a Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 22 Jun 2018 12:24:20 -0700 Subject: [PATCH 27/27] Change version to 0.18.0 (since feature flags changed) --- Cargo.toml | 2 +- README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 1103faf..3978ee5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "curve25519-dalek" -version = "0.17.0" +version = "0.18.0" authors = ["Isis Lovecruft ", "Henry de Valence "] readme = "README.md" diff --git a/README.md b/README.md index 0943fd7..ddf4283 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ make doc-internal To import `curve25519-dalek`, add the following to the dependencies section of your project's `Cargo.toml`: ```toml -curve25519-dalek = "^0.17" +curve25519-dalek = "^0.18" ``` Then import the crate as: ```rust,no_run