Implement improved algorithms from paper

Paper: Complete Addition Law for Montgomery Curves

Author: daxpedda

ed448-goldilocks/src/montgomery/ops.rs

@@ -8,37 +8,44 @@ use subtle::{Choice, ConditionallySelectable};
 
 use super::{MontgomeryPoint, MontgomeryScalar, ProjectiveMontgomeryPoint};
 
+#[allow(non_upper_case_globals)]
+const s: FieldElement = FieldElement(ConstMontyType::new(&U448::from_u64(3)));
+
 impl Add<&ProjectiveMontgomeryPoint> for &ProjectiveMontgomeryPoint {
     type Output = ProjectiveMontgomeryPoint;
 
-    // Copied from https://github.com/armfazh/redox-ecc/blob/5a8c09c5ef9fe6a8d2c749d05eca011c6d661599/src/montgomery/point.rs#L80-L104.
+    // See Complete Addition Law for Montgomery Curves - Algorithm 1.
+    // With "Trade-Off Technique". 
     fn add(self, rhs: &ProjectiveMontgomeryPoint) -> Self::Output {
-        const S: FieldElement = FieldElement(ConstMontyType::new(&U448::from_u64(3)));
-
         let (x1, y1, z1) = (self.U, self.V, self.W);
         let (x2, y2, z2) = (rhs.U, rhs.V, rhs.W);
-        let (a_ec, s_ec) = (FieldElement::J, S);
-        let (t0, t1, t2) = (x1 * x2, y1 * y2, z1 * z2);
-        let (t3, t4) = (x1 * y2, x2 * y1);
-        let (t5, t6) = (y1 * z2, y2 * z1);
-        let (t7, t8) = (x1 * z2, x2 * z1);
-        let t9 = t7 + t8;
-        let ta = t9 + (t0 * a_ec);
-        let rr = t5 + t6;
-        let tt = ta - t1;
-        let vv = t9 * a_ec + t0.double() + t0 + t2;
-        let ss = (t3 - t4) * s_ec + t0 - t2;
-        let uu = (t7 - t8) * s_ec - t3 - t4;
-        let ww = (t5 - t6) * s_ec + ta + t1;
-        let x3 = rr * ss - tt * uu;
-        let y3 = tt * ww - vv * ss;
-        let z3 = vv * uu - rr * ww;
 
-        ProjectiveMontgomeryPoint {
-            U: x3,
-            V: y3,
-            W: z3,
-        }
+        let t0 = x1 * x2;
+        let t1 = y1 * y2;
+        let t2 = z1 * z2;
+        let t3 = x1 * y2;
+        let t4 = x2 * y1;
+        let t5 = y1 * z2;
+        let t6 = y2 * z1;
+        let t7 = x1 * z2;
+        let t8 = x2 * z1;
+        let t9 = t7 + t8;
+        let t10 = t9 + FieldElement::J * t0;
+        let R = t5 + t6;
+        let T = t10 - t1;
+        let V = FieldElement::J * t9 + t0.double() + t0 + t2;
+        let S = s * (t3 - t4) + t0 - t2;
+        let U = s * (t7 - t8) - t3 - t4;
+        let W = s * (t5 - t6) + t10 + t1;
+        let C = (R + T) * (S - U);
+        let D = (R - T) * (S + U);
+        let E = (T + V) * (W - S);
+        let F = (T - V) * (W + S);
+        let X = C + D;
+        let Y = E + F;
+        let Z = (U - W).double() * (R + V) + C - D + E - F;
+
+        ProjectiveMontgomeryPoint { U: X, V: Y, W: Z }
     }
 }
 
