From 879ec566b993c8c063c72f0a29c53739efed3158 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Sun, 20 Jul 2025 00:28:33 +0200
Subject: [PATCH 01/10] Add optimized Edwards addition and doubling algorithms

---
 ed448-goldilocks/src/edwards/extended.rs | 60 ++++++++++++------------
 1 file changed, 30 insertions(+), 30 deletions(-)

diff --git a/ed448-goldilocks/src/edwards/extended.rs b/ed448-goldilocks/src/edwards/extended.rs
index 9b0da4adb..c6af6a79c 100644
--- a/ed448-goldilocks/src/edwards/extended.rs
+++ b/ed448-goldilocks/src/edwards/extended.rs
@@ -597,41 +597,41 @@ impl EdwardsPoint {
     }
 
     /// Add two points
-    //https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf (3.1)
-    // These formulas are unified, so for now we can use it for doubling. Will refactor later for speed
+    // (3.1) https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
     pub fn add(&self, other: &EdwardsPoint) -> Self {
-        let aXX = self.X * other.X; // aX1X2
-        let dTT = FieldElement::EDWARDS_D * self.T * other.T; // dT1T2
-        let ZZ = self.Z * other.Z; // Z1Z2
-        let YY = self.Y * other.Y;
-
-        let X = {
-            let x_1 = (self.X * other.Y) + (self.Y * other.X);
-            let x_2 = ZZ - dTT;
-            x_1 * x_2
-        };
-        let Y = {
-            let y_1 = YY - aXX;
-            let y_2 = ZZ + dTT;
-            y_1 * y_2
-        };
-
-        let T = {
-            let t_1 = YY - aXX;
-            let t_2 = (self.X * other.Y) + (self.Y * other.X);
-            t_1 * t_2
-        };
-
-        let Z = { (ZZ - dTT) * (ZZ + dTT) };
-
-        EdwardsPoint { X, Y, Z, T }
+        let A = self.X * other.X;
+        let B = self.Y * other.Y;
+        let C = self.T * other.T * FieldElement::EDWARDS_D;
+        let D = self.Z * other.Z;
+        let E = (self.X + self.Y) * (other.X + other.Y) - A - B;
+        let F = D - C;
+        let G = D + C;
+        let H = B - A;
+        Self {
+            X: E * F,
+            Y: G * H,
+            Z: F * G,
+            T: E * H,
+        }
     }
 
     /// Double this point
-    // XXX: See comment on addition, the formula is unified, so this will do for now
-    //https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf (3.1)
+    // (3.3) https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
     pub fn double(&self) -> Self {
-        self.add(self)
+        let A = self.X.square();
+        let B = self.Y.square();
+        let C = self.Z.square().double();
+        let D = A;
+        let E = (self.X + self.Y).square() - A - B;
+        let G = D + B;
+        let F = G - C;
+        let H = D - B;
+        Self {
+            X: E * F,
+            Y: G * H,
+            Z: F * G,
+            T: E * H,
+        }
     }
 
     /// Check if this point is on the curve

From b120a04238e8de2117dd823fbff56cd0d8acac39 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Sat, 26 Jul 2025 23:57:33 +0200
Subject: [PATCH 02/10] Reduce unnecessary conversions to `ExtensiblePoint` and
 back

---
 .../src/curve/scalar_mul/double_and_add.rs    |   7 +-
 .../src/curve/scalar_mul/variable_base.rs     |  13 +-
 .../src/curve/scalar_mul/window/wnaf.rs       |  16 +-
 .../src/curve/twedwards/affine.rs             |  11 +-
 .../src/curve/twedwards/extended.rs           | 169 +++++++++++-------
 .../src/curve/twedwards/extensible.rs         | 155 ++++++----------
 .../src/curve/twedwards/projective.rs         |  18 +-
 ed448-goldilocks/src/decaf/ops.rs             |   6 +-
 ed448-goldilocks/src/decaf/points.rs          |  26 +--
 ed448-goldilocks/src/edwards/extended.rs      |  14 +-
 10 files changed, 217 insertions(+), 218 deletions(-)

diff --git a/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs b/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs
index aac590611..917ba16d5 100644
--- a/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs
@@ -1,9 +1,10 @@
 use crate::curve::twedwards::extended::ExtendedPoint;
+use crate::curve::twedwards::extensible::ExtensiblePoint;
 use subtle::{Choice, ConditionallySelectable};
 
 /// Traditional double and add algorithm
