ed448-goldilocks: optimize Edwards isogeny map (#1316)

See also #1349

Author: daxpedda

ed448-goldilocks/src/curve/twedwards/extended.rs

@@ -104,25 +104,33 @@ impl ExtendedPoint {
         AffinePoint { x, y }
     }
 
-    /// Edwards_Isogeny is derived from the doubling formula
-    /// XXX: There is a duplicate method in the twisted edwards module to compute the dual isogeny
-    /// XXX: Not much point trying to make it generic I think. So what we can do is optimise each respective isogeny method for a=1 or a = -1 (currently, I just made it really slow and simple)
-    fn edwards_isogeny(&self, a: FieldElement) -> EdwardsExtendedPoint {
+    /// Uses a 2-isogeny to map the point to the Ed448-Goldilocks
+    // (2. Algorithm 2) https://eprint.iacr.org/2008/522.pdf
+    pub fn to_untwisted(self) -> EdwardsExtendedPoint {
+        // x = 2xy / (y^2 + a*x^2)
+        // y = (y^2 - a*x^2) / (2 - y^2 - a*x^2)
+
         // Convert to affine now, then derive extended version later
         let affine = self.to_affine();
         let x = affine.x;
         let y = affine.y;
 
+        let yy = y.square();
+        let xx = x.square();
+        let axx = -xx;
+        let yy_plus_axx = yy + axx;
+
         // Compute x
-        let xy = x * y;
-        let x_numerator = xy.double();
-        let x_denom = y.square() - (a * x.square());
-        let new_x = x_numerator * x_denom.invert();
+        let x_numerator = (x * y).double();
+        let x_denom = yy - axx;
 
         // Compute y
-        let y_numerator = y.square() + (a * x.square());
-        let y_denom = (FieldElement::ONE + FieldElement::ONE) - y.square() - (a * x.square());
-        let new_y = y_numerator * y_denom.invert();
+        let y_numerator = yy_plus_axx;
+        let y_denom = FieldElement::TWO - yy_plus_axx;
+
+        let common_denom = (x_denom * y_denom).invert();
+        let new_x = x_numerator * y_denom * common_denom;
+        let new_y = y_numerator * x_denom * common_denom;
 
         EdwardsExtendedPoint {
             X: new_x,
@@ -132,11 +140,6 @@ impl ExtendedPoint {
         }
     }
 
-    /// Uses a 2-isogeny to map the point to the Ed448-Goldilocks
-    pub fn to_untwisted(self) -> EdwardsExtendedPoint {
-        self.edwards_isogeny(FieldElement::MINUS_ONE)
-    }
-
     /// Checks if the point is on the curve
     pub(crate) fn is_on_curve(&self) -> Choice {
         let XY = self.X * self.Y;

ed448-goldilocks/src/edwards/extended.rs

@@ -669,38 +669,31 @@ impl EdwardsPoint {
         AffinePoint { x, y }
     }
 
-    /// Edwards_Isogeny is derived from the doubling formula
-    /// XXX: There is a duplicate method in the twisted edwards module to compute the dual isogeny
-    /// XXX: Not much point trying to make it generic I think. So what we can do is optimise each respective isogeny method for a=1 or a = -1 (currently, I just made it really slow and simple)
-    fn edwards_isogeny(&self, a: FieldElement) -> TwistedExtendedPoint {
-        // Convert to affine now, then derive extended version later
-        let affine = self.to_affine();
-        let x = affine.x;
-        let y = affine.y;
-
-        // Compute x
-        let xy = x * y;
-        let x_numerator = xy.double();
-        let x_denom = y.square() - (a * x.square());
-        let new_x = x_numerator * x_denom.invert();
-
-        // Compute y
-        let y_numerator = y.square() + (a * x.square());
-        let y_denom = (FieldElement::ONE + FieldElement::ONE) - y.square() - (a * x.square());
-        let new_y = y_numerator * y_denom.invert();
+    // Copied from https://github.com/otrv4/libgoldilocks/blob/d07cb5b423995bae1155702aa949846c95d855c1/src/goldilocks.c#L980-L994.
+    pub(crate) fn to_twisted(self) -> TwistedExtendedPoint {
+        let c = self.X.square();
+        let a = self.Y.square();
+        let d = c + a;
+        let t = self.Y + self.X;
+        let b = t.square();
+        let b = b - d;
+        let t = a - c;
+        let x = self.Z.square();
+        let z = x.double();
+        let a = z - d;
+        let new_x = a * b;
+        let new_z = t * a;
+        let new_y = t * d;
+        let new_t = b * d;
 
         TwistedExtendedPoint {
             X: new_x,
             Y: new_y,
-            Z: FieldElement::ONE,
-            T: new_x * new_y,
+            Z: new_z,
+            T: new_t,
         }
     }
 
-    pub(crate) fn to_twisted(self) -> TwistedExtendedPoint {
-        self.edwards_isogeny(FieldElement::ONE)
-    }
-
     /// Compute the negation of this point's `x`-coordinate.
     pub fn negate(&self) -> Self {
         EdwardsPoint {

ed448-goldilocks/src/field/element.rs

@@ -246,6 +246,7 @@ impl FieldElement {
         "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000262a8",
     )));
     pub const ONE: Self = Self(ConstMontyType::new(&U448::ONE));
+    pub const TWO: Self = Self(ConstMontyType::new(&U448::from_u64(2)));
     pub const TWISTED_D: Self = Self(ConstMontyType::new(&U448::from_be_hex(
         "fffffffffffffffffffffffffffffffffffffffffffffffffffffffeffffffffffffffffffffffffffffffffffffffffffffffffffff6755",
     )));