Add support for divisors for Ed25519

crypto/dalek-ff-group/src/lib.rs

@@ -344,7 +344,7 @@ macro_rules! dalek_group {
   ) => {
     /// Wrapper around the dalek Point type. For Ed25519, this is restricted to the prime subgroup.
     #[derive(Clone, Copy, PartialEq, Eq, Debug, Zeroize)]
-    pub struct $Point($DPoint);
+    pub struct $Point(pub $DPoint);
     deref_borrow!($Point, $DPoint);
     constant_time!($Point, $DPoint);
     math_neg!($Point, Scalar, $DPoint::add, $DPoint::sub, $DPoint::mul);

crypto/divisors/Cargo.toml

@@ -25,4 +25,6 @@ group = "0.13"
 [dev-dependencies]
 rand_core = { version = "0.6", features = ["getrandom"] }
 
+hex = "0.4"
+dalek-ff-group = { path = "../dalek-ff-group", features = ["std"] }
 pasta_curves = { version = "0.5", default-features = false, features = ["bits", "alloc"], git = "https://github.com/kayabaNerve/pasta_curves.git" }

crypto/divisors/src/lib.rs

@@ -19,24 +19,25 @@ where
   /// An element of the field this curve is defined over.
   type FieldElement: PrimeField;
 
-  /// The A in the curve equation y^2 = x^3 + A x^2 + B x + C.
-  const A: u64;
-  /// The B in the curve equation y^2 = x^3 + A x^2 + B x + C.
-  const B: u64;
-  /// The C in the curve equation y^2 = x^3 + A x^2 + B x + C.
-  const C: u64;
-
-  /// y^2 - x^3 - A x^2 - B x - C
+  /// The A in the curve equation y^2 = x^3 + A x + B.
+  fn a() -> Self::FieldElement;
+  /// The B in the curve equation y^2 = x^3 + A x + B.
+  fn b() -> Self::FieldElement;
+
+  /// y^2 - x^3 - A x - B
   fn divisor_modulus() -> Poly<Self::FieldElement> {
     Poly {
       y_coefficients: vec![Self::FieldElement::ZERO, Self::FieldElement::ONE],
       yx_coefficients: vec![],
       x_coefficients: vec![
-        -Self::FieldElement::from(Self::B),
-        -Self::FieldElement::from(Self::A),
+        // A x
+        -Self::a(),
+        // x^2
+        Self::FieldElement::ZERO,
+        // x^3
         -Self::FieldElement::ONE,
       ],
-      zero_coefficient: -Self::FieldElement::from(Self::C),
+      zero_coefficient: -Self::b(),
     }
   }
 
@@ -47,7 +48,9 @@ where
 // Calculate the slope and intercept for two points.
 fn slope_intercept<C: DivisorCurve>(a: C, b: C) -> (C::FieldElement, C::FieldElement) {
   let (ax, ay) = C::to_xy(a);
+  debug_assert_eq!(C::divisor_modulus().eval(ax, ay), C::FieldElement::ZERO);
   let (bx, by) = C::to_xy(b);
+  debug_assert_eq!(C::divisor_modulus().eval(bx, by), C::FieldElement::ZERO);
   let slope = (by - ay) *
     Option::<C::FieldElement>::from((bx - ax).invert())
       .expect("trying to get slope/intercept of points sharing an x coordinate");

crypto/divisors/src/tests/mod.rs

@@ -1,6 +1,10 @@
 use rand_core::OsRng;
 
-use group::{ff::Field, Group, Curve};
+use group::{
+  ff::{Field, PrimeField},
+  Group, GroupEncoding, Curve,
+};
+use dalek_ff_group::EdwardsPoint;
 use pasta_curves::{
   arithmetic::{Coordinates, CurveAffine},
   Ep, Fp,
@@ -11,9 +15,12 @@ use crate::{DivisorCurve, Poly, new_divisor};
 impl DivisorCurve for Ep {
   type FieldElement = Fp;
 
-  const A: u64 = 0;
-  const B: u64 = 0;
-  const C: u64 = 5;
+  fn a() -> Self::FieldElement {
+    Self::FieldElement::ZERO
+  }
+  fn b() -> Self::FieldElement {
+    Self::FieldElement::from(5u64)
+  }
 
   fn to_xy(point: Self) -> (Self::FieldElement, Self::FieldElement) {
     Option::<Coordinates<_>>::from(point.to_affine().coordinates())
@@ -22,39 +29,89 @@ impl DivisorCurve for Ep {
   }
 }
 
-type F = <Ep as DivisorCurve>::FieldElement;
+impl DivisorCurve for EdwardsPoint {
+  type FieldElement = dalek_ff_group::FieldElement;
+
+  // Wei25519 a/b
+  // https://www.ietf.org/archive/id/draft-ietf-lwig-curve-representations-02.pdf E.3
+  fn a() -> Self::FieldElement {
+    let mut be_bytes =
+      hex::decode("2aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa984914a144").unwrap();
+    be_bytes.reverse();
+    let le_bytes = be_bytes;
+    Self::FieldElement::from_repr(le_bytes.try_into().unwrap()).unwrap()
+  }
+  fn b() -> Self::FieldElement {
+    let mut be_bytes =
+      hex::decode("7b425ed097b425ed097b425ed097b425ed097b425ed097b4260b5e9c7710c864").unwrap();
+    be_bytes.reverse();
+    let le_bytes = be_bytes;
 
-#[test]
-fn test_divisor() {
+    Self::FieldElement::from_repr(le_bytes.try_into().unwrap()).unwrap()
+  }
+
+  // https://www.ietf.org/archive/id/draft-ietf-lwig-curve-representations-02.pdf E.2
+  fn to_xy(point: Self) -> (Self::FieldElement, Self::FieldElement) {
+    // Extract the y coordinate from the compressed point
+    let mut edwards_y = point.to_bytes();
+    let x_is_odd = edwards_y[31] >> 7;
+    edwards_y[31] &= (1 << 7) - 1;
+    let edwards_y = Self::FieldElement::from_repr(edwards_y).unwrap();
+
+    // Recover the x coordinate
+    let edwards_y_sq = edwards_y * edwards_y;
+    let D =
+      -Self::FieldElement::from(121665u64) * Self::FieldElement::from(121666u64).invert().unwrap();
+    let mut edwards_x = ((edwards_y_sq - Self::FieldElement::ONE) *
+      ((D * edwards_y_sq) + Self::FieldElement::ONE).invert().unwrap())
+    .sqrt()
+    .unwrap();
+    if u8::from(bool::from(edwards_x.is_odd())) != x_is_odd {
+      edwards_x = -edwards_x;
+    }
+
+    // Calculate the x and y coordinates for Wei25519
+    let edwards_y_plus_one = Self::FieldElement::ONE + edwards_y;
+    let one_minus_edwards_y = Self::FieldElement::ONE - edwards_y;
+    let wei_x = (edwards_y_plus_one * one_minus_edwards_y.invert().unwrap()) +
+      (Self::FieldElement::from(486662u64) * Self::FieldElement::from(3u64).invert().unwrap());
+    let c =
+      (-(Self::FieldElement::from(486662u64) + Self::FieldElement::from(2u64))).sqrt().unwrap();
+    let wei_y = c * edwards_y_plus_one * (one_minus_edwards_y * edwards_x).invert().unwrap();
+    (wei_x, wei_y)
+  }
+}
+
+fn test_divisor<C: DivisorCurve>() {
   for i in 1 ..= 255 {
     dbg!(i);
 
     // Select points
     let mut points = vec![];
     for _ in 0 .. i {
-      points.push(Ep::random(&mut OsRng));
+      points.push(C::random(&mut OsRng));
     }
-    points.push(-points.iter().sum::<Ep>());
+    points.push(-points.iter().sum::<C>());
 
     // Create the divisor
-    let divisor = new_divisor::<Ep>(&points).unwrap();
+    let divisor = new_divisor::<C>(&points).unwrap();
 
     // Decide challgenges
-    let c0 = Ep::random(&mut OsRng);
-    let c1 = Ep::random(&mut OsRng);
+    let c0 = C::random(&mut OsRng);
+    let c1 = C::random(&mut OsRng);
     let c2 = -(c0 + c1);
-    let (slope, intercept) = crate::slope_intercept::<Ep>(c0, c1);
+    let (slope, intercept) = crate::slope_intercept::<C>(c0, c1);
 
     // Perform the original check
     {
       let eval = |c| {
-        let (x, y) = Ep::to_xy(c);
+        let (x, y) = C::to_xy(c);
         divisor.eval(x, y)
       };
 
-      let mut rhs = F::ONE;
+      let mut rhs = C::FieldElement::ONE;
       for point in &points {
-        let (x, y) = Ep::to_xy(*point);
+        let (x, y) = C::to_xy(*point);
         rhs *= intercept - (y - (slope * x));
       }
       assert_eq!(eval(c0) * eval(c1) * eval(c2), rhs);
@@ -66,15 +123,15 @@ fn test_divisor() {
         let dx = Poly {
           y_coefficients: vec![],
           yx_coefficients: vec![],
-          x_coefficients: vec![F::ZERO, F::from(3)],
-          zero_coefficient: F::from(<Ep as DivisorCurve>::A),
+          x_coefficients: vec![C::FieldElement::ZERO, C::FieldElement::from(3)],
+          zero_coefficient: C::a(),
         };
 
         let dy = Poly {
-          y_coefficients: vec![F::from(2)],
+          y_coefficients: vec![C::FieldElement::from(2)],
           yx_coefficients: vec![],
           x_coefficients: vec![],
-          zero_coefficient: F::ZERO,
+          zero_coefficient: C::FieldElement::ZERO,
         };
 
         let dz = (dy.clone() * -slope) + &dx;
@@ -86,23 +143,23 @@ fn test_divisor() {
 
       {
         let sanity_eval = |c| {
-          let (x, y) = Ep::to_xy(c);
+          let (x, y) = C::to_xy(c);
           dx_over_dz.0.eval(x, y) * dx_over_dz.1.eval(x, y).invert().unwrap()
         };
         let sanity = sanity_eval(c0) + sanity_eval(c1) + sanity_eval(c2);
         // This verifies the dx/dz polynomial is correct
-        assert_eq!(sanity, F::ZERO);
+        assert_eq!(sanity, C::FieldElement::ZERO);
       }
 
       // Logarithmic derivative check
       let test = |divisor: Poly<_>| {
         let (dx, dy) = divisor.differentiate();
 
         let lhs = |c| {
-          let (x, y) = Ep::to_xy(c);
+          let (x, y) = C::to_xy(c);
 
-          let n_0 = (F::from(3) * (x * x)) + F::from(<Ep as DivisorCurve>::A);
-          let d_0 = (F::from(2) * y).invert().unwrap();
+          let n_0 = (C::FieldElement::from(3) * (x * x)) + C::a();
+          let d_0 = (C::FieldElement::from(2) * y).invert().unwrap();
           let nd_0 = n_0 * d_0;
 
           let n_1 = dy.eval(x, y);
@@ -122,9 +179,9 @@ fn test_divisor() {
         };
         let lhs = lhs(c0) + lhs(c1) + lhs(c2);
 
-        let mut rhs = F::ZERO;
+        let mut rhs = C::FieldElement::ZERO;
         for point in &points {
-          let (x, y) = <Ep as DivisorCurve>::to_xy(*point);
+          let (x, y) = <C as DivisorCurve>::to_xy(*point);
           rhs += (intercept - (y - (slope * x))).invert().unwrap();
         }
 
@@ -137,27 +194,69 @@ fn test_divisor() {
   }
 }
 
-#[test]
-fn test_same_point() {
-  let mut points = vec![Ep::random(&mut OsRng)];
+fn test_same_point<C: DivisorCurve>() {
+  let mut points = vec![C::random(&mut OsRng)];
   points.push(points[0]);
-  points.push(-points.iter().sum::<Ep>());
+  points.push(-points.iter().sum::<C>());
 
-  let divisor = new_divisor::<Ep>(&points).unwrap();
+  let divisor = new_divisor::<C>(&points).unwrap();
   let eval = |c| {
-    let (x, y) = Ep::to_xy(c);
+    let (x, y) = C::to_xy(c);
     divisor.eval(x, y)
   };
 
-  let c0 = Ep::random(&mut OsRng);
-  let c1 = Ep::random(&mut OsRng);
+  let c0 = C::random(&mut OsRng);
+  let c1 = C::random(&mut OsRng);
   let c2 = -(c0 + c1);
-  let (slope, intercept) = crate::slope_intercept::<Ep>(c0, c1);
+  let (slope, intercept) = crate::slope_intercept::<C>(c0, c1);
 
-  let mut rhs = F::ONE;
+  let mut rhs = <C as DivisorCurve>::FieldElement::ONE;
   for point in points {
-    let (x, y) = Ep::to_xy(point);
+    let (x, y) = C::to_xy(point);
     rhs *= intercept - (y - (slope * x));
   }
   assert_eq!(eval(c0) * eval(c1) * eval(c2), rhs);
 }
+
+#[test]
+fn test_divisor_pallas() {
+  test_divisor::<Ep>();
+  test_same_point::<Ep>();
+}
+
+#[test]
+fn test_divisor_ed25519() {
+  // Since we're implementing Wei25519 ourselves, check the isomorphism works as expected
+  {
+    let incomplete_add = |p1, p2| {
+      let (x1, y1) = EdwardsPoint::to_xy(p1);
+      let (x2, y2) = EdwardsPoint::to_xy(p2);
+
+      // mmadd-1998-cmo
+      let u = y2 - y1;
+      let uu = u * u;
+      let v = x2 - x1;
+      let vv = v * v;
+      let vvv = v * vv;
+      let R = vv * x1;
+      let A = uu - vvv - R.double();
+      let x3 = v * A;
+      let y3 = (u * (R - A)) - (vvv * y1);
+      let z3 = vvv;
+
+      // Normalize from XYZ to XY
+      let x3 = x3 * z3.invert().unwrap();
+      let y3 = y3 * z3.invert().unwrap();
+
+      // Edwards addition -> Wei25519 coordinates should be equivalent to Wei25519 addition
+      assert_eq!(EdwardsPoint::to_xy(p1 + p2), (x3, y3));
+    };
+
+    for _ in 0 .. 256 {
+      incomplete_add(EdwardsPoint::random(&mut OsRng), EdwardsPoint::random(&mut OsRng));
+    }
+  }
+
+  test_divisor::<EdwardsPoint>();
+  test_same_point::<EdwardsPoint>();
+}