fn operate_with(&self, other: &Self) -> Self {
        if self.is_neutral_element() {
            return other.clone();
        }

        if other.is_neutral_element() {
            return self.clone();
        }

        let [x1, y1, z1] = self.coordinates();
        let [x2, y2, z2] = other.coordinates();

        let z1_sq = z1.square(); // Z1^2
        let z2_sq = z2.square(); // Z2^2

        let u1 = x1 * &z2_sq; // U1 = X1 * Z2^2
        let u2 = x2 * &z1_sq; // U2 = X2 * Z1^2

        let z1_cu = z1 * &z1_sq; // Z1^3
        let z2_cu = z2 * &z2_sq; // Z2^3

        let s1 = y1 * &z2_cu; // S1 = Y1 * Z2^3
        let s2 = y2 * &z1_cu; // S2 = Y2 * Z1^3

        if u1 == u2 {
            if s1 == s2 {
                return self.double(); // P + P = 2P
            } else {
                return Self::neutral_element(); // P + (-P) = 0
            }
        }
        // H = U2 - U1
        let h = u2 - &u1;
        // I = (2 * H)^2
        let i = h.double().square();
        // J = H * I
        let j = -(&h * &i);

        // R = 2 * (S2 - S1)
        let r = (s2 - &s1).double();

        // V = U1 * I
        let v = u1 * &i;

        // X3 = R^2 + J - 2 * V
        let x3 = r.square() + &j - v.double();

        // Y3 = R * (V - X3) + 2 * S1 * J
        let y3 = r * (v - &x3) + (s1 * &j.double());

        // Z3 = 2 * Z1 * Z2 * H
        let z3 = z1 * z2;
        let z3 = z3.double() * h;

        debug_assert_eq!(
            E::defining_equation_jacobian(&x3, &y3, &z3),
            FieldElement::<E::BaseField>::zero()
        );
        // SAFETY: The values `x_3, y_3, z_3` are computed correctly to be on the curve.
        // The assertion above verifies that the resulting point is valid.
        Self::new([x3, y3, z3]).unwrap()
    }
Computes the addition of `self` and `other`.
https://github.com/mratsim/constantine/blob/65147ed815d96fa94a05d307c1d9980877b7d0e8/constantine/math/elliptic/ec_shortweierstrass_jacobian.md
crates/math/src/elliptic_curve/short_weierstrass/point.rs