Add Jacobian coordinates for Short Weierstrass (#917)

* first commit, add some operations for jacobian coordinates

* save work

* save work

* implemented the Jacobian coordintes and tested it

* do as clippy says

* add tests and correct some code based on PR comments

* small refactor

* fix neutral element

* Disable failing test temporarily

* add point validation, update functions and tests

* fix operate_with for jacobian coordinates

---------

Co-authored-by: Diego K <[email protected]>

math/src/elliptic_curve/point.rs

@@ -42,7 +42,7 @@ impl<E: IsEllipticCurve> ProjectivePoint<E> {
         let [x, y, z] = self.coordinates();
         // If it's the point at infinite
         if z == &FieldElement::zero() {
-            // We make sure all the points in the infinite have the same values
+            // We make sure all the points at infinite have the same values
             return Self::new([
                 FieldElement::zero(),
                 FieldElement::one(),
@@ -63,7 +63,80 @@ impl<E: IsEllipticCurve> PartialEq for ProjectivePoint<E> {
 }
 
 impl<E: IsEllipticCurve> Eq for ProjectivePoint<E> {}
+#[derive(Debug, Clone)]
+
+pub struct JacobianPoint<E: IsEllipticCurve> {
+    pub value: [FieldElement<E::BaseField>; 3],
+}
+
+impl<E: IsEllipticCurve> JacobianPoint<E> {
+    /// Creates an elliptic curve point giving the Jacobian [x: y: z] coordinates.
+    pub const fn new(value: [FieldElement<E::BaseField>; 3]) -> Self {
+        Self { value }
+    }
 
+    /// Returns the `x` coordinate of the point.
+    pub fn x(&self) -> &FieldElement<E::BaseField> {
+        &self.value[0]
+    }
+
+    /// Returns the `y` coordinate of the point.
+    pub fn y(&self) -> &FieldElement<E::BaseField> {
+        &self.value[1]
+    }
+
+    /// Returns the `z` coordinate of the point.
+    pub fn z(&self) -> &FieldElement<E::BaseField> {
+        &self.value[2]
+    }
+
+    /// Returns a tuple [x, y, z] with the coordinates of the point.
+    pub fn coordinates(&self) -> &[FieldElement<E::BaseField>; 3] {
+        &self.value
+    }
+
+    pub fn to_affine(&self) -> Self {
+        let [x, y, z] = self.coordinates();
+        // If it's the point at infinite
+        if z == &FieldElement::zero() {
+            // We make sure all the points at infinite have the same values
+            return Self::new([
+                FieldElement::one(),
+                FieldElement::one(),
+                FieldElement::zero(),
+            ]);
+        };
+        let inv_z = z.inv().unwrap();
+        let inv_z_square = inv_z.square();
+        let inv_z_cube = &inv_z_square * &inv_z;
+        JacobianPoint::new([x * inv_z_square, y * inv_z_cube, FieldElement::one()])
+    }
+}
+
+impl<E: IsEllipticCurve> PartialEq for JacobianPoint<E> {
+    fn eq(&self, other: &Self) -> bool {
+        // In Jacobian coordinates, the equality of two points is defined as:
+        // X1 * Z2^2 == X2 * Z1^2 y Y1 * Z2^3 == Y2 * Z1^3
+
+        let [px, py, pz] = self.coordinates();
+        let [qx, qy, qz] = other.coordinates();
+
+        let zp_sq = pz.square();
+        let zq_sq = qz.square();
+
+        let zp_cu = &zp_sq * pz;
+        let zq_cu = &zq_sq * qz;
+
+        let xp_zq_sq = px * zq_sq;
+        let xq_zp_sq = qx * zp_sq;
+
+        let yp_zq_cu = py * zq_cu;
+        let yq_zp_cu = qy * zp_cu;
+
+        (xp_zq_sq == xq_zp_sq) && (yp_zq_cu == yq_zp_cu)
+    }
+}
+impl<E: IsEllipticCurve> Eq for JacobianPoint<E> {}
 #[cfg(test)]
 mod tests {
     use crate::cyclic_group::IsGroup;

math/src/elliptic_curve/short_weierstrass/curves/bls12_381/curve.rs

@@ -13,6 +13,8 @@ use crate::{
 pub const SUBGROUP_ORDER: U256 =
     U256::from_hex_unchecked("73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001");
 
+pub const CURVE_COFACTOR: U256 = U256::from_hex_unchecked("0x396c8c005555e1568c00aaab0000aaab");
+
 pub type BLS12381FieldElement = FieldElement<BLS12381PrimeField>;
 pub type BLS12381TwistCurveFieldElement = FieldElement<Degree2ExtensionField>;