fn operate_with(&self, other: &Self) -> Self {
        // One of them is the neutral element.
        if self.is_neutral_element() {
            other.clone()
        } else if other.is_neutral_element() {
            self.clone()
        } else {
            let [x1, y1, _] = self.to_affine().coordinates().clone();
            let [x2, y2, _] = other.to_affine().coordinates().clone();
            // In this case P == -Q
            if x2 == x1 && &y2 + &y1 == FieldElement::zero() {
                Self::neutral_element()
            // The points are the same P == Q
            } else if self == other {
                // P = Q = (x, y)
                // y cant be zero here because if y = 0 then
                // P = Q = (x, 0) and P = -Q, which is the
                // previous case.
                let one = FieldElement::from(1);
                let (a, b) = (E::a(), E::b());

                let x1a = &a * &x1;
                let x1_square = &x1 * &x1;
                let num = &x1_square + &x1_square + x1_square + &x1a + x1a + &one;
                let den = (&b + &b) * &y1;

                // We are using that den != 0 because b and y1 aren't zero.
                // b != 0 because the cofficient b of a montgomery elliptic curve has to be different from zero.
                // y1 != 0 because if not, it woould be the case from above: x2 = x1 and y2 + y1 = 0.
                let div = unsafe { (num / den).unwrap_unchecked() };

                let new_x = &div * &div * &b - (&x1 + x2) - a;
                let new_y = div * (x1 - &new_x) - y1;

                // SAFETY:
                // - The Montgomery addition formula guarantees a **valid** curve point.
                // - `unwrap()` is safe because the input points are **valid**.
                let point = Self::new([new_x, new_y, one]);
                point.unwrap()
            // In the rest of the cases we have x1 != x2
            } else {
                let num = &y2 - &y1;
                let den = &x2 - &x1;

                let div = unsafe { (num / den).unwrap_unchecked() };

                let new_x = &div * &div * E::b() - (&x1 + &x2) - E::a();
                let new_y = div * (x1 - &new_x) - y1;

                // SAFETY:
                // - The result of the Montgomery addition formula is **guaranteed** to be a valid point.
                // - `unwrap()` is safe because we **control** the inputs.
                let point = Self::new([new_x, new_y, FieldElement::one()]);
                point.unwrap()
            }
        }
    }
Computes the addition of `self` and `other`.

This implementation follows the addition law for Montgomery curves as described in:
**Moonmath Manual, Definition 5.2.2.1, Page 94**.

# Safety

- This function assumes that both `self` and `other` are **valid** points on the curve.
- The resulting point is **guaranteed** to be valid due to the **Montgomery curve addition formula**.
- `unwrap()` is used because the formula ensures the result remains a valid curve point.
crates/math/src/elliptic_curve/montgomery/point.rs