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

        let [px, py, pz] = self.coordinates();
        let [qx, qy, _qz] = other.coordinates();
        let u = qy * pz;
        let v = qx * pz;

        if u == *py {
            if v != *px || *py == FieldElement::zero() {
                // SAFETY: The point (0, 1, 0) is defined as the point at infinity.
                return Self::new([
                    FieldElement::zero(),
                    FieldElement::one(),
                    FieldElement::zero(),
                ])
                .unwrap();
            } else {
                return self.double();
            }
        }

        let u = &u - py;
        let v = &v - px;
        let vv = &v * &v;
        let uu = &u * &u;
        let vvv = &v * &vv;
        let r = &vv * px;
        let a = &uu * pz - &vvv - &r - &r;

        let x = &v * &a;
        let y = &u * (&r - &a) - &vvv * py;
        let z = &vvv * pz;

        debug_assert_eq!(
            E::defining_equation_projective(&x, &y, &z),
            FieldElement::<E::BaseField>::zero()
        );
        // SAFETY: The values `x, y, z` are computed correctly to be on the curve.
        // The assertion above verifies that the resulting point is valid.
        let point = Self::new([x, y, z]);
        point.unwrap()
    }
More efficient than operate_with, but must ensure that other is in affine form
crates/math/src/elliptic_curve/short_weierstrass/point.rs