@@ -51,8 +58,38 @@ define_add_variants!(
 impl Add<&MontgomeryPoint> for &ProjectiveMontgomeryPoint {
     type Output = ProjectiveMontgomeryPoint;
 
-    fn add(self, other: &MontgomeryPoint) -> ProjectiveMontgomeryPoint {
-        *self + *other
+    // See Complete Addition Law for Montgomery Curves - Algorithm 1.
+    // With "Trade-Off Technique". 
+    fn add(self, rhs: &MontgomeryPoint) -> ProjectiveMontgomeryPoint {
+        let (x1, y1, z1) = (self.U, self.V, self.W);
+        let (x2, y2) = (rhs.x, rhs.y);
+
+        let t0 = x1 * x2;
+        let t1 = y1 * y2;
+        let t2 = z1;
+        let t3 = x1 * y2;
+        let t4 = x2 * y1;
+        let t5 = y1;
+        let t6 = y2 * z1;
+        let t7 = x1;
+        let t8 = x2 * z1;
+        let t9 = t7 + t8;
+        let t10 = t9 + FieldElement::J * t0;
+        let R = t5 + t6;
+        let T = t10 - t1;
+        let V = FieldElement::J * t9 + t0.double() + t0 + t2;
+        let S = s * (t3 - t4) + t0 - t2;
+        let U = s * (t7 - t8) - t3 - t4;
+        let W = s * (t5 - t6) + t10 + t1;
+        let C = (R + T) * (S - U);
+        let D = (R - T) * (S + U);
+        let E = (T + V) * (W - S);
+        let F = (T - V) * (W + S);
+        let X = C + D;
+        let Y = E + F;
+        let Z = (U - W).double() * (R + V) + C - D + E - F;
+
+        ProjectiveMontgomeryPoint { U: X, V: Y, W: Z }
     }
 }
 
@@ -66,7 +103,7 @@ impl Add<&ProjectiveMontgomeryPoint> for &MontgomeryPoint {
     type Output = ProjectiveMontgomeryPoint;
 
     fn add(self, other: &ProjectiveMontgomeryPoint) -> ProjectiveMontgomeryPoint {
-        *other + *self
+        other + self
     }
 }
 
@@ -111,9 +148,9 @@ impl Mul<&MontgomeryScalar> for &ProjectiveMontgomeryPoint {
         let mut p = ProjectiveMontgomeryPoint::IDENTITY;
         let bits = scalar.bits();
 
-        for s in (0..448).rev() {
+        for index in (0..448).rev() {
             p = p + p;
-            p.conditional_assign(&(p + self), Choice::from(bits[s] as u8));
+            p.conditional_assign(&(p + self), Choice::from(bits[index] as u8));
         }
 
         p
@@ -264,3 +301,20 @@ where
         iter.fold(Self::IDENTITY, |acc, item| acc + item.borrow())
     }
 }
+
+#[cfg(test)]
+mod test {
+    use elliptic_curve::Group;
+    use rand_core::OsRng;
+
+    use super::*;
+
+    #[test]
+    fn mixed_addition() {
+        let p1 = ProjectiveMontgomeryPoint::try_from_rng(&mut OsRng).unwrap();
+        let p2 = ProjectiveMontgomeryPoint::try_from_rng(&mut OsRng).unwrap();
+        let p3 = p1 + p2;
+
+        assert_eq!(p3.to_affine(), (p1.to_affine() + p2).to_affine());
+    }
+}

ed448-goldilocks/src/montgomery/point.rs

@@ -25,8 +25,8 @@ use super::{
 /// A point in Montgomery form including the y-coordinate.
 #[derive(Copy, Clone, Debug, Default)]
 pub struct MontgomeryPoint {
-    x: FieldElement,
-    y: FieldElement,
+    pub(super) x: FieldElement,
+    pub(super) y: FieldElement,
 }
 
 impl MontgomeryPoint {
@@ -160,10 +160,6 @@ impl ProjectiveMontgomeryPoint {
         Self { U, V, W }
     }
 
-    fn double(&self) -> Self {
-        self + self
-    }
-
     /// Convert the point to its form without the y-coordinate
     pub fn to_projective_x(&self) -> ProjectiveMontgomeryXpoint {
         ProjectiveMontgomeryXpoint::conditional_select(
@@ -281,7 +277,46 @@ impl Group for ProjectiveMontgomeryPoint {
     }
 
     fn double(&self) -> Self {
-        self.double()
+        self + self
+
+        // TODO: figure out why this doesn't work.
+        // See Complete Addition Law for Montgomery Curves - Algorithm 3.
+        // const A_MINUS_1: FieldElement = FieldElement(ConstMontyType::new(&U448::from_u64(156325)));
+
+        // let (x, y, z) = (self.U, self.V, self.W);
+
+        // let t0 = x.square();
+        // let t1 = y.square();
+        // let t2 = z.square();
+        // let t3 = (x + y).square();
+        // let t4 = (y + z).square();
+        // let t5 = (x + z).square();
+        // let t6 = t1 + t2;
+        // let t7 = (t0 - t2).double();
+        // let t8 = A_MINUS_1 * t0;
+        // let t9 = t0 - t1;
+        // let t10 = FieldElement::J * (t5 - t2) + t0 + t9;
+        // let t11 = t5 + t8;
+        // let t12 = t10 - t6.double();
+        // let SU_ = t3 - t6;
+        // let SU = -SU_ + t7;
+        // let RT_ = t4 - t11;
+        // let RT = t4 + t12;
+        // let WS_ = t11 - t9;
+        // let WS = WS_ + t7;
+        // let TV_ = t12 - t10 + t8;
+        // let TV = t5 + t10;
+        // let UW_ = SU - WS;
+        // let RV = RT_ + TV;
+        // let C = RT * SU_;
+        // let D = RT_ * SU;
+        // let E = TV * WS_;
+        // let F = TV_ * WS;
+        // let X = C + D;
+        // let Y = E + F;
+        // let Z = UW_.double() * RV + C - D + E - F;
+
+        // Self { U: X, V: Y, W: Z }
     }
 }