diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index a4dd1ac..0ebbc52 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -14,7 +14,7 @@ #![allow(bad_style)] use std::convert::From; -use std::ops::{Add, Mul, Neg}; +use std::ops::{Add, Sub, Mul, Neg}; use stdsimd::simd::{u32x8, i32x8}; @@ -70,6 +70,7 @@ impl<'a> Neg for &'a ExtendedPoint { fn neg(self) -> ExtendedPoint { let mut neg = *self; + // (X Y Z T) -> (-X Y Z -T) neg.0.mask_negate(0b10100101); neg } @@ -216,36 +217,114 @@ impl<'a, 'b> Add<&'b ExtendedPoint> for &'a ExtendedPoint { let mut t0 = FieldElement32x4::zero(); let mut t1 = FieldElement32x4::zero(); + // set t0 = (X1 Y1 X2 Y2) for i in 0..5 { t0.0[i] = _mm256_permute2x128_si256(P.0[i].into(), Q.0[i].into(), 32).into(); } + // set t0 = (Y1-X1 Y1+X1 Y2-X2 Y2+X2) = (S0 S1 S2 S3) t0.diff_sum(); + // set t1 = (S0 S1 Z1 T1) + // set t0 = (S2 S3 Z2 T2) for i in 0..5 { t1.0[i] = _mm256_blend_epi32(t0.0[i].into(), P.0[i].into(), 0b11110000).into(); t0.0[i] = _mm256_permute2x128_si256(t0.0[i].into(), Q.0[i].into(), 49).into(); } + // set t2 = (S0*S2 S1*S3 Z1*Z2 T1*T2) = (S4 S5 S6 S7) let mut t2 = &t0 * &t1; - - t2.scale_by_curve_constants(); - - for i in 0..5 { - let swapped = _mm256_shuffle_epi32(t2.0[i].into(), 0b10_11_00_01); - t2.0[i] = _mm256_blend_epi32(t2.0[i].into(), swapped, 0b11110000).into(); - } + // set t2 = (S8 S9 S10 S11) + t2.scale_by_curve_constants(true); + + // set t2 = (S8 S9 S11 S10) + t2.swap_CD(); + + // set t2 = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) t2.diff_sum(); 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) for i in 0..5 { t0.0[i] = _mm256_permutevar8x32_epi32(t2.0[i], c0); t1.0[i] = _mm256_permutevar8x32_epi32(t2.0[i], c1); } + // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + ExtendedPoint(&t0 * &t1) + } + } +} + +impl<'a, 'b> Sub<&'b ExtendedPoint> for &'a ExtendedPoint { + type Output = ExtendedPoint; + + /// Uses a slight tweak of the parallel unified formulas of HWCD'08 + fn sub(self, other: &'b ExtendedPoint) -> ExtendedPoint { + unsafe { + use stdsimd::vendor::_mm256_permute2x128_si256; + use stdsimd::vendor::_mm256_permutevar8x32_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + use stdsimd::vendor::_mm256_shuffle_epi32; + + let P: &FieldElement32x4 = &self.0; + let Q: &FieldElement32x4 = &other.0; + + let mut t0 = FieldElement32x4::zero(); + let mut t1 = FieldElement32x4::zero(); + + // set t0 = (X1 Y1 X2 Y2) + for i in 0..5 { + t0.0[i] = _mm256_permute2x128_si256(P.0[i].into(), Q.0[i].into(), 32).into(); + } + + // Since we're subtracting instead of adding, we want to add the point (-X2 Y2 Z2 -T2). + // Set (X2' Y2' Z2' T2') = (-X2 Y2 Z2 -T2) + // + // so S2 = Y2 - X2' = Y2 - (-X2) = Y2 + X2 + // and S3 = Y2 + X2' = Y2 + (-X2) = Y2 - X2 + + // set t0 = (Y1-X1 Y1+X1 Y2-X2 Y2+X2) = (S0 S1 S3 S2) + t0.diff_sum(); + + // set t0 = (S0 S1 S2 S3) + t0.swap_CD(); + + // set t1 = (S0 S1 Z1 T1) + // set t0 = (S2 S3 Z2 T2) = (S2 S3 Z2' -T2') + for i in 0..5 { + t1.0[i] = _mm256_blend_epi32(t0.0[i].into(), P.0[i].into(), 0b11110000).into(); + t0.0[i] = _mm256_permute2x128_si256(t0.0[i].into(), Q.0[i].into(), 49).into(); + } + + // set t2 = (S0*S2 S1*S3 Z1*Z2 T1*T2 ) + // = (S0*S2 S1*S3 Z1*Z2' -T1*T2') = (S4 S5 S6 -S7) + let mut t2 = &t0 * &t1; + + // set t2 = (S8 S9 S10 S11) + t2.scale_by_curve_constants(false); + + // set t2 = (S8 S9 S11 S10) + t2.swap_CD(); + + // set t2 = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) + t2.diff_sum(); + + 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) + for i in 0..5 { + t0.0[i] = _mm256_permutevar8x32_epi32(t2.0[i], c0); + t1.0[i] = _mm256_permutevar8x32_epi32(t2.0[i], c1); + } + + // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) ExtendedPoint(&t0 * &t1) } } @@ -502,7 +581,7 @@ pub mod vartime { Q = &Q + &odd_multiple[( naf[i]/2) as usize]; } else if naf[i] < 0 { // XXX impl Sub - Q = &Q + &(-&odd_multiple[(-naf[i]/2) as usize]); + Q = &Q - &odd_multiple[(-naf[i]/2) as usize]; } } } @@ -577,40 +656,65 @@ mod test { } fn addition_test_helper(P: edwards::ExtendedPoint, Q: edwards::ExtendedPoint) { - let R1: edwards::ExtendedPoint = serial_add(P.into(), Q.into()).into(); - let R2: edwards::ExtendedPoint = (&ExtendedPoint::from(P) + &ExtendedPoint::from(Q)).into(); + // Test the serial implementation of the parallel addition formulas + let R_serial: edwards::ExtendedPoint = serial_add(P.into(), Q.into()).into(); + // Test the vector implementation of the parallel addition formulas + let R_vector: edwards::ExtendedPoint = (&ExtendedPoint::from(P) + &ExtendedPoint::from(Q)).into(); + // Test the vector implementation of the parallel subtraction formulas + let S_vector: edwards::ExtendedPoint = (&ExtendedPoint::from(P) - &ExtendedPoint::from(Q)).into(); + println!("Testing point addition:"); println!("P = {:?}", P); println!("Q = {:?}", Q); - println!("(serial) R1 = {:?}", R1); - println!("(vector) R2 = {:?}", R2); - println!("P + Q = {:?}", &P + &Q); - assert_eq!(R1.compress(), (&P + &Q).compress()); - assert_eq!(R2.compress(), (&P + &Q).compress()); + println!("R = P + Q = {:?}", &P + &Q); + println!("R_serial = {:?}", R_serial); + println!("R_vector = {:?}", R_vector); + println!("S = P - Q = {:?}", &P - &Q); + println!("S_vector = {:?}", S_vector); + assert_eq!(R_serial.compress(), (&P + &Q).compress()); + assert_eq!(R_vector.compress(), (&P + &Q).compress()); + assert_eq!(S_vector.compress(), (&P - &Q).compress()); println!("OK!\n"); } + #[test] + fn sub_vs_add_minus() { + let P: ExtendedPoint = edwards::ExtendedPoint::identity().into(); + let Q: ExtendedPoint = edwards::ExtendedPoint::identity().into(); + + let mQ = -&Q; + + println!("sub"); + let R1: edwards::ExtendedPoint = (&P - &Q).into(); + println!("add neg"); + let R2: edwards::ExtendedPoint = (&P + &mQ).into(); + + assert_eq!(R2.compress(), edwards::ExtendedPoint::identity().compress()); + assert_eq!(R1.compress(), edwards::ExtendedPoint::identity().compress()); + } + + #[test] fn vector_addition_vs_serial_addition_vs_edwards_extendedpoint() { use constants; use scalar::Scalar; - println!("Testing id + id"); + println!("Testing id +- id"); let P = edwards::ExtendedPoint::identity(); let Q = edwards::ExtendedPoint::identity(); addition_test_helper(P, Q); - println!("Testing id + B"); + println!("Testing id +- B"); let P = edwards::ExtendedPoint::identity(); let Q = constants::ED25519_BASEPOINT_POINT; addition_test_helper(P, Q); - println!("Testing B + B"); + println!("Testing B +- B"); let P = constants::ED25519_BASEPOINT_POINT; let Q = constants::ED25519_BASEPOINT_POINT; addition_test_helper(P, Q); - println!("Testing B + kB"); + println!("Testing B +- kB"); let P = constants::ED25519_BASEPOINT_POINT; let Q = &constants::ED25519_BASEPOINT_TABLE * &Scalar::from_u64(8475983829); addition_test_helper(P, Q); diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index b911c1c..b76d454 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -114,6 +114,18 @@ impl FieldElement32x4 { self.reduce32(); } + // Given `self = (A,B,C,D)`, set `self = (A,B,D,C)` + pub fn swap_CD(&mut self) { + unsafe { + use stdsimd::vendor::_mm256_shuffle_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + for i in 0..5 { + let swapped = _mm256_shuffle_epi32(self.0[i].into(), 0b10_11_00_01); + self.0[i] = _mm256_blend_epi32(self.0[i].into(), swapped, 0b11110000).into(); + } + } + } + // Given `self = (A,B,C,D)`, set `self = (B - A, B + A, D - C, D + C)`. pub fn diff_sum(&mut self) { /// (v0 v1 v2 v3 v4 v5 v6 v7) -> (v1 v0 v3 v2 v5 v4 v7 v6) @@ -188,7 +200,17 @@ impl FieldElement32x4 { FieldElement32x4(out) } - pub fn scale_by_curve_constants(&mut self) { + /// Let `self` \\(= (A, B, C, D) \\). + /// + /// If `negate_121665 = true`, compute + /// + /// $$( 121666A, 121666B, 2\cdot 121666C, -2\cdot 121665 D).$$ + /// + /// If `negate_121665 = false`, compute + /// + /// $$( 121666A, 121666B, 2\cdot 121666C, 2\cdot 121665 D).$$ + /// + pub fn scale_by_curve_constants(&mut self, negate_121665: bool) { let mut b = [u64x4::splat(0); 10]; let consts = u32x8::new(121666, 0, 121666, 0, 2*121666, 0, 2*121665, 0); @@ -203,32 +225,57 @@ impl FieldElement32x4 { let (b0, b1) = unpack_pair(self.0[0]); let b0 = _mm256_mul_epu32(b0, consts); // need a new binding since now let b1 = _mm256_mul_epu32(b1, consts); // b0 has type u64x4 - b[0] = _mm256_blend_epi32(b0.into(), (low__p20 - b0).into(), 0b11_00_00_00).into(); - b[1] = _mm256_blend_epi32(b1.into(), (odd__p20 - b1).into(), 0b11_00_00_00).into(); + if negate_121665 { + b[0] = _mm256_blend_epi32(b0.into(), (low__p20 - b0).into(), 0b11_00_00_00).into(); + b[1] = _mm256_blend_epi32(b1.into(), (odd__p20 - b1).into(), 0b11_00_00_00).into(); + } else { + b[0] = b0; + b[1] = b1; + } let (b2, b3) = unpack_pair(self.0[1]); let b2 = _mm256_mul_epu32(b2, consts); let b3 = _mm256_mul_epu32(b3, consts); - b[2] = _mm256_blend_epi32(b2.into(), (even_p20 - b2).into(), 0b11_00_00_00).into(); - b[3] = _mm256_blend_epi32(b3.into(), (odd__p20 - b3).into(), 0b11_00_00_00).into(); + if negate_121665 { + b[2] = _mm256_blend_epi32(b2.into(), (even_p20 - b2).into(), 0b11_00_00_00).into(); + b[3] = _mm256_blend_epi32(b3.into(), (odd__p20 - b3).into(), 0b11_00_00_00).into(); + } else { + b[2] = b2; + b[3] = b3; + } let (b4, b5) = unpack_pair(self.0[2]); let b4 = _mm256_mul_epu32(b4, consts); let b5 = _mm256_mul_epu32(b5, consts); - b[4] = _mm256_blend_epi32(b4.into(), (even_p20 - b4).into(), 0b11_00_00_00).into(); - b[5] = _mm256_blend_epi32(b5.into(), (odd__p20 - b5).into(), 0b11_00_00_00).into(); + if negate_121665 { + b[4] = _mm256_blend_epi32(b4.into(), (even_p20 - b4).into(), 0b11_00_00_00).into(); + b[5] = _mm256_blend_epi32(b5.into(), (odd__p20 - b5).into(), 0b11_00_00_00).into(); + } else { + b[4] = b4; + b[5] = b5; + } let (b6, b7) = unpack_pair(self.0[3]); let b6 = _mm256_mul_epu32(b6, consts); let b7 = _mm256_mul_epu32(b7, consts); - b[6] = _mm256_blend_epi32(b6.into(), (even_p20 - b6).into(), 0b11_00_00_00).into(); - b[7] = _mm256_blend_epi32(b7.into(), (odd__p20 - b7).into(), 0b11_00_00_00).into(); + if negate_121665 { + b[6] = _mm256_blend_epi32(b6.into(), (even_p20 - b6).into(), 0b11_00_00_00).into(); + b[7] = _mm256_blend_epi32(b7.into(), (odd__p20 - b7).into(), 0b11_00_00_00).into(); + } else { + b[6] = b6; + b[7] = b7; + } let (b8, b9) = unpack_pair(self.0[4]); let b8 = _mm256_mul_epu32(b8, consts); let b9 = _mm256_mul_epu32(b9, consts); - b[8] = _mm256_blend_epi32(b8.into(), (even_p20 - b8).into(), 0b11_00_00_00).into(); - b[9] = _mm256_blend_epi32(b9.into(), (odd__p20 - b9).into(), 0b11_00_00_00).into(); + if negate_121665 { + b[8] = _mm256_blend_epi32(b8.into(), (even_p20 - b8).into(), 0b11_00_00_00).into(); + b[9] = _mm256_blend_epi32(b9.into(), (odd__p20 - b9).into(), 0b11_00_00_00).into(); + } else { + b[8] = b8; + b[9] = b9; + } } *self = FieldElement32x4::reduce64(b); @@ -517,13 +564,22 @@ mod test { #[test] fn scale_by_curve_constants() { let mut x = FieldElement32x4::splat(&FieldElement64::one()); - x.scale_by_curve_constants(); + x.scale_by_curve_constants(true); 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])); + + let mut y = FieldElement32x4::splat(&FieldElement64::one()); + y.scale_by_curve_constants(false); + + let ys = y.split(); + assert_eq!(ys[0], FieldElement64([ 121666,0,0,0,0])); + assert_eq!(ys[1], FieldElement64([ 121666,0,0,0,0])); + assert_eq!(ys[2], FieldElement64([2*121666,0,0,0,0])); + assert_eq!(ys[3], FieldElement64([2*121665,0,0,0,0])); } #[test]