-pub(crate) fn double_and_add(point: &ExtendedPoint, s_bits: [bool; 448]) -> ExtendedPoint {
-    let mut result = ExtendedPoint::IDENTITY;
+pub(crate) fn double_and_add(point: &ExtendedPoint, s_bits: [bool; 448]) -> ExtensiblePoint {
+    let mut result = ExtensiblePoint::IDENTITY;
 
     // NB, we reverse here, so we are going from MSB to LSB
     // XXX: Would be great if subtle had a From<u32> for Choice. But maybe that is not it's purpose?
@@ -12,7 +13,7 @@ pub(crate) fn double_and_add(point: &ExtendedPoint, s_bits: [bool; 448]) -> Exte
 
         let mut p = ExtendedPoint::IDENTITY;
         p.conditional_assign(point, Choice::from(bit as u8));
-        result = result.add(&p);
+        result = result.to_extended().add_extended(&p);
     }
 
     result
diff --git a/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs b/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
index bdbca2b79..4c420f52f 100644
--- a/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
@@ -5,7 +5,7 @@ use crate::EdwardsScalar;
 use crate::curve::twedwards::{extended::ExtendedPoint, extensible::ExtensiblePoint};
 use subtle::{Choice, ConditionallyNegatable};
 
-pub fn variable_base(point: &ExtendedPoint, s: &EdwardsScalar) -> ExtendedPoint {
+pub fn variable_base(point: &ExtendedPoint, s: &EdwardsScalar) -> ExtensiblePoint {
     let mut result = ExtensiblePoint::IDENTITY;
 
     // Recode Scalar
@@ -28,10 +28,10 @@ pub fn variable_base(point: &ExtendedPoint, s: &EdwardsScalar) -> ExtendedPoint
         let mut neg_P = lookup.select(abs_value);
         neg_P.conditional_negate(Choice::from((sign) as u8));
 
-        result = result.add_projective_niels(&neg_P);
+        result = result.to_extended().add_projective_niels(&neg_P);
     }
 
-    result.to_extended()
+    result
 }
 
 #[cfg(test)]
@@ -55,7 +55,7 @@ mod test {
         assert_eq!(got, got2);
 
         // Lets see if this is conserved over the isogenies
-        let edwards_point = twisted_point.to_untwisted();
+        let edwards_point = twisted_point.to_extensible().to_untwisted();
         let got_untwisted_point = edwards_point.scalar_mul(&scalar);
         let expected_untwisted_point = got.to_untwisted();
         assert_eq!(got_untwisted_point, expected_untwisted_point);
@@ -69,9 +69,8 @@ mod test {
         let exp = variable_base(&x, &EdwardsScalar::from(1u8));
         assert!(x == exp);
         // Test that 2 * (P + P) = 4 * P
-        let x_ext = x.to_extensible();
-        let expected_two_x = x_ext.add_extensible(&x_ext).double();
+        let expected_two_x = x.add_extended(&x).double();
         let got = variable_base(&x, &EdwardsScalar::from(4u8));
-        assert!(expected_two_x.to_extended() == got);
+        assert!(expected_two_x == got);
     }
 }
diff --git a/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs b/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs
index 1c804fded..73c1c651f 100644
--- a/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs
@@ -6,13 +6,14 @@ pub struct LookupTable([ProjectiveNielsPoint; 8]);
 
 /// Precomputes odd multiples of the point passed in
 impl From<&ExtendedPoint> for LookupTable {
-    fn from(point: &ExtendedPoint) -> LookupTable {
-        let P = point.to_extensible();
-
+    fn from(P: &ExtendedPoint) -> LookupTable {
         let mut table = [P.to_projective_niels(); 8];
 
         for i in 1..8 {
-            table[i] = P.add_projective_niels(&table[i - 1]).to_projective_niels();
+            table[i] = P
+                .add_projective_niels(&table[i - 1])
+                .to_extended()
+                .to_projective_niels();
         }
 
         LookupTable(table)
@@ -42,11 +43,8 @@ fn test_lookup() {
     let mut expected_point = ExtendedPoint::IDENTITY;
     for i in 0..8 {
         let selected_point = points.select(i);
-        assert_eq!(selected_point.to_extended(), expected_point);
+        assert_eq!(selected_point.to_extensible(), expected_point);
 
-        expected_point = expected_point
-            .to_extensible()
-            .add_extended(&p)
-            .to_extended();
+        expected_point = expected_point.add_extended(&p).to_extended();
     }
 }
diff --git a/ed448-goldilocks/src/curve/twedwards/affine.rs b/ed448-goldilocks/src/curve/twedwards/affine.rs
index 055e49014..578e00c2c 100644
--- a/ed448-goldilocks/src/curve/twedwards/affine.rs
+++ b/ed448-goldilocks/src/curve/twedwards/affine.rs
@@ -121,13 +121,14 @@ impl AffineNielsPoint {
             && (self.td == other.td)
     }
 
-    /// Converts an AffineNielsPoint to an ExtendedPoint
-    pub(crate) fn to_extended(self) -> ExtendedPoint {
-        ExtendedPoint {
+    /// Converts an AffineNielsPoint to an ExtensiblePoint
+    pub(crate) fn to_extensible(self) -> ExtensiblePoint {
+        ExtensiblePoint {
             X: self.y_plus_x - self.y_minus_x,
             Y: self.y_minus_x + self.y_plus_x,
             Z: FieldElement::ONE,
-            T: self.y_plus_x * self.y_minus_x,
+            T1: self.y_plus_x,
+            T2: self.y_minus_x,
         }
     }
 }
@@ -140,7 +141,7 @@ mod tests {
     #[test]
     fn test_negation() {
         use crate::TWISTED_EDWARDS_BASE_POINT;
-        let a = TWISTED_EDWARDS_BASE_POINT.to_affine();
+        let a = TWISTED_EDWARDS_BASE_POINT.to_extensible().to_affine();
         assert!(a.is_on_curve());
 
         let neg_a = a.negate();
diff --git a/ed448-goldilocks/src/curve/twedwards/extended.rs b/ed448-goldilocks/src/curve/twedwards/extended.rs
index 609c391b0..db15aa8e5 100644
--- a/ed448-goldilocks/src/curve/twedwards/extended.rs
+++ b/ed448-goldilocks/src/curve/twedwards/extended.rs
@@ -1,9 +1,9 @@
 #![allow(non_snake_case)]
 #![allow(dead_code)]
 
-use crate::curve::twedwards::affine::AffinePoint;
+use crate::curve::twedwards::affine::AffineNielsPoint;
 use crate::curve::twedwards::extensible::ExtensiblePoint;
-use crate::edwards::EdwardsPoint as EdwardsExtendedPoint;
+use crate::curve::twedwards::projective::ProjectiveNielsPoint;
 use crate::field::FieldElement;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq};
 
@@ -43,6 +43,11 @@ impl PartialEq for ExtendedPoint {
         self.ct_eq(other).into()
     }
 }
+impl PartialEq<ExtensiblePoint> for ExtendedPoint {
+    fn eq(&self, other: &ExtensiblePoint) -> bool {
+        self.to_extensible().ct_eq(other).into()
+    }
+}
 impl Eq for ExtendedPoint {}
 
 impl Default for ExtendedPoint {
@@ -69,14 +74,90 @@ impl ExtendedPoint {
         T: FieldElement::ZERO,
     };
 
-    /// Doubles an extended point
-    pub(crate) fn double(&self) -> ExtendedPoint {
-        self.to_extensible().double().to_extended()
+    /// Adds an extensible point to an extended point
+    /// Returns an extensible point
+    /// (3.1) https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
+    pub fn add_extended(&self, other: &ExtendedPoint) -> ExtensiblePoint {
+        let A = self.X * other.X;
+        let B = self.Y * other.Y;
+        let C = self.T * other.T * FieldElement::TWISTED_D;
+        let D = self.Z * other.Z;
+        let E = (self.X + self.Y) * (other.X + other.Y) - A - B;
+        let F = D - C;
+        let G = D + C;
+        let H = B + A;
+        ExtensiblePoint {
+            X: E * F,
+            Y: G * H,
+            T1: E,
+            T2: H,
+            Z: F * G,
+        }
+    }
+
+    /// Subtracts an extensible point from an extended point
+    /// Returns an extensible point
+    /// This is a direct modification of the addition formula to the negation of `other`
+    pub fn sub_extended(&self, other: &ExtendedPoint) -> ExtensiblePoint {
+        let A = self.X * other.X;
+        let B = self.Y * other.Y;
+        let C = self.T * other.T * FieldElement::TWISTED_D;
+        let D = self.Z * other.Z;
+        let E = (self.X + self.Y) * (other.Y - other.X) + A - B;
+        let F = D + C;
+        let G = D - C;
+        let H = B - A;
+        ExtensiblePoint {
+            X: E * F,
+            Y: G * H,
+            T1: E,
+            T2: H,
+            Z: F * G,
+        }
+    }
+
+    /// Adds an extensible point to an AffineNiels point
+    /// Returns an Extensible point
+    pub fn add_affine_niels(&self, other: AffineNielsPoint) -> ExtensiblePoint {
+        let A = other.y_minus_x * (self.Y - self.X);
+        let B = other.y_plus_x * (self.X + self.Y);
+        let C = other.td * self.T;
+        let D = B + A;
+        let E = B - A;
+        let F = self.Z - C;
+        let G = self.Z + C;
+        ExtensiblePoint {
+            X: E * F,
+            Y: G * D,
+            Z: F * G,
+            T1: E,
+            T2: D,
+        }
     }
 
-    /// Adds an extended point to itself
-    pub(crate) fn add(&self, other: &ExtendedPoint) -> ExtendedPoint {
-        self.to_extensible().add_extended(other).to_extended()
+    /// Adds an extensible point to a ProjectiveNiels point
+    /// Returns an extensible point
+    /// (3.1)[Last set of formulas] https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
+    /// This differs from the formula above by a factor of 2. Saving 1 Double
+    /// Cost 8M
+    pub fn add_projective_niels(&self, other: &ProjectiveNielsPoint) -> ExtensiblePoint {
+        // This is the only step which makes it different than adding an AffineNielsPoint
+        let Z = self.Z * other.Z;
+
+        let A = (self.Y - self.X) * other.Y_minus_X;
+        let B = (self.Y + self.X) * other.Y_plus_X;
+        let C = other.Td * self.T;
+        let D = B + A;
+        let E = B - A;
+        let F = Z - C;
+        let G = Z + C;
+        ExtensiblePoint {
+            X: E * F,
+            Y: G * D,
+            Z: F * G,
+            T1: E,
+            T2: D,
+        }
     }
 
     /// Converts an ExtendedPoint to an ExtensiblePoint
@@ -90,53 +171,16 @@ impl ExtendedPoint {
         }
     }
 
-    /// Converts an extended point to Affine co-ordinates
-    pub(crate) fn to_affine(self) -> AffinePoint {
-        // Points to consider:
-        // - All points where Z=0, translate to (0,0)
-        // - The identity point has z=1, so it is not a problem
-
-        let INV_Z = self.Z.invert();
-
-        let x = self.X * INV_Z;
-        let y = self.Y * INV_Z;
-
-        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 {
-        // 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();
-
-        EdwardsExtendedPoint {
-            X: new_x,
-            Y: new_y,
-            Z: FieldElement::ONE,
-            T: new_x * new_y,
+    /// Converts an Extensible point to a ProjectiveNiels Point
+    pub fn to_projective_niels(self) -> ProjectiveNielsPoint {
+        ProjectiveNielsPoint {
+            Y_plus_X: self.X + self.Y,
+            Y_minus_X: self.Y - self.X,
+            Z: self.Z.double(),
+            Td: self.T * FieldElement::TWO_TIMES_TWISTED_D,
         }
     }
 
-    /// 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;
@@ -178,6 +222,7 @@ impl ExtendedPoint {
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::curve::twedwards::affine::AffinePoint;
     use crate::{GOLDILOCKS_BASE_POINT, TWISTED_EDWARDS_BASE_POINT};
 
     fn hex_to_field(hex: &'static str) -> FieldElement {
@@ -196,7 +241,7 @@ mod tests {
         let y = hex_to_field(
             "ae05e9634ad7048db359d6205086c2b0036ed7a035884dd7b7e36d728ad8c4b80d6565833a2a3098bbbcb2bed1cda06bdaeafbcdea9386ed",
         );
-        let a = AffinePoint { x, y }.to_extended();
+        let a = AffinePoint { x, y }.to_extensible();
         let twist_a = a.to_untwisted().to_twisted();
         assert_eq!(twist_a, a.double().double())
     }
@@ -220,28 +265,28 @@ mod tests {
     #[test]
     fn test_point_add() {
         let a = TWISTED_EDWARDS_BASE_POINT;
-        let b = a.double();
+        let b = a.to_extensible().double().to_extended();
 
         // A + B = B + A = C
-        let c_1 = a.to_extensible().add_extended(&b).to_extended();
-        let c_2 = b.to_extensible().add_extended(&a).to_extended();
+        let c_1 = a.add_extended(&b).to_extended();
+        let c_2 = b.add_extended(&a).to_extended();
         assert!(c_1 == c_2);
 
         // Adding identity point should not change result
-        let c = c_1.to_extensible().add_extended(&ExtendedPoint::IDENTITY);
-        assert!(c.to_extended() == c_1);
+        let c = c_1.add_extended(&ExtendedPoint::IDENTITY);
+        assert!(c == c_1);
     }
 
     #[test]
     fn test_point_sub() {
         let a = TWISTED_EDWARDS_BASE_POINT;
-        let b = a.double();
+        let b = a.to_extensible().double().to_extended();
 
         // A - B = C
-        let c_1 = a.to_extensible().sub_extended(&b).to_extended();
+        let c_1 = a.sub_extended(&b).to_extended();
 
         // -B + A = C
-        let c_2 = b.negate().to_extensible().add_extended(&a).to_extended();
+        let c_2 = b.negate().add_extended(&a).to_extended();
         assert!(c_1 == c_2);
     }
 
@@ -250,6 +295,6 @@ mod tests {
         let a = TWISTED_EDWARDS_BASE_POINT;
         let neg_a = a.negate();
 
-        assert!(a.to_extensible().add_extended(&neg_a) == ExtensiblePoint::IDENTITY);
+        assert!(a.add_extended(&neg_a) == ExtensiblePoint::IDENTITY);
     }
 }
diff --git a/ed448-goldilocks/src/curve/twedwards/extensible.rs b/ed448-goldilocks/src/curve/twedwards/extensible.rs
index ff3dc2d07..fe2e5dc8a 100644
--- a/ed448-goldilocks/src/curve/twedwards/extensible.rs
+++ b/ed448-goldilocks/src/curve/twedwards/extensible.rs
@@ -1,9 +1,8 @@
 #![allow(non_snake_case)]
 #![allow(dead_code)]
 
-use crate::curve::twedwards::{
-    affine::AffineNielsPoint, extended::ExtendedPoint, projective::ProjectiveNielsPoint,
-};
+use crate::curve::twedwards::{affine::AffinePoint, extended::ExtendedPoint};
+use crate::edwards::EdwardsPoint as EdwardsExtendedPoint;
 use crate::field::FieldElement;
 use subtle::{Choice, ConstantTimeEq};
 
@@ -12,6 +11,7 @@ use subtle::{Choice, ConstantTimeEq};
 // Where x = X/Z , y = Y/Z , T1 * T2 = T
 // XXX: I think we have too many point representations,
 // But let's not remove any yet
+#[derive(Copy, Clone, Debug)]
 pub struct ExtensiblePoint {
     pub(crate) X: FieldElement,
     pub(crate) Y: FieldElement,
@@ -36,6 +36,11 @@ impl PartialEq for ExtensiblePoint {
         self.ct_eq(other).into()
     }
 }
+impl PartialEq<ExtendedPoint> for ExtensiblePoint {
+    fn eq(&self, other: &ExtendedPoint) -> bool {
+        self.ct_eq(&other.to_extensible()).into()
+    }
+}
 impl Eq for ExtensiblePoint {}
 
 impl ExtensiblePoint {
@@ -67,99 +72,8 @@ impl ExtensiblePoint {
         }
     }
 
-    /// Adds two extensible points together by converting the other point to a ExtendedPoint
-    pub fn add_extensible(&self, other: &ExtensiblePoint) -> ExtensiblePoint {
-        self.add_extended(&other.to_extended())
-    }
-
-    /// Adds an extensible point to an extended point
-    /// Returns an extensible point
-    /// (3.1) https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
-    pub fn add_extended(&self, other: &ExtendedPoint) -> ExtensiblePoint {
-        let A = self.X * other.X;
-        let B = self.Y * other.Y;
-        let C = self.T1 * self.T2 * other.T * FieldElement::TWISTED_D;
-        let D = self.Z * other.Z;
-        let E = (self.X + self.Y) * (other.X + other.Y) - A - B;
-        let F = D - C;
-        let G = D + C;
-        let H = B + A;
-        ExtensiblePoint {
-            X: E * F,
-            Y: G * H,
-            T1: E,
-            T2: H,
-            Z: F * G,
-        }
-    }
-
-    /// Subtracts an extensible point from an extended point
-    /// Returns an extensible point
-    /// This is a direct modification of the addition formula to the negation of `other`
-    pub fn sub_extended(&self, other: &ExtendedPoint) -> ExtensiblePoint {
-        let A = self.X * other.X;
-        let B = self.Y * other.Y;
-        let C = self.T1 * self.T2 * other.T * FieldElement::TWISTED_D;
-        let D = self.Z * other.Z;
-        let E = (self.X + self.Y) * (other.Y - other.X) + A - B;
-        let F = D + C;
-        let G = D - C;
-        let H = B - A;
-        ExtensiblePoint {
-            X: E * F,
-            Y: G * H,
-            T1: E,
-            T2: H,
-            Z: F * G,
-        }
-    }
-
-    /// Adds an extensible point to an AffineNiels point
-    /// Returns an Extensible point
-    pub fn add_affine_niels(&self, other: AffineNielsPoint) -> ExtensiblePoint {
-        let A = other.y_minus_x * (self.Y - self.X);
-        let B = other.y_plus_x * (self.X + self.Y);
-        let C = other.td * self.T1 * self.T2;
-        let D = B + A;
-        let E = B - A;
-        let F = self.Z - C;
-        let G = self.Z + C;
-        ExtensiblePoint {
-            X: E * F,
-            Y: G * D,
-            Z: F * G,
-            T1: E,
-            T2: D,
-        }
-    }
-
-    /// Adds an extensible point to a ProjectiveNiels point
-    /// Returns an extensible point
-    /// (3.1)[Last set of formulas] https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
-    /// This differs from the formula above by a factor of 2. Saving 1 Double
-    /// Cost 8M
-    pub fn add_projective_niels(&self, other: &ProjectiveNielsPoint) -> ExtensiblePoint {
-        // This is the only step which makes it different than adding an AffineNielsPoint
-        let Z = self.Z * other.Z;
-
-        let A = (self.Y - self.X) * other.Y_minus_X;
-        let B = (self.Y + self.X) * other.Y_plus_X;
-        let C = other.Td * self.T1 * self.T2;
-        let D = B + A;
-        let E = B - A;
-        let F = Z - C;
-        let G = Z + C;
-        ExtensiblePoint {
-            X: E * F,
-            Y: G * D,
-            Z: F * G,
-            T1: E,
-            T2: D,
-        }
-    }
-
     /// Converts an extensible point to an extended point
-    pub fn to_extended(&self) -> ExtendedPoint {
+    pub fn to_extended(self) -> ExtendedPoint {
         ExtendedPoint {
             X: self.X,
             Y: self.Y,
@@ -168,13 +82,50 @@ impl ExtensiblePoint {
         }
     }
 
-    /// Converts an Extensible point to a ProjectiveNiels Point
-    pub fn to_projective_niels(&self) -> ProjectiveNielsPoint {
-        ProjectiveNielsPoint {
-            Y_plus_X: self.X + self.Y,
-            Y_minus_X: self.Y - self.X,
-            Z: self.Z + self.Z,
-            Td: self.T1 * self.T2 * FieldElement::TWO_TIMES_TWISTED_D,
+    /// Converts an extended point to Affine co-ordinates
+    pub(crate) fn to_affine(self) -> AffinePoint {
+        // Points to consider:
+        // - All points where Z=0, translate to (0,0)
+        // - The identity point has z=1, so it is not a problem
+
+        let INV_Z = self.Z.invert();
+
+        let x = self.X * INV_Z;
+        let y = self.Y * INV_Z;
+
+        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 {
+        // 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();
+
+        EdwardsExtendedPoint {
+            X: new_x,
+            Y: new_y,
+            Z: FieldElement::ONE,
+            T: new_x * new_y,
         }
     }
+
+    /// Uses a 2-isogeny to map the point to the Ed448-Goldilocks
+    pub fn to_untwisted(self) -> EdwardsExtendedPoint {
+        self.edwards_isogeny(FieldElement::MINUS_ONE)
+    }
 }
diff --git a/ed448-goldilocks/src/curve/twedwards/projective.rs b/ed448-goldilocks/src/curve/twedwards/projective.rs
index d5c9cafb7..7e2ab725c 100644
--- a/ed448-goldilocks/src/curve/twedwards/projective.rs
+++ b/ed448-goldilocks/src/curve/twedwards/projective.rs
@@ -1,6 +1,7 @@
 #![allow(non_snake_case)]
 
-use crate::curve::twedwards::{extended::ExtendedPoint, extensible::ExtensiblePoint};
+use crate::curve::twedwards::extended::ExtendedPoint;
+use crate::curve::twedwards::extensible::ExtensiblePoint;
 use crate::field::FieldElement;
 use subtle::{Choice, ConditionallyNegatable, ConditionallySelectable};
 
@@ -22,7 +23,7 @@ pub struct ProjectiveNielsPoint {
 
 impl PartialEq for ProjectiveNielsPoint {
     fn eq(&self, other: &ProjectiveNielsPoint) -> bool {
-        self.to_extended().eq(&other.to_extended())
+        self.to_extensible().eq(&other.to_extensible())
     }
 }
 impl Eq for ProjectiveNielsPoint {}
@@ -46,17 +47,18 @@ impl ConditionallyNegatable for ProjectiveNielsPoint {
 
 impl ProjectiveNielsPoint {
     pub fn identity() -> ProjectiveNielsPoint {
-        ExtensiblePoint::IDENTITY.to_projective_niels()
+        ExtendedPoint::IDENTITY.to_projective_niels()
     }
 
-    pub fn to_extended(self) -> ExtendedPoint {
+    pub fn to_extensible(self) -> ExtensiblePoint {
         let A = self.Y_plus_X - self.Y_minus_X;
         let B = self.Y_plus_X + self.Y_minus_X;
-        ExtendedPoint {
+        ExtensiblePoint {
             X: self.Z * A,
             Y: self.Z * B,
             Z: self.Z.square(),
-            T: B * A,
+            T1: B,
+            T2: A,
         }
     }
 }
@@ -68,10 +70,10 @@ mod tests {
     fn test_conditional_negate() {
         let bp = ExtendedPoint::GENERATOR;
 
-        let mut bp_neg = bp.to_extensible().to_projective_niels();
+        let mut bp_neg = bp.to_projective_niels();
         bp_neg.conditional_negate(1.into());
 
-        let expect_identity = bp_neg.to_extended().add(&bp);
+        let expect_identity = bp_neg.to_extensible().to_extended().add_extended(&bp);
         assert_eq!(ExtendedPoint::IDENTITY, expect_identity);
     }
 }
diff --git a/ed448-goldilocks/src/decaf/ops.rs b/ed448-goldilocks/src/decaf/ops.rs
index 97f82f2a2..1ba168537 100644
--- a/ed448-goldilocks/src/decaf/ops.rs
+++ b/ed448-goldilocks/src/decaf/ops.rs
@@ -14,7 +14,7 @@ impl Mul<&DecafScalar> for &DecafPoint {
 
     fn mul(self, scalar: &DecafScalar) -> DecafPoint {
         // XXX: We can do better than double and add
-        DecafPoint(double_and_add(&self.0, scalar.bits()))
+        DecafPoint(double_and_add(&self.0, scalar.bits()).to_extended())
     }
 }
 
@@ -37,7 +37,7 @@ impl Add<&DecafPoint> for &DecafPoint {
     type Output = DecafPoint;
 
     fn add(self, other: &DecafPoint) -> DecafPoint {
-        DecafPoint(self.0.to_extensible().add_extended(&other.0).to_extended())
+        DecafPoint(self.0.add_extended(&other.0).to_extended())
     }
 }
 
@@ -101,7 +101,7 @@ impl Sub<&DecafPoint> for &DecafPoint {
     type Output = DecafPoint;
 
     fn sub(self, other: &DecafPoint) -> DecafPoint {
-        DecafPoint(self.0.to_extensible().sub_extended(&other.0).to_extended())
+        DecafPoint(self.0.sub_extended(&other.0).to_extended())
     }
 }
 
diff --git a/ed448-goldilocks/src/decaf/points.rs b/ed448-goldilocks/src/decaf/points.rs
index d4e1dd5c2..ae1dd8712 100644
--- a/ed448-goldilocks/src/decaf/points.rs
+++ b/ed448-goldilocks/src/decaf/points.rs
@@ -192,7 +192,7 @@ impl Group for DecafPoint {
     }
 
     fn double(&self) -> Self {
-        Self(self.0.double())
+        Self(self.0.to_extensible().double().to_extended())
     }
 }
 
@@ -240,31 +240,31 @@ impl CurveGroup for DecafPoint {
     type AffineRepr = DecafAffinePoint;
 
     fn to_affine(&self) -> Self::AffineRepr {
-        DecafAffinePoint(self.0.to_affine())
+        DecafAffinePoint(self.0.to_extensible().to_affine())
     }
 }
 
 impl From<EdwardsPoint> for DecafPoint {
     fn from(point: EdwardsPoint) -> Self {
-        Self(point.to_twisted())
+        Self(point.to_twisted().to_extended())
     }
 }
 
 impl From<&EdwardsPoint> for DecafPoint {
     fn from(point: &EdwardsPoint) -> Self {
-        Self(point.to_twisted())
+        Self(point.to_twisted().to_extended())
     }
 }
 
 impl From<DecafPoint> for EdwardsPoint {
     fn from(point: DecafPoint) -> Self {
-        point.0.to_untwisted()
+        point.0.to_extensible().to_untwisted()
     }
 }
 
 impl From<&DecafPoint> for EdwardsPoint {
     fn from(point: &DecafPoint) -> Self {
-        point.0.to_untwisted()
+        point.0.to_extensible().to_untwisted()
     }
 }
 
@@ -282,13 +282,13 @@ impl From<&DecafAffinePoint> for DecafPoint {
 
 impl From<DecafPoint> for DecafAffinePoint {
     fn from(point: DecafPoint) -> Self {
-        DecafAffinePoint(point.0.to_affine())
+        DecafAffinePoint(point.0.to_extensible().to_affine())
     }
 }
 
 impl From<&DecafPoint> for DecafAffinePoint {
     fn from(point: &DecafPoint) -> Self {
-        DecafAffinePoint(point.0.to_affine())
+        DecafAffinePoint(point.0.to_extensible().to_affine())
     }
 }
 
@@ -307,12 +307,12 @@ impl DecafPoint {
 
     /// Add two points
     pub fn add(&self, other: &DecafPoint) -> DecafPoint {
-        DecafPoint(self.0.to_extensible().add_extended(&other.0).to_extended())
+        DecafPoint(self.0.add_extended(&other.0).to_extended())
     }
 
     /// Subtract two points
     pub fn sub(&self, other: &DecafPoint) -> DecafPoint {
-        DecafPoint(self.0.to_extensible().sub_extended(&other.0).to_extended())
+        DecafPoint(self.0.sub_extended(&other.0).to_extended())
     }
 
     /// Compress this point
@@ -369,7 +369,7 @@ impl DecafPoint {
         let u1 = FieldElement::from_bytes(&hi);
         let q0 = u0.map_to_curve_decaf448();
         let q1 = u1.map_to_curve_decaf448();
-        Self(q0.add(&q1))
+        Self(q0.add_extended(&q1).to_extended())
     }
 }
 
@@ -600,8 +600,8 @@ mod test {
 
         let P = TWISTED_EDWARDS_BASE_POINT;
 
-        let P2 = P.double();
-        let P3 = P2.to_extensible().add_extended(&P).to_extended();
+        let P2 = P.to_extensible().double().to_extended();
+        let P3 = P2.add_extended(&P).to_extended();
 
         // Encode and decode to make them Decaf points
         let Decaf_P = DecafPoint(P).compress().decompress().unwrap();
diff --git a/ed448-goldilocks/src/edwards/extended.rs b/ed448-goldilocks/src/edwards/extended.rs
index c6af6a79c..2b07d5b3b 100644
--- a/ed448-goldilocks/src/edwards/extended.rs
+++ b/ed448-goldilocks/src/edwards/extended.rs
@@ -4,7 +4,7 @@ use core::iter::Sum;
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use crate::curve::scalar_mul::variable_base;
-use crate::curve::twedwards::extended::ExtendedPoint as TwistedExtendedPoint;
+use crate::curve::twedwards::extensible::ExtensiblePoint as TwistedExtensiblePoint;
 use crate::field::FieldElement;
 use crate::*;
 use elliptic_curve::{
@@ -548,7 +548,8 @@ impl EdwardsPoint {
         scalar_div_four.div_by_four();
 
         // Use isogeny and dual isogeny to compute phi^-1((s/4) * phi(P))
-        let partial_result = variable_base(&self.to_twisted(), &scalar_div_four).to_untwisted();
+        let partial_result =
+            variable_base(&self.to_twisted().to_extended(), &scalar_div_four).to_untwisted();
         // Add partial result to (scalar mod 4) * P
         partial_result.add(&self.scalar_mod_four(scalar))
     }
@@ -664,7 +665,7 @@ impl EdwardsPoint {
     /// 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 {
+    fn edwards_isogeny(&self, a: FieldElement) -> TwistedExtensiblePoint {
         // Convert to affine now, then derive extended version later
         let affine = self.to_affine();
         let x = affine.x;
@@ -681,15 +682,16 @@ impl EdwardsPoint {
         let y_denom = (FieldElement::ONE + FieldElement::ONE) - y.square() - (a * x.square());
         let new_y = y_numerator * y_denom.invert();
 
-        TwistedExtendedPoint {
+        TwistedExtensiblePoint {
             X: new_x,
             Y: new_y,
             Z: FieldElement::ONE,
-            T: new_x * new_y,
+            T1: new_x,
+            T2: new_y,
         }
     }
 
-    pub(crate) fn to_twisted(self) -> TwistedExtendedPoint {
+    pub(crate) fn to_twisted(self) -> TwistedExtensiblePoint {
         self.edwards_isogeny(FieldElement::ONE)
     }
 

From ff8aeba64315f3d42d4300a016f3921042e82a41 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Sun, 27 Jul 2025 01:09:49 +0200
Subject: [PATCH 03/10] Make `ProjectiveNielsPoint::identity()` and associated
 constant

---
 .../src/curve/scalar_mul/window/wnaf.rs       |  2 +-
 .../src/curve/twedwards/projective.rs         | 31 +++++++++++++++----
 ed448-goldilocks/src/field/element.rs         |  1 +
 3 files changed, 27 insertions(+), 7 deletions(-)

diff --git a/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs b/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs
index 73c1c651f..5f846880e 100644
--- a/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul/window/wnaf.rs
@@ -23,7 +23,7 @@ impl From<&ExtendedPoint> for LookupTable {
 impl LookupTable {
     /// Selects a projective niels point from a lookup table in constant time
     pub fn select(&self, index: u32) -> ProjectiveNielsPoint {
-        let mut result = ProjectiveNielsPoint::identity();
+        let mut result = ProjectiveNielsPoint::IDENTITY;
 
         for i in 1..9 {
             let swap = index.ct_eq(&(i as u32));
diff --git a/ed448-goldilocks/src/curve/twedwards/projective.rs b/ed448-goldilocks/src/curve/twedwards/projective.rs
index 7e2ab725c..ae0597eb6 100644
--- a/ed448-goldilocks/src/curve/twedwards/projective.rs
+++ b/ed448-goldilocks/src/curve/twedwards/projective.rs
@@ -1,19 +1,18 @@
 #![allow(non_snake_case)]
 
-use crate::curve::twedwards::extended::ExtendedPoint;
 use crate::curve::twedwards::extensible::ExtensiblePoint;
 use crate::field::FieldElement;
 use subtle::{Choice, ConditionallyNegatable, ConditionallySelectable};
 
 impl Default for ProjectiveNielsPoint {
     fn default() -> ProjectiveNielsPoint {
-        ProjectiveNielsPoint::identity()
+        ProjectiveNielsPoint::IDENTITY
     }
 }
 
 // Its a variant of Niels, where a Z coordinate is added for unmixed readdition
 // ((y+x)/2, (y-x)/2, dxy, Z)
-#[derive(Copy, Clone)]
+#[derive(Copy, Clone, Debug)]
 pub struct ProjectiveNielsPoint {
     pub(crate) Y_plus_X: FieldElement,
     pub(crate) Y_minus_X: FieldElement,
@@ -46,9 +45,12 @@ impl ConditionallyNegatable for ProjectiveNielsPoint {
 }
 
 impl ProjectiveNielsPoint {
-    pub fn identity() -> ProjectiveNielsPoint {
-        ExtendedPoint::IDENTITY.to_projective_niels()
-    }
+    pub const IDENTITY: ProjectiveNielsPoint = ProjectiveNielsPoint {
+        Y_plus_X: FieldElement::ONE,
+        Y_minus_X: FieldElement::ONE,
+        Td: FieldElement::ZERO,
+        Z: FieldElement::TWO,
+    };
 
     pub fn to_extensible(self) -> ExtensiblePoint {
         let A = self.Y_plus_X - self.Y_minus_X;
@@ -65,6 +67,23 @@ impl ProjectiveNielsPoint {
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::curve::twedwards::extended::ExtendedPoint;
+
+    #[test]
+    fn identity() {
+        // Internally are compared by converting to `ExtendedPoint`.
+        // Here the right-side identity point is converted to Niel's
+        // and then both sides are converted to twisted-curve form.
+        assert_eq!(
+            ProjectiveNielsPoint::IDENTITY,
+            ExtendedPoint::IDENTITY.to_projective_niels(),
+        );
+        // Here only the left-side identity point is converted.
+        assert_eq!(
+            ProjectiveNielsPoint::IDENTITY.to_extensible(),
+            ExtendedPoint::IDENTITY,
+        );
+    }
 
     #[test]
     fn test_conditional_negate() {
diff --git a/ed448-goldilocks/src/field/element.rs b/ed448-goldilocks/src/field/element.rs
index 2a1d9d59b..7e96a5367 100644
--- a/ed448-goldilocks/src/field/element.rs
+++ b/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",
     )));

From 46a6de6a3efcb295df5c9f611aaa7e7f091269b2 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Sun, 20 Jul 2025 00:37:54 +0200
Subject: [PATCH 04/10] Use Decaf448 specific addition algorithm

---
 ed448-goldilocks/src/curve/twedwards/extended.rs | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/ed448-goldilocks/src/curve/twedwards/extended.rs b/ed448-goldilocks/src/curve/twedwards/extended.rs
index db15aa8e5..ace31cd97 100644
--- a/ed448-goldilocks/src/curve/twedwards/extended.rs
+++ b/ed448-goldilocks/src/curve/twedwards/extended.rs
@@ -78,11 +78,11 @@ impl ExtendedPoint {
     /// Returns an extensible point
     /// (3.1) https://iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
     pub fn add_extended(&self, other: &ExtendedPoint) -> ExtensiblePoint {
-        let A = self.X * other.X;
-        let B = self.Y * other.Y;
-        let C = self.T * other.T * FieldElement::TWISTED_D;
-        let D = self.Z * other.Z;
-        let E = (self.X + self.Y) * (other.X + other.Y) - A - B;
+        let A = (self.Y - self.X) * (other.Y - other.X);
+        let B = (self.Y + self.X) * (other.Y + other.X);
+        let C = FieldElement::TWO_TIMES_TWISTED_D * self.T * other.T;
+        let D = (self.Z * other.Z).double();
+        let E = B - A;
         let F = D - C;
         let G = D + C;
         let H = B + A;

From c23ba25ec8431e29ca627f2c5d9330dfb4be9fac Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Fri, 18 Jul 2025 15:59:39 +0200
Subject: [PATCH 05/10] Reuse Edwards windowed scalar multiplication for Decaf

---
 ed448-goldilocks/src/curve/scalar_mul.rs      |  2 --
 .../src/curve/scalar_mul/double_and_add.rs    | 20 ---------------
 .../src/curve/scalar_mul/variable_base.rs     | 25 ++++++++++++++++---
 ed448-goldilocks/src/decaf/ops.rs             |  6 ++---
 4 files changed, 25 insertions(+), 28 deletions(-)
 delete mode 100644 ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs

diff --git a/ed448-goldilocks/src/curve/scalar_mul.rs b/ed448-goldilocks/src/curve/scalar_mul.rs
index b6ff71a1e..f25bfa715 100644
--- a/ed448-goldilocks/src/curve/scalar_mul.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul.rs
@@ -1,7 +1,5 @@
-pub(crate) mod double_and_add;
 // pub(crate) mod double_base;
 pub(crate) mod variable_base;
 pub(crate) mod window;
 
-pub(crate) use double_and_add::double_and_add;
 pub(crate) use variable_base::variable_base;
diff --git a/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs b/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs
deleted file mode 100644
index 917ba16d5..000000000
--- a/ed448-goldilocks/src/curve/scalar_mul/double_and_add.rs
+++ /dev/null
@@ -1,20 +0,0 @@
-use crate::curve::twedwards::extended::ExtendedPoint;
-use crate::curve::twedwards::extensible::ExtensiblePoint;
-use subtle::{Choice, ConditionallySelectable};
-
-/// Traditional double and add algorithm
-pub(crate) fn double_and_add(point: &ExtendedPoint, s_bits: [bool; 448]) -> ExtensiblePoint {
-    let mut result = ExtensiblePoint::IDENTITY;
-
-    // NB, we reverse here, so we are going from MSB to LSB
-    // XXX: Would be great if subtle had a From<u32> for Choice. But maybe that is not it's purpose?
-    for bit in s_bits.into_iter().rev() {
-        result = result.double();
-
-        let mut p = ExtendedPoint::IDENTITY;
-        p.conditional_assign(point, Choice::from(bit as u8));
-        result = result.to_extended().add_extended(&p);
-    }
-
-    result
-}
diff --git a/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs b/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
index 4c420f52f..27309deff 100644
--- a/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
@@ -1,11 +1,12 @@
 #![allow(non_snake_case)]
 
 use super::window::wnaf::LookupTable;
-use crate::EdwardsScalar;
+use crate::Scalar;
 use crate::curve::twedwards::{extended::ExtendedPoint, extensible::ExtensiblePoint};
+use crate::field::CurveWithScalar;
 use subtle::{Choice, ConditionallyNegatable};
 
-pub fn variable_base(point: &ExtendedPoint, s: &EdwardsScalar) -> ExtensiblePoint {
+pub fn variable_base<C: CurveWithScalar>(point: &ExtendedPoint, s: &Scalar<C>) -> ExtensiblePoint {
     let mut result = ExtensiblePoint::IDENTITY;
 
     // Recode Scalar
@@ -37,12 +38,30 @@ pub fn variable_base(point: &ExtendedPoint, s: &EdwardsScalar) -> ExtensiblePoin
 #[cfg(test)]
 mod test {
     use super::*;
+    use crate::EdwardsScalar;
     use crate::TWISTED_EDWARDS_BASE_POINT;
-    use crate::curve::scalar_mul::double_and_add;
     use elliptic_curve::bigint::U448;
+    use subtle::ConditionallySelectable;
 
     #[test]
     fn test_scalar_mul() {
+        /// Traditional double and add algorithm
+        fn double_and_add(point: &ExtendedPoint, s_bits: [bool; 448]) -> ExtensiblePoint {
+            let mut result = ExtensiblePoint::IDENTITY;
+
+            // NB, we reverse here, so we are going from MSB to LSB
+            // XXX: Would be great if subtle had a From<u32> for Choice. But maybe that is not it's purpose?
+            for bit in s_bits.into_iter().rev() {
+                result = result.double();
+
+                let mut p = ExtendedPoint::IDENTITY;
+                p.conditional_assign(point, Choice::from(bit as u8));
+                result = result.to_extended().add_extended(&p);
+            }
+
+            result
+        }
+
         // XXX: In the future use known multiples from Sage in bytes form?
         let twisted_point = TWISTED_EDWARDS_BASE_POINT;
         let scalar = EdwardsScalar::new(U448::from_be_hex(
diff --git a/ed448-goldilocks/src/decaf/ops.rs b/ed448-goldilocks/src/decaf/ops.rs
index 1ba168537..01d59e85b 100644
--- a/ed448-goldilocks/src/decaf/ops.rs
+++ b/ed448-goldilocks/src/decaf/ops.rs
@@ -1,4 +1,5 @@
-use crate::{DecafAffinePoint, DecafScalar, curve::scalar_mul::double_and_add};
+use crate::curve::scalar_mul::variable_base;
+use crate::{DecafAffinePoint, DecafScalar};
 use core::{
     borrow::Borrow,
     iter::Sum,
@@ -13,8 +14,7 @@ impl Mul<&DecafScalar> for &DecafPoint {
     type Output = DecafPoint;
 
     fn mul(self, scalar: &DecafScalar) -> DecafPoint {
-        // XXX: We can do better than double and add
-        DecafPoint(double_and_add(&self.0, scalar.bits()).to_extended())
+        DecafPoint(variable_base(&self.0, scalar).to_extended())
     }
 }
 

From 9732c75f731b9af28bf8406a900cd7d5e496ce4a Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Fri, 18 Jul 2025 16:20:40 +0200
Subject: [PATCH 06/10] Simplify double-and-add code

---
 .../src/curve/scalar_mul/variable_base.rs           |  8 ++++----
 ed448-goldilocks/src/curve/twedwards/extensible.rs  | 13 ++++++++++++-
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs b/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
index 27309deff..6f2c21553 100644
--- a/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
+++ b/ed448-goldilocks/src/curve/scalar_mul/variable_base.rs
@@ -53,10 +53,10 @@ mod test {
             // XXX: Would be great if subtle had a From<u32> for Choice. But maybe that is not it's purpose?
             for bit in s_bits.into_iter().rev() {
                 result = result.double();
-
-                let mut p = ExtendedPoint::IDENTITY;
-                p.conditional_assign(point, Choice::from(bit as u8));
-                result = result.to_extended().add_extended(&p);
+                result.conditional_assign(
+                    &result.to_extended().add_extended(point),
+                    Choice::from(u8::from(bit)),
+                );
             }
 
             result
diff --git a/ed448-goldilocks/src/curve/twedwards/extensible.rs b/ed448-goldilocks/src/curve/twedwards/extensible.rs
index fe2e5dc8a..dd790040b 100644
--- a/ed448-goldilocks/src/curve/twedwards/extensible.rs
+++ b/ed448-goldilocks/src/curve/twedwards/extensible.rs
@@ -4,7 +4,7 @@
 use crate::curve::twedwards::{affine::AffinePoint, extended::ExtendedPoint};
 use crate::edwards::EdwardsPoint as EdwardsExtendedPoint;
 use crate::field::FieldElement;
-use subtle::{Choice, ConstantTimeEq};
+use subtle::{Choice, ConditionallySelectable, ConstantTimeEq};
 
 /// This is the representation that we will do most of the group operations on.
 // In affine (x,y) is the extensible point (X, Y, Z, T1, T2)
@@ -20,6 +20,17 @@ pub struct ExtensiblePoint {
     pub(crate) T2: FieldElement,
 }
 
+impl ConditionallySelectable for ExtensiblePoint {
+    fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self {
+        Self {
+            X: FieldElement::conditional_select(&a.X, &b.X, choice),
+            Y: FieldElement::conditional_select(&a.Y, &b.Y, choice),
+            Z: FieldElement::conditional_select(&a.Z, &b.Z, choice),
+            T1: FieldElement::conditional_select(&a.T1, &b.T1, choice),
+            T2: FieldElement::conditional_select(&a.T2, &b.T2, choice),
+        }
+    }
+}
 impl ConstantTimeEq for ExtensiblePoint {
     fn ct_eq(&self, other: &Self) -> Choice {
         let XZ = self.X * other.Z;

From 220ea813b90958a36cd476e48e809bcd11f112cd Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Tue, 29 Jul 2025 16:34:49 +0200
Subject: [PATCH 07/10] Use `ConstMontyForm::invert()` in place of `pow()`

---
 ed448-goldilocks/src/field/element.rs | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/ed448-goldilocks/src/field/element.rs b/ed448-goldilocks/src/field/element.rs
index 7e96a5367..2c3a1632e 100644
--- a/ed448-goldilocks/src/field/element.rs
+++ b/ed448-goldilocks/src/field/element.rs
@@ -265,10 +265,7 @@ impl FieldElement {
     /// Inverts a field element
     /// Previous chain length: 462, new length 460
     pub fn invert(&self) -> Self {
-        const INV_EXP: U448 = U448::from_be_hex(
-            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffffffffffffffffffffffffffffffffffffffffffffffffffffd",
-        );
-        Self(self.0.pow(&INV_EXP))
+        Self(self.0.invert().unwrap_or(ConstMontyType::default()))
     }
 
     pub fn square(&self) -> Self {

From 42c0cd584dd7d97619ab1458bc0cb62f4afd20e3 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Tue, 29 Jul 2025 23:14:13 +0200
Subject: [PATCH 08/10] ed448-goldilocks: add basic benchmark suite

---
 Cargo.lock                        |   1 +
 ed448-goldilocks/Cargo.toml       |   5 +
 ed448-goldilocks/benches/bench.rs | 181 ++++++++++++++++++++++++++++++
 3 files changed, 187 insertions(+)
 create mode 100644 ed448-goldilocks/benches/bench.rs

diff --git a/Cargo.lock b/Cargo.lock
index c022e4b19..7665d9357 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -392,6 +392,7 @@ dependencies = [
 name = "ed448-goldilocks"
 version = "0.14.0-pre.2"
 dependencies = [
+ "criterion",
  "ed448",
  "elliptic-curve",
  "hash2curve",
diff --git a/ed448-goldilocks/Cargo.toml b/ed448-goldilocks/Cargo.toml
index 8096392ee..99be7d896 100644
--- a/ed448-goldilocks/Cargo.toml
+++ b/ed448-goldilocks/Cargo.toml
@@ -37,9 +37,14 @@ signing = ["dep:ed448", "dep:signature"]
 serde = ["dep:serdect", "ed448?/serde_bytes"]
 
 [dev-dependencies]
+criterion = { version = "0.7", default-features = false, features = ["cargo_bench_support"] }
 hex-literal = "1"
 hex = "0.4"
 rand_core = { version = "0.9", features = ["os_rng"] }
 rand_chacha = "0.9"
 serde_bare = "0.5"
 serde_json = "1.0"
+
+[[bench]]
+harness = false
+name = "bench"
diff --git a/ed448-goldilocks/benches/bench.rs b/ed448-goldilocks/benches/bench.rs
new file mode 100644
index 000000000..96c029b90
--- /dev/null
+++ b/ed448-goldilocks/benches/bench.rs
@@ -0,0 +1,181 @@
+use criterion::{BatchSize, Criterion, criterion_group, criterion_main};
+use ed448_goldilocks::{
+    CompressedDecaf, CompressedEdwardsY, Decaf448, DecafPoint, DecafScalar, EdwardsPoint,
+    EdwardsScalar, MontgomeryPoint,
+};
+use elliptic_curve::{Field, Group};
+use hash2curve::{ExpandMsgXof, GroupDigest};
+use rand_core::{OsRng, TryRngCore};
+use sha3::Shake256;
+
+pub fn ed448(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Ed448");
+
+    group.bench_function("scalar multiplication", |b| {
+        b.iter_batched(
+            || {
+                let point = EdwardsPoint::try_from_rng(&mut OsRng).unwrap();
+                let scalar = EdwardsScalar::try_from_rng(&mut OsRng).unwrap();
+                (point, scalar)
+            },
+            |(point, scalar)| point * scalar,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("point addition", |b| {
+        b.iter_batched(
+            || {
+                let p1 = EdwardsPoint::try_from_rng(&mut OsRng).unwrap();
+                let p2 = EdwardsPoint::try_from_rng(&mut OsRng).unwrap();
+                (p1, p2)
+            },
+            |(p1, p2)| p1 + p2,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("hash_to_curve", |b| {
+        b.iter_batched(
+            || {
+                let mut msg = [0; 64];
+                OsRng.try_fill_bytes(&mut msg).unwrap();
+                msg
+            },
+            |msg| EdwardsPoint::hash_with_defaults(&msg),
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("encode_to_curve", |b| {
+        b.iter_batched(
+            || {
+                let mut msg = [0; 64];
+                OsRng.try_fill_bytes(&mut msg).unwrap();
+                msg
+            },
+            |msg| EdwardsPoint::encode_with_defaults(&msg),
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("compress", |b| {
+        b.iter_batched(
+            || EdwardsPoint::try_from_rng(&mut OsRng).unwrap(),
+            |point| point.compress().0,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("decompress", |b| {
+        b.iter_batched(
+            || EdwardsPoint::try_from_rng(&mut OsRng).unwrap().compress().0,
+            |bytes| CompressedEdwardsY(bytes).decompress().unwrap(),
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.finish();
+}
+
+pub fn decaf448(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Decaf448");
+
+    group.bench_function("scalar multiplication", |b| {
+        b.iter_batched(
+            || {
+                let point = DecafPoint::try_from_rng(&mut OsRng).unwrap();
+                let scalar = DecafScalar::try_from_rng(&mut OsRng).unwrap();
+                (point, scalar)
+            },
+            |(point, scalar)| point * scalar,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("point addition", |b| {
+        b.iter_batched(
+            || {
+                let p1 = DecafPoint::try_from_rng(&mut OsRng).unwrap();
+                let p2 = DecafPoint::try_from_rng(&mut OsRng).unwrap();
+                (p1, p2)
+            },
+            |(p1, p2)| p1 + p2,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("hash_to_curve", |b| {
+        b.iter_batched(
+            || {
+                let mut msg = [0; 64];
+                OsRng.try_fill_bytes(&mut msg).unwrap();
+                msg
+            },
+            |msg| {
+                Decaf448::hash_from_bytes::<ExpandMsgXof<Shake256>>(
+                    &[&msg],
+                    &[b"decaf448_XOF:SHAKE256_D448MAP_RO_"],
+                )
+            },
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("encode_to_curve", |b| {
+        b.iter_batched(
+            || {
+                let mut msg = [0; 64];
+                OsRng.try_fill_bytes(&mut msg).unwrap();
+                msg
+            },
+            |msg| {
+                Decaf448::encode_from_bytes::<ExpandMsgXof<Shake256>>(
+                    &[&msg],
+                    &[b"decaf448_XOF:SHAKE256_D448MAP_NU_"],
+                )
+            },
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("compress", |b| {
+        b.iter_batched(
+            || DecafPoint::try_from_rng(&mut OsRng).unwrap(),
+            |point| point.compress().0,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("decompress", |b| {
+        b.iter_batched(
+            || DecafPoint::try_from_rng(&mut OsRng).unwrap().compress().0,
+            |bytes| CompressedDecaf(bytes).decompress().unwrap(),
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.finish();
+}
+
+pub fn curve448(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Curve448");
+
+    group.bench_function("scalar multiplication", |b| {
+        b.iter_batched(
+            || {
+                let mut point = MontgomeryPoint::default();
+                OsRng.try_fill_bytes(&mut point.0).unwrap();
+                let scalar = EdwardsScalar::try_from_rng(&mut OsRng).unwrap();
+                (point, scalar)
+            },
+            |(point, scalar)| &point * &scalar,
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.finish();
+}
+
+criterion_group!(benches, ed448, decaf448, curve448);
+criterion_main!(benches);

From f291a56b190d543b82aec6ff37a86fb836052f25 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Wed, 30 Jul 2025 18:26:51 +0200
Subject: [PATCH 09/10] Compute Ed448 scalar multiplication in untwisted form

---
 ed448-goldilocks/src/edwards.rs          |   2 +
 ed448-goldilocks/src/edwards/extended.rs |  37 +------
 ed448-goldilocks/src/edwards/mul.rs      | 119 +++++++++++++++++++++++
 ed448-goldilocks/src/field/scalar.rs     |   7 --
 4 files changed, 123 insertions(+), 42 deletions(-)
 create mode 100644 ed448-goldilocks/src/edwards/mul.rs

diff --git a/ed448-goldilocks/src/edwards.rs b/ed448-goldilocks/src/edwards.rs
index 4d5245ca8..59097ae56 100644
--- a/ed448-goldilocks/src/edwards.rs
+++ b/ed448-goldilocks/src/edwards.rs
@@ -11,7 +11,9 @@
 /// If this is a problem, one can use a different isogeny strategy (Decaf)
 pub(crate) mod affine;
 pub(crate) mod extended;
+mod mul;
 mod scalar;
 pub use affine::AffinePoint;
 pub use extended::{CompressedEdwardsY, EdwardsPoint};
+use mul::scalar_mul;
 pub use scalar::{EdwardsScalar, EdwardsScalarBytes, WideEdwardsScalarBytes};
diff --git a/ed448-goldilocks/src/edwards/extended.rs b/ed448-goldilocks/src/edwards/extended.rs
index 2b07d5b3b..903bed4d0 100644
--- a/ed448-goldilocks/src/edwards/extended.rs
+++ b/ed448-goldilocks/src/edwards/extended.rs
@@ -3,7 +3,7 @@ use core::fmt::{Display, Formatter, LowerHex, Result as FmtResult, UpperHex};
 use core::iter::Sum;
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-use crate::curve::scalar_mul::variable_base;
+use super::scalar_mul;
 use crate::curve::twedwards::extensible::ExtensiblePoint as TwistedExtensiblePoint;
 use crate::field::FieldElement;
 use crate::*;
@@ -543,40 +543,7 @@ impl EdwardsPoint {
 
     /// Generic scalar multiplication to compute s*P
     pub fn scalar_mul(&self, scalar: &EdwardsScalar) -> Self {
-        // Compute floor(s/4)
-        let mut scalar_div_four = *scalar;
-        scalar_div_four.div_by_four();
-
-        // Use isogeny and dual isogeny to compute phi^-1((s/4) * phi(P))
-        let partial_result =
-            variable_base(&self.to_twisted().to_extended(), &scalar_div_four).to_untwisted();
-        // Add partial result to (scalar mod 4) * P
-        partial_result.add(&self.scalar_mod_four(scalar))
-    }
-
-    /// Returns (scalar mod 4) * P in constant time
-    pub(crate) fn scalar_mod_four(&self, scalar: &EdwardsScalar) -> Self {
-        // Compute compute (scalar mod 4)
-        let s_mod_four = scalar[0] & 3;
-
-        // Compute all possible values of (scalar mod 4) * P
-        let zero_p = EdwardsPoint::IDENTITY;
-        let one_p = self;
-        let two_p = one_p.double();
-        let three_p = two_p.add(self);
-
-        // Under the reasonable assumption that `==` is constant time
-        // Then the whole function is constant time.
-        // This should be cheaper than calling double_and_add or a scalar mul operation
-        // as the number of possibilities are so small.
-        // XXX: This claim has not been tested (although it sounds intuitive to me)
-        let mut result = EdwardsPoint::IDENTITY;
-        result.conditional_assign(&zero_p, Choice::from((s_mod_four == 0) as u8));
-        result.conditional_assign(one_p, Choice::from((s_mod_four == 1) as u8));
-        result.conditional_assign(&two_p, Choice::from((s_mod_four == 2) as u8));
-        result.conditional_assign(&three_p, Choice::from((s_mod_four == 3) as u8));
-
-        result
+        scalar_mul(self, scalar)
     }
 
     /// Standard compression; store Y and sign of X
diff --git a/ed448-goldilocks/src/edwards/mul.rs b/ed448-goldilocks/src/edwards/mul.rs
new file mode 100644
index 000000000..397482b70
--- /dev/null
+++ b/ed448-goldilocks/src/edwards/mul.rs
@@ -0,0 +1,119 @@
+use super::{EdwardsPoint, EdwardsScalar};
+use crate::field::FieldElement;
+use subtle::{Choice, ConditionallyNegatable, ConditionallySelectable, ConstantTimeEq};
+
+pub(super) fn scalar_mul(point: &EdwardsPoint, scalar: &EdwardsScalar) -> EdwardsPoint {
+    let mut result = ExtensiblePoint::IDENTITY;
+
+    // Recode Scalar
+    let scalar = scalar.to_radix_16();
+
+    let lookup = LookupTable::from(point);
+
+    for i in (0..113).rev() {
+        result = result.double();
+        result = result.double();
+        result = result.double();
+        result = result.double();
+
+        // The mask is the top bit, will be 1 for negative numbers, 0 for positive numbers
+        let mask = scalar[i] >> 7;
+        let sign = mask & 0x1;
+        // Use the mask to get the absolute value of scalar
+        let abs_value = ((scalar[i] + mask) ^ mask) as u32;
+
+        let mut neg_P = lookup.select(abs_value);
+        neg_P.conditional_negate(Choice::from((sign) as u8));
+
+        result = (EdwardsPoint::from(result) + neg_P).into();
+    }
+
+    result.into()
+}
+
+struct ExtensiblePoint {
+    X: FieldElement,
+    Y: FieldElement,
+    Z: FieldElement,
+    T1: FieldElement,
+    T2: FieldElement,
+}
+
+impl ExtensiblePoint {
+    const IDENTITY: ExtensiblePoint = ExtensiblePoint {
+        X: FieldElement::ZERO,
+        Y: FieldElement::ONE,
+        Z: FieldElement::ONE,
+        T1: FieldElement::ZERO,
+        T2: FieldElement::ONE,
+    };
+
+    fn double(&self) -> Self {
+        let A = self.X.square();
+        let B = self.Y.square();
+        let C = self.Z.square().double();
+        let D = A;
+        let E = (self.X + self.Y).square() - A - B;
+        let G = D + B;
+        let F = G - C;
+        let H = D - B;
+        Self {
+            X: E * F,
+            Y: G * H,
+            Z: F * G,
+            T1: E,
+            T2: H,
+        }
+    }
+}
+
+impl From<ExtensiblePoint> for EdwardsPoint {
+    fn from(value: ExtensiblePoint) -> Self {
+        Self {
+            X: value.X,
+            Y: value.Y,
+            Z: value.Z,
+            T: value.T1 * value.T2,
+        }
+    }
+}
+
+impl From<EdwardsPoint> for ExtensiblePoint {
+    fn from(value: EdwardsPoint) -> Self {
+        Self {
+            X: value.X,
+            Y: value.Y,
+            Z: value.Z,
+            T1: value.T,
+            T2: FieldElement::ONE,
+        }
+    }
+}
+
+pub struct LookupTable([EdwardsPoint; 8]);
+
+/// Precomputes odd multiples of the point passed in
+impl From<&EdwardsPoint> for LookupTable {
+    fn from(P: &EdwardsPoint) -> LookupTable {
+        let mut table = [*P; 8];
+
+        for i in 1..8 {
+            table[i] = P + table[i - 1];
+        }
+
+        LookupTable(table)
+    }
+}
+
+impl LookupTable {
+    /// Selects a projective niels point from a lookup table in constant time
+    pub fn select(&self, index: u32) -> EdwardsPoint {
+        let mut result = EdwardsPoint::IDENTITY;
+
+        for i in 1..9 {
+            let swap = index.ct_eq(&(i as u32));
+            result.conditional_assign(&self.0[i - 1], swap);
+        }
+        result
+    }
+}
diff --git a/ed448-goldilocks/src/field/scalar.rs b/ed448-goldilocks/src/field/scalar.rs
index 1407022d8..b76d22aa5 100644
--- a/ed448-goldilocks/src/field/scalar.rs
+++ b/ed448-goldilocks/src/field/scalar.rs
@@ -666,13 +666,6 @@ impl<C: CurveWithScalar> Scalar<C> {
         self.scalar.is_zero()
     }
 
-    /// Divides a scalar by four without reducing mod p
-    /// This is used in the 2-isogeny when mapping points from Ed448-Goldilocks
-    /// to Twisted-Goldilocks
-    pub(crate) fn div_by_four(&mut self) {
-        self.scalar >>= 2;
-    }
-
     // This method was modified from Curve25519-Dalek codebase. [scalar.rs]
     // We start with 14 u32s and convert them to 56 u8s.
     // We then use the code copied from Dalek to convert the 56 u8s to radix-16 and re-center the coefficients to be between [-16,16)

From d1de36fa99970b838ef9db412c4910cbb0d885a2 Mon Sep 17 00:00:00 2001
From: daxpedda <daxpedda@gmail.com>
Date: Fri, 1 Aug 2025 12:57:18 +0200
Subject: [PATCH 10/10] Add `MixedAdditionPoint`

---
 ed448-goldilocks/src/edwards/mul.rs | 87 +++++++++++++++++++++++++++--
 1 file changed, 81 insertions(+), 6 deletions(-)

diff --git a/ed448-goldilocks/src/edwards/mul.rs b/ed448-goldilocks/src/edwards/mul.rs
index 397482b70..846e0d2df 100644
--- a/ed448-goldilocks/src/edwards/mul.rs
+++ b/ed448-goldilocks/src/edwards/mul.rs
@@ -1,5 +1,6 @@
 use super::{EdwardsPoint, EdwardsScalar};
 use crate::field::FieldElement;
+use core::ops::Add;
 use subtle::{Choice, ConditionallyNegatable, ConditionallySelectable, ConstantTimeEq};
 
 pub(super) fn scalar_mul(point: &EdwardsPoint, scalar: &EdwardsScalar) -> EdwardsPoint {
@@ -25,7 +26,7 @@ pub(super) fn scalar_mul(point: &EdwardsPoint, scalar: &EdwardsScalar) -> Edward
         let mut neg_P = lookup.select(abs_value);
         neg_P.conditional_negate(Choice::from((sign) as u8));
 
-        result = (EdwardsPoint::from(result) + neg_P).into();
+        result = &EdwardsPoint::from(result) + &neg_P;
     }
 
     result.into()
@@ -90,15 +91,89 @@ impl From<EdwardsPoint> for ExtensiblePoint {
     }
 }
 
-pub struct LookupTable([EdwardsPoint; 8]);
+#[derive(Clone, Copy)]
+struct MixedAdditionPoint {
+    X: FieldElement,
+    Y: FieldElement,
+    Z: FieldElement,
+    Td: FieldElement,
+}
+
+impl MixedAdditionPoint {
+    const IDENTITY: Self = Self {
+        X: FieldElement::ZERO,
+        Y: FieldElement::ONE,
+        Z: FieldElement::ONE,
+        Td: FieldElement::ZERO,
+    };
+}
+
+impl From<&EdwardsPoint> for MixedAdditionPoint {
+    fn from(value: &EdwardsPoint) -> Self {
+        Self {
+            X: value.X,
+            Y: value.Y,
+            Z: value.Z,
+            Td: value.T * FieldElement::EDWARDS_D,
+        }
+    }
+}
+
+impl From<EdwardsPoint> for MixedAdditionPoint {
+    fn from(value: EdwardsPoint) -> Self {
+        (&value).into()
+    }
+}
+
+impl Add<&MixedAdditionPoint> for &EdwardsPoint {
+    type Output = ExtensiblePoint;
+
+    fn add(self, rhs: &MixedAdditionPoint) -> Self::Output {
+        let A = self.X * rhs.X;
+        let B = self.Y * rhs.Y;
+        let C = self.T * rhs.Td;
+        let D = self.Z * rhs.Z;
+        let E = (self.X + self.Y) * (rhs.X + rhs.Y) - A - B;
+        let F = D - C;
+        let G = D + C;
+        let H = B - A;
+        ExtensiblePoint {
+            X: E * F,
+            Y: G * H,
+            Z: F * G,
+            T1: E,
+            T2: H,
+        }
+    }
+}
+
+impl ConditionallySelectable for MixedAdditionPoint {
+    fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self {
+        Self {
+            X: FieldElement::conditional_select(&a.X, &b.X, choice),
+            Y: FieldElement::conditional_select(&a.Y, &b.Y, choice),
+            Z: FieldElement::conditional_select(&a.Z, &b.Z, choice),
+            Td: FieldElement::conditional_select(&a.Td, &b.Td, choice),
+        }
+    }
+}
+
+impl ConditionallyNegatable for MixedAdditionPoint {
+    fn conditional_negate(&mut self, choice: Choice) {
+        self.X.conditional_negate(choice);
+        self.Td.conditional_negate(choice);
+    }
+}
+
+struct LookupTable([MixedAdditionPoint; 8]);
 
 /// Precomputes odd multiples of the point passed in
 impl From<&EdwardsPoint> for LookupTable {
     fn from(P: &EdwardsPoint) -> LookupTable {
-        let mut table = [*P; 8];
+        let mut table = [MixedAdditionPoint::from(P); 8];
 
         for i in 1..8 {
-            table[i] = P + table[i - 1];
+            table[i] = EdwardsPoint::from(P + &table[i - 1]).into();
         }
 
         LookupTable(table)
@@ -107,8 +182,8 @@ impl From<&EdwardsPoint> for LookupTable {
 
 impl LookupTable {
     /// Selects a projective niels point from a lookup table in constant time
-    pub fn select(&self, index: u32) -> EdwardsPoint {
-        let mut result = EdwardsPoint::IDENTITY;
+    fn select(&self, index: u32) -> MixedAdditionPoint {
+        let mut result = MixedAdditionPoint::IDENTITY;
 
         for i in 1..9 {
             let swap = index.ct_eq(&(i as u32));