Implement point subtraction

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();
+            // set t2 = (S8 S9 S10 S11)
+            t2.scale_by_curve_constants(true);
 
-            for i in 0..5 {
+            // set t2 = (S8 S9 S11 S10)
-                let swapped = _mm256_shuffle_epi32(t2.0[i].into(), 0b10_11_00_01);
+            t2.swap_CD();
-                t2.0[i] = _mm256_blend_epi32(t2.0[i].into(), swapped, 0b11110000).into();
-            }
 
+            // 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();
+        // Test the serial implementation of the parallel addition formulas
-        let R2: edwards::ExtendedPoint = (&ExtendedPoint::from(P) + &ExtendedPoint::from(Q)).into();
+        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!("R = P + Q = {:?}", &P + &Q);
-        println!("(vector) R2 = {:?}", R2);
+        println!("R_serial = {:?}", R_serial);
-        println!("P + Q = {:?}", &P + &Q);
+        println!("R_vector = {:?}", R_vector);
-        assert_eq!(R1.compress(), (&P + &Q).compress());
+        println!("S = P - Q = {:?}", &P - &Q);
-        assert_eq!(R2.compress(), (&P + &Q).compress());
+        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);

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();
+            if negate_121665 {
-            b[1] = _mm256_blend_epi32(b1.into(), (odd__p20 - b1).into(), 0b11_00_00_00).into();
+                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();
+            if negate_121665 {
-            b[3] = _mm256_blend_epi32(b3.into(), (odd__p20 - b3).into(), 0b11_00_00_00).into();
+                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();
+            if negate_121665 {
-            b[5] = _mm256_blend_epi32(b5.into(), (odd__p20 - b5).into(), 0b11_00_00_00).into();
+                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();
+            if negate_121665 {
-            b[7] = _mm256_blend_epi32(b7.into(), (odd__p20 - b7).into(), 0b11_00_00_00).into();
+                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();
+            if negate_121665 {
-            b[9] = _mm256_blend_epi32(b9.into(), (odd__p20 - b9).into(), 0b11_00_00_00).into();
+                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]