diff --git a/crypto/src/hash/monolith/mod.rs b/crypto/src/hash/monolith/mod.rs
index 240d5418d..acf99284b 100644
--- a/crypto/src/hash/monolith/mod.rs
+++ b/crypto/src/hash/monolith/mod.rs
@@ -120,8 +120,8 @@ impl<const WIDTH: usize, const NUM_FULL_ROUNDS: usize> MonolithMersenne31<WIDTH,
     // S-box lookups
     fn bars(&self, state: &mut [u32]) {
         for state in state.iter_mut().take(NUM_BARS) {
-            *state = (self.lookup2[(*state >> 16) as u16 as usize] as u32) << 16
-                | self.lookup1[*state as u16 as usize] as u32;
+            *state = ((self.lookup2[(*state >> 16) as u16 as usize] as u32) << 16)
+                | (self.lookup1[*state as u16 as usize] as u32);
         }
     }
 
diff --git a/math/benches/criterion_elliptic_curve.rs b/math/benches/criterion_elliptic_curve.rs
index f9bffc781..85db985fc 100644
--- a/math/benches/criterion_elliptic_curve.rs
+++ b/math/benches/criterion_elliptic_curve.rs
@@ -10,6 +10,6 @@ use elliptic_curves::{
 criterion_group!(
     name = elliptic_curve_benches;
     config = Criterion::default().with_profiler(PProfProfiler::new(100, Output::Flamegraph(None)));
-    targets = bn_254_elliptic_curve_benchmarks,bls12_381_elliptic_curve_benchmarks,bls12_377_elliptic_curve_benchmarks
+    targets = bn_254_elliptic_curve_benchmarks, bls12_377_elliptic_curve_benchmarks, bls12_381_elliptic_curve_benchmarks
 );
 criterion_main!(elliptic_curve_benches);
diff --git a/math/benches/elliptic_curves/bls12_377.rs b/math/benches/elliptic_curves/bls12_377.rs
index 0708c8951..b304597c5 100644
--- a/math/benches/elliptic_curves/bls12_377.rs
+++ b/math/benches/elliptic_curves/bls12_377.rs
@@ -15,7 +15,7 @@ pub fn bls12_377_elliptic_curve_benchmarks(c: &mut Criterion) {
     let a = BLS12377Curve::generator().operate_with_self(a_val);
     let b = BLS12377Curve::generator().operate_with_self(b_val);
 
-    let mut group = c.benchmark_group("BLS12-381 Ops");
+    let mut group = c.benchmark_group("BLS12-377 Ops");
     group.significance_level(0.1).sample_size(10000);
     group.throughput(criterion::Throughput::Elements(1));
 
diff --git a/math/benches/elliptic_curves/bn_254.rs b/math/benches/elliptic_curves/bn_254.rs
index 01b39f2b2..e14ec83f9 100644
--- a/math/benches/elliptic_curves/bn_254.rs
+++ b/math/benches/elliptic_curves/bn_254.rs
@@ -31,14 +31,14 @@ pub fn bn_254_elliptic_curve_benchmarks(c: &mut Criterion) {
     let a_g1 = BN254Curve::generator().operate_with_self(a_val);
     let b_g1 = BN254Curve::generator().operate_with_self(b_val);
 
-    let a_g2 = BN254TwistCurve::generator().operate_with_self(a_val);
+    let a_g2 = BN254TwistCurve::generator().operate_with_self(b_val);
     let b_g2 = BN254TwistCurve::generator().operate_with_self(b_val);
     let f_12 = Fp12E::from_coefficients(&[
         "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
     ]);
     let f_2 = Fp2E::new([FpE::from(a_val as u64), FpE::from(b_val as u64)]);
 
-    let miller_loop_output = miller_optimized(&a_g1, &a_g2);
+    let miller_loop_output = miller_optimized(&a_g1.to_affine(), &a_g2.to_affine());
 
     let mut group = c.benchmark_group("BN254 Ops");
 
diff --git a/math/benches/fields/baby_bear.rs b/math/benches/fields/baby_bear.rs
index fbe6f8adb..48d15a438 100644
--- a/math/benches/fields/baby_bear.rs
+++ b/math/benches/fields/baby_bear.rs
@@ -15,7 +15,7 @@ use lambdaworks_math::field::{
 
 use p3_baby_bear::BabyBear;
 use p3_field::extension::BinomialExtensionField;
-use p3_field::{Field, FieldAlgebra};
+use p3_field::{Field, PrimeCharacteristicRing};
 
 use rand::random;
 use rand::Rng;
@@ -83,18 +83,25 @@ pub fn rand_babybear_u32_fp4_elements(num: usize) -> Vec<(Fp4Eu32, Fp4Eu32)> {
     }
     result
 }
-
+fn random_baby_bear<R: Rng>(rng: &mut R) -> BabyBear {
+    BabyBear::new(rng.gen::<u32>())
+}
 fn rand_babybear_elements_p3(num: usize) -> Vec<(BabyBear, BabyBear)> {
     let mut rng = rand::thread_rng();
     (0..num)
-        .map(|_| (rng.gen::<BabyBear>(), rng.gen::<BabyBear>()))
+        .map(|_| (random_baby_bear(&mut rng), random_baby_bear(&mut rng)))
         .collect()
 }
 
 fn rand_babybear_fp4_elements_p3(num: usize) -> Vec<(EF4, EF4)> {
     let mut rng = rand::thread_rng();
     (0..num)
-        .map(|_| (rng.gen::<EF4>(), rng.gen::<EF4>()))
+        .map(|_| {
+            (
+                EF4::from(random_baby_bear(&mut rng)),
+                EF4::from(random_baby_bear(&mut rng)),
+            )
+        })
         .collect()
 }
 
diff --git a/math/benches/fields/mersenne31.rs b/math/benches/fields/mersenne31.rs
index e8d99d1c2..f479382ec 100644
--- a/math/benches/fields/mersenne31.rs
+++ b/math/benches/fields/mersenne31.rs
@@ -15,7 +15,6 @@ pub type Fp2E = FieldElement<Degree2ExtensionField>;
 pub type Fp4E = FieldElement<Degree4ExtensionField>;
 
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::rand_mersenne31_field_elements"]
 pub fn rand_field_elements(num: usize) -> Vec<(F, F)> {
     let mut result = Vec::with_capacity(num);
diff --git a/math/benches/fields/mersenne31_montgomery.rs b/math/benches/fields/mersenne31_montgomery.rs
index a3298a0d1..a96ddde83 100644
--- a/math/benches/fields/mersenne31_montgomery.rs
+++ b/math/benches/fields/mersenne31_montgomery.rs
@@ -22,7 +22,6 @@ pub type F = FieldElement<Mersenne31MontgomeryPrimeField>;
 const NUM_LIMBS: usize = 1;
 
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::rand_mersenne31_mont_field_elements"]
 pub fn rand_field_elements(num: usize) -> Vec<(F, F)> {
     let mut result = Vec::with_capacity(num);
diff --git a/math/benches/fields/stark252.rs b/math/benches/fields/stark252.rs
index 39c4ef838..7a16f8ffc 100644
--- a/math/benches/fields/stark252.rs
+++ b/math/benches/fields/stark252.rs
@@ -21,7 +21,6 @@ use rand::random;
 pub type F = FieldElement<Stark252PrimeField>;
 
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::rand_field_elements"]
 pub fn rand_field_elements(num: usize) -> Vec<(F, F)> {
     let mut result = Vec::with_capacity(num);
diff --git a/math/benches/utils/stark252_utils.rs b/math/benches/utils/stark252_utils.rs
index 249da0b85..25cdd3040 100644
--- a/math/benches/utils/stark252_utils.rs
+++ b/math/benches/utils/stark252_utils.rs
@@ -14,14 +14,12 @@ pub type FE = FieldElement<F>;
 
 // NOTE: intentional duplicate to help IAI skip setup code
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::bitrev_permute"]
 pub fn bitrev_permute(input: &mut [FE]) {
     in_place_bit_reverse_permute(input);
 }
 
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::rand_field_elements"]
 pub fn rand_field_elements(order: u64) -> Vec<FE> {
     let mut result = Vec::with_capacity(1 << order);
@@ -33,14 +31,12 @@ pub fn rand_field_elements(order: u64) -> Vec<FE> {
 }
 
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::rand_poly"]
 pub fn rand_poly(order: u64) -> Polynomial<FE> {
     Polynomial::new(&rand_field_elements(order))
 }
 
 #[inline(never)]
-#[no_mangle]
 #[export_name = "util::get_twiddles"]
 pub fn twiddles(order: u64, config: RootsConfig) -> Vec<FE> {
     get_twiddles(order, config).unwrap()
diff --git a/math/src/elliptic_curve/edwards/curves/bandersnatch/curve.rs b/math/src/elliptic_curve/edwards/curves/bandersnatch/curve.rs
index 85211713a..a8c413399 100644
--- a/math/src/elliptic_curve/edwards/curves/bandersnatch/curve.rs
+++ b/math/src/elliptic_curve/edwards/curves/bandersnatch/curve.rs
@@ -12,10 +12,26 @@ impl IsEllipticCurve for BandersnatchCurve {
     type BaseField = BaseBandersnatchFieldElement;
     type PointRepresentation = EdwardsProjectivePoint<Self>;
 
-    // Values are from https://github.com/arkworks-rs/curves/blob/5a41d7f27a703a7ea9c48512a4148443ec6c747e/ed_on_bls12_381_bandersnatch/src/curves/mod.rs#L120
-    // Converted to Hex
+    /// Returns the generator point of the Bandersnatch curve.
+    ///
+    /// The generator point is defined with coordinates `(x, y, 1)`, where `x` and `y`
+    /// are precomputed constants that belong to the curve.
+    ///
+    /// # Safety
+    ///
+    /// - The generator values are taken from the [Arkworks implementation](https://github.com/arkworks-rs/curves/blob/5a41d7f27a703a7ea9c48512a4148443ec6c747e/ed_on_bls12_381_bandersnatch/src/curves/mod.rs#L120)
+    ///   and have been converted to hexadecimal.
+    /// - `unwrap()` does not panic because:
+    ///   - The generator point is **known to be valid** on the curve.
+    ///   - The function only uses **hardcoded** and **verified** constants.
+    /// - This function should **never** be modified unless the new generator is fully verified.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point coordinates (x, y) are taken from a well-tested,
+        //   verified implementation.
+        // - The constructor will only fail if the values are invalid, which is
+        //   impossible given that they are constants taken from a trusted source.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::new_base(
                 "29C132CC2C0B34C5743711777BBE42F32B79C022AD998465E1E71866A252AE18",
             ),
@@ -23,7 +39,8 @@ impl IsEllipticCurve for BandersnatchCurve {
                 "2A6C669EDA123E0F157D8B50BADCD586358CAD81EEE464605E3167B6CC974166",
             ),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/edwards/curves/ed448_goldilocks.rs b/math/src/elliptic_curve/edwards/curves/ed448_goldilocks.rs
index e9fe3d32e..0b6bd046c 100644
--- a/math/src/elliptic_curve/edwards/curves/ed448_goldilocks.rs
+++ b/math/src/elliptic_curve/edwards/curves/ed448_goldilocks.rs
@@ -13,13 +13,27 @@ impl IsEllipticCurve for Ed448Goldilocks {
     type BaseField = P448GoldilocksPrimeField;
     type PointRepresentation = EdwardsProjectivePoint<Self>;
 
-    /// Taken from https://www.rfc-editor.org/rfc/rfc7748#page-6
+    /// Returns the generator point of the Ed448-Goldilocks curve.
+    ///
+    /// This generator is taken from [RFC 7748](https://www.rfc-editor.org/rfc/rfc7748#page-6).
+    ///
+    /// # Safety
+    ///
+    /// - The generator coordinates `(x, y, 1)` are well-known, predefined constants.
+    /// - `unwrap()` is used because the values are **known to be valid** points
+    ///   on the Ed448-Goldilocks curve.
+    /// - This function must **not** be modified unless new constants are mathematically verified.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - These values are taken from RFC 7748 and are known to be valid.
+        // - `unwrap()` is safe because `new()` will only fail if the point is
+        //   invalid, which is **not possible** with hardcoded, verified values.
+        let point= Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::from_hex("4f1970c66bed0ded221d15a622bf36da9e146570470f1767ea6de324a3d3a46412ae1af72ab66511433b80e18b00938e2626a82bc70cc05e").unwrap(),
             FieldElement::<Self::BaseField>::from_hex("693f46716eb6bc248876203756c9c7624bea73736ca3984087789c1e05a0c2d73ad3ff1ce67c39c4fdbd132c4ed7c8ad9808795bf230fa14").unwrap(),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/edwards/curves/tiny_jub_jub.rs b/math/src/elliptic_curve/edwards/curves/tiny_jub_jub.rs
index 0fccdda53..df4106051 100644
--- a/math/src/elliptic_curve/edwards/curves/tiny_jub_jub.rs
+++ b/math/src/elliptic_curve/edwards/curves/tiny_jub_jub.rs
@@ -14,12 +14,26 @@ impl IsEllipticCurve for TinyJubJubEdwards {
     type BaseField = U64PrimeField<13>;
     type PointRepresentation = EdwardsProjectivePoint<Self>;
 
+    /// Returns the generator point of the TinyJubJub Edwards curve.
+    ///
+    /// This generator is taken from **Moonmath Manual (page 97)**.
+    ///
+    /// # Safety
+    ///
+    /// - The generator coordinates `(8, 5, 1)` are **predefined** and belong to the TinyJubJub curve.
+    /// - `unwrap()` is used because the generator is a **verified valid point**,
+    ///   meaning there is **no risk** of runtime failure.
+    /// - This function must **not** be modified unless the new generator is mathematically verified.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point `(8, 5, 1)` is **mathematically valid** on the curve.
+        // - `unwrap()` is safe because we **know** the point satisfies the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::from(8),
             FieldElement::from(5),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/edwards/point.rs b/math/src/elliptic_curve/edwards/point.rs
index 57985312e..2439b7216 100644
--- a/math/src/elliptic_curve/edwards/point.rs
+++ b/math/src/elliptic_curve/edwards/point.rs
@@ -12,10 +12,27 @@ use super::traits::IsEdwards;
 #[derive(Clone, Debug)]
 pub struct EdwardsProjectivePoint<E: IsEllipticCurve>(ProjectivePoint<E>);
 
-impl<E: IsEllipticCurve> EdwardsProjectivePoint<E> {
+impl<E: IsEllipticCurve + IsEdwards> EdwardsProjectivePoint<E> {
     /// Creates an elliptic curve point giving the projective [x: y: z] coordinates.
-    pub fn new(value: [FieldElement<E::BaseField>; 3]) -> Self {
-        Self(ProjectivePoint::new(value))
+    pub fn new(value: [FieldElement<E::BaseField>; 3]) -> Result<Self, EllipticCurveError> {
+        let (x, y, z) = (&value[0], &value[1], &value[2]);
+
+        // The point at infinity is (0, 1, 1).
+        // We convert every (0, y, y) into the infinity.
+        if x == &FieldElement::<E::BaseField>::zero() && z == y {
+            return Ok(Self(ProjectivePoint::new([
+                FieldElement::<E::BaseField>::zero(),
+                FieldElement::<E::BaseField>::one(),
+                FieldElement::<E::BaseField>::one(),
+            ])));
+        }
+        if z != &FieldElement::<E::BaseField>::zero()
+            && E::defining_equation_projective(x, y, z) == FieldElement::<E::BaseField>::zero()
+        {
+            Ok(Self(ProjectivePoint::new(value)))
+        } else {
+            Err(EllipticCurveError::InvalidPoint)
+        }
     }
 
     /// Returns the `x` coordinate of the point.
@@ -56,26 +73,33 @@ impl<E: IsEdwards> FromAffine<E::BaseField> for EdwardsProjectivePoint<E> {
     fn from_affine(
         x: FieldElement<E::BaseField>,
         y: FieldElement<E::BaseField>,
-    ) -> Result<Self, crate::elliptic_curve::traits::EllipticCurveError> {
-        if E::defining_equation(&x, &y) != FieldElement::zero() {
-            Err(EllipticCurveError::InvalidPoint)
-        } else {
-            let coordinates = [x, y, FieldElement::one()];
-            Ok(EdwardsProjectivePoint::new(coordinates))
-        }
+    ) -> Result<Self, EllipticCurveError> {
+        let coordinates = [x, y, FieldElement::one()];
+        EdwardsProjectivePoint::new(coordinates)
     }
 }
 
 impl<E: IsEllipticCurve> Eq for EdwardsProjectivePoint<E> {}
 
 impl<E: IsEdwards> IsGroup for EdwardsProjectivePoint<E> {
-    /// The point at infinity.
+    /// Returns the point at infinity (neutral element) in projective coordinates.
+    ///
+    /// # Safety
+    ///
+    /// - The values `[0, 1, 1]` are the **canonical representation** of the neutral element
+    ///   in the Edwards curve, meaning they are guaranteed to be a valid point.
+    /// - `unwrap()` is used because this point is **known** to be valid, so
+    ///   there is no need for additional runtime checks.
     fn neutral_element() -> Self {
-        Self::new([
+        // SAFETY:
+        // - `[0, 1, 1]` is a mathematically verified neutral element in Edwards curves.
+        // - `unwrap()` is safe because this point is **always valid**.
+        let point = Self::new([
             FieldElement::zero(),
             FieldElement::one(),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 
     fn is_neutral_element(&self) -> bool {
@@ -83,8 +107,16 @@ impl<E: IsEdwards> IsGroup for EdwardsProjectivePoint<E> {
         px == &FieldElement::zero() && py == pz
     }
 
-    /// Computes the addition of `self` and `other`.
-    /// Taken from "Moonmath" (Eq 5.38, page 97)
+    /// Computes the addition of `self` and `other` using the Edwards curve addition formula.
+    ///
+    /// This implementation follows Equation (5.38) from "Moonmath" (page 97).
+    ///
+    /// # Safety
+    ///
+    /// - The function assumes both `self` and `other` are valid points on the curve.
+    /// - The resulting coordinates are computed using a well-defined formula that
+    ///   maintains the elliptic curve invariants.
+    /// - `unwrap()` is safe because the formula guarantees the result is valid.
     fn operate_with(&self, other: &Self) -> Self {
         // This avoids dropping, which in turn saves us from having to clone the coordinates.
         let (s_affine, o_affine) = (self.to_affine(), other.to_affine());
@@ -103,13 +135,24 @@ impl<E: IsEdwards> IsGroup for EdwardsProjectivePoint<E> {
         let num_s2 = &y1y2 - E::a() * &x1x2;
         let den_s2 = &one - &dx1x2y1y2;
 
-        Self::new([&num_s1 / &den_s1, &num_s2 / &den_s2, one])
+        // SAFETY: The creation of the result point is safe because the inputs are always points that belong to the curve.
+        let point = Self::new([&num_s1 / &den_s1, &num_s2 / &den_s2, one]);
+        point.unwrap()
     }
 
     /// Returns the additive inverse of the projective point `p`
+    ///  
+    /// # Safety
+    ///
+    /// - Negating the x-coordinate of a valid Edwards point results in another valid point.
+    /// - `unwrap()` is safe because negation does not break the curve equation.
     fn neg(&self) -> Self {
         let [px, py, pz] = self.coordinates();
-        Self::new([-px, py.clone(), pz.clone()])
+        // SAFETY:
+        // - The negation formula for Edwards curves is well-defined.
+        // - The result remains a valid curve point.
+        let point = Self::new([-px, py.clone(), pz.clone()]);
+        point.unwrap()
     }
 }
 
@@ -156,17 +199,20 @@ mod tests {
             FieldElement::from(5),
             FieldElement::from(5),
             FieldElement::from(1),
-        ]);
+        ])
+        .unwrap();
         let q = EdwardsProjectivePoint::<TinyJubJubEdwards>::new([
             FieldElement::from(8),
             FieldElement::from(5),
             FieldElement::from(1),
-        ]);
+        ])
+        .unwrap();
         let expected = EdwardsProjectivePoint::<TinyJubJubEdwards>::new([
             FieldElement::from(0),
             FieldElement::from(1),
             FieldElement::from(1),
-        ]);
+        ])
+        .unwrap();
         assert_eq!(p.operate_with(&q), expected);
     }
 
diff --git a/math/src/elliptic_curve/edwards/traits.rs b/math/src/elliptic_curve/edwards/traits.rs
index 9890c051c..6c7152152 100644
--- a/math/src/elliptic_curve/edwards/traits.rs
+++ b/math/src/elliptic_curve/edwards/traits.rs
@@ -7,6 +7,8 @@ pub trait IsEdwards: IsEllipticCurve + Clone + Debug {
 
     fn d() -> FieldElement<Self::BaseField>;
 
+    // Edwards equation in affine coordinates:
+    // ax^2 + y^2 - 1 = d * x^2 * y^2
     fn defining_equation(
         x: &FieldElement<Self::BaseField>,
         y: &FieldElement<Self::BaseField>,
@@ -15,4 +17,16 @@ pub trait IsEdwards: IsEllipticCurve + Clone + Debug {
             - FieldElement::<Self::BaseField>::one()
             - Self::d() * x.pow(2_u16) * y.pow(2_u16)
     }
+
+    // Edwards equation in projective coordinates.
+    // a * x^2 * z^2 + y^2 * z^2 - z^4 = d * x^2 * y^2
+    fn defining_equation_projective(
+        x: &FieldElement<Self::BaseField>,
+        y: &FieldElement<Self::BaseField>,
+        z: &FieldElement<Self::BaseField>,
+    ) -> FieldElement<Self::BaseField> {
+        Self::a() * x.square() * z.square() + y.square() * z.square()
+            - z.square().square()
+            - Self::d() * x.square() * y.square()
+    }
 }
diff --git a/math/src/elliptic_curve/montgomery/curves/tiny_jub_jub.rs b/math/src/elliptic_curve/montgomery/curves/tiny_jub_jub.rs
index c059fee8c..9b130e703 100644
--- a/math/src/elliptic_curve/montgomery/curves/tiny_jub_jub.rs
+++ b/math/src/elliptic_curve/montgomery/curves/tiny_jub_jub.rs
@@ -14,12 +14,27 @@ impl IsEllipticCurve for TinyJubJubMontgomery {
     type BaseField = U64PrimeField<13>;
     type PointRepresentation = MontgomeryProjectivePoint<Self>;
 
+    /// Returns the generator point of the TinyJubJub Montgomery curve.
+    ///
+    /// This generator is taken from **Moonmath Manual (page 91)**.
+    ///
+    /// # Safety
+    ///
+    /// - The generator coordinates `(3, 5, 1)` are **predefined** and are **valid** points
+    ///   on the TinyJubJub Montgomery curve.
+    /// - `unwrap()` is used because the generator is **guaranteed** to satisfy
+    ///   the Montgomery curve equation.
+    /// - This function must **not** be modified unless the new generator is mathematically verified.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point `(3, 5, 1)` is **mathematically verified**.
+        // - `unwrap()` is safe because the input values **guarantee** validity.
+        let point = Self::PointRepresentation::new([
             FieldElement::from(3),
             FieldElement::from(5),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/montgomery/point.rs b/math/src/elliptic_curve/montgomery/point.rs
index f7da2e186..7e860ccdf 100644
--- a/math/src/elliptic_curve/montgomery/point.rs
+++ b/math/src/elliptic_curve/montgomery/point.rs
@@ -12,10 +12,28 @@ use super::traits::IsMontgomery;
 #[derive(Clone, Debug)]
 pub struct MontgomeryProjectivePoint<E: IsEllipticCurve>(ProjectivePoint<E>);
 
-impl<E: IsEllipticCurve> MontgomeryProjectivePoint<E> {
+impl<E: IsEllipticCurve + IsMontgomery> MontgomeryProjectivePoint<E> {
     /// Creates an elliptic curve point giving the projective [x: y: z] coordinates.
-    pub fn new(value: [FieldElement<E::BaseField>; 3]) -> Self {
-        Self(ProjectivePoint::new(value))
+    pub fn new(value: [FieldElement<E::BaseField>; 3]) -> Result<Self, EllipticCurveError> {
+        let (x, y, z) = (&value[0], &value[1], &value[2]);
+
+        if z != &FieldElement::<E::BaseField>::zero()
+            && E::defining_equation_projective(x, y, z) == FieldElement::<E::BaseField>::zero()
+        {
+            Ok(Self(ProjectivePoint::new(value)))
+        // The point at infinity is (0, 1, 0)
+        // We convert every (0, _, 0) into the infinity.
+        } else if x == &FieldElement::<E::BaseField>::zero()
+            && z == &FieldElement::<E::BaseField>::zero()
+        {
+            Ok(Self(ProjectivePoint::new([
+                FieldElement::<E::BaseField>::zero(),
+                FieldElement::<E::BaseField>::one(),
+                FieldElement::<E::BaseField>::zero(),
+            ])))
+        } else {
+            Err(EllipticCurveError::InvalidPoint)
+        }
     }
 
     /// Returns the `x` coordinate of the point.
@@ -56,13 +74,9 @@ impl<E: IsMontgomery> FromAffine<E::BaseField> for MontgomeryProjectivePoint<E>
     fn from_affine(
         x: FieldElement<E::BaseField>,
         y: FieldElement<E::BaseField>,
-    ) -> Result<Self, crate::elliptic_curve::traits::EllipticCurveError> {
-        if E::defining_equation(&x, &y) != FieldElement::zero() {
-            Err(EllipticCurveError::InvalidPoint)
-        } else {
-            let coordinates = [x, y, FieldElement::one()];
-            Ok(MontgomeryProjectivePoint::new(coordinates))
-        }
+    ) -> Result<Self, EllipticCurveError> {
+        let coordinates = [x, y, FieldElement::one()];
+        MontgomeryProjectivePoint::new(coordinates)
     }
 }
 
@@ -70,12 +84,22 @@ impl<E: IsEllipticCurve> Eq for MontgomeryProjectivePoint<E> {}
 
 impl<E: IsMontgomery> IsGroup for MontgomeryProjectivePoint<E> {
     /// The point at infinity.
+    ///    
+    /// # Safety
+    ///
+    /// - The point `(0, 1, 0)` is a well-defined **neutral element** for Montgomery curves.
+    /// - `unwrap_unchecked()` is used because this point is **always valid**.
     fn neutral_element() -> Self {
-        Self::new([
+        // SAFETY:
+        // - `(0, 1, 0)` is **mathematically valid** as the neutral element.
+        // - `unwrap_unchecked()` is safe because this is **a known valid point**.
+        let point = Self::new([
             FieldElement::zero(),
             FieldElement::one(),
             FieldElement::zero(),
-        ])
+        ]);
+        debug_assert!(point.is_ok());
+        point.unwrap()
     }
 
     fn is_neutral_element(&self) -> bool {
@@ -84,7 +108,15 @@ impl<E: IsMontgomery> IsGroup for MontgomeryProjectivePoint<E> {
     }
 
     /// Computes the addition of `self` and `other`.
-    /// Taken from "Moonmath" (Definition 5.2.2.1, page 94)
+    ///
+    /// 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.
     fn operate_with(&self, other: &Self) -> Self {
         // One of them is the neutral element.
         if self.is_neutral_element() {
@@ -115,7 +147,11 @@ impl<E: IsMontgomery> IsGroup for MontgomeryProjectivePoint<E> {
                 let new_x = &div * &div * &b - (&x1 + x2) - a;
                 let new_y = div * (x1 - &new_x) - y1;
 
-                Self::new([new_x, new_y, one])
+                // 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;
@@ -125,15 +161,29 @@ impl<E: IsMontgomery> IsGroup for MontgomeryProjectivePoint<E> {
                 let new_x = &div * &div * E::b() - (&x1 + &x2) - E::a();
                 let new_y = div * (x1 - &new_x) - y1;
 
-                Self::new([new_x, new_y, FieldElement::one()])
+                // 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()
             }
         }
     }
 
     /// Returns the additive inverse of the projective point `p`
+    ///
+    /// # Safety
+    ///
+    /// - The negation formula preserves the curve equation.
+    /// - `unwrap()` is safe because negation **does not** create invalid points.
     fn neg(&self) -> Self {
         let [px, py, pz] = self.coordinates();
-        Self::new([px.clone(), -py, pz.clone()])
+        // SAFETY:
+        // - Negating `y` maintains the curve structure.
+        // - `unwrap()` is safe because negation **is always valid**.
+        let point = Self::new([px.clone(), -py, pz.clone()]);
+        debug_assert!(point.is_ok());
+        point.unwrap()
     }
 }
 
@@ -182,17 +232,20 @@ mod tests {
             FieldElement::from(9),
             FieldElement::from(2),
             FieldElement::from(1),
-        ]);
+        ])
+        .unwrap();
         let q = MontgomeryProjectivePoint::<TinyJubJubMontgomery>::new([
             FieldElement::from(7),
             FieldElement::from(12),
             FieldElement::from(1),
-        ]);
+        ])
+        .unwrap();
         let expected = MontgomeryProjectivePoint::<TinyJubJubMontgomery>::new([
             FieldElement::from(10),
             FieldElement::from(3),
             FieldElement::from(1),
-        ]);
+        ])
+        .unwrap();
         assert_eq!(p.operate_with(&q), expected);
     }
 
diff --git a/math/src/elliptic_curve/montgomery/traits.rs b/math/src/elliptic_curve/montgomery/traits.rs
index eda4b318e..9f8181e87 100644
--- a/math/src/elliptic_curve/montgomery/traits.rs
+++ b/math/src/elliptic_curve/montgomery/traits.rs
@@ -8,12 +8,23 @@ pub trait IsMontgomery: IsEllipticCurve + Clone + Debug {
 
     fn b() -> FieldElement<Self::BaseField>;
 
-    /// Evaluates the short Weierstrass equation at (x, y z).
-    /// Used for checking if [x: y: z] belongs to the elliptic curve.
+    /// Evaluates the equation at (x, y).
+    /// Used for checking if the point belongs to the elliptic curve.
+    /// Equation: by^2 = x^3 + ax^2 + x.
     fn defining_equation(
         x: &FieldElement<Self::BaseField>,
         y: &FieldElement<Self::BaseField>,
     ) -> FieldElement<Self::BaseField> {
-        (Self::b() * y.pow(2_u16)) - (x.pow(3_u16) + Self::a() * x.pow(2_u16) + x)
+        (Self::b() * y.square()) - (x.pow(3_u16) + Self::a() * x.square() + x)
+    }
+
+    /// Evaluates the equation at the projective point (x, y, z).
+    /// Projective equation: zby^2 = x^3 + zax^2 + z^2x
+    fn defining_equation_projective(
+        x: &FieldElement<Self::BaseField>,
+        y: &FieldElement<Self::BaseField>,
+        z: &FieldElement<Self::BaseField>,
+    ) -> FieldElement<Self::BaseField> {
+        z * Self::b() * y.square() - x.pow(3_u16) - z * Self::a() * x.square() - z.square() * x
     }
 }
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs
index 7197e8c05..7a11eeffa 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs
@@ -28,13 +28,25 @@ impl IsEllipticCurve for BLS12377Curve {
     type BaseField = BLS12377PrimeField;
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
-    // generator values are taken from https://neuromancer.sk/std/bls/BLS12-377
+    /// Returns the generator point of the BLS12-377 curve.
+    ///
+    /// Generator values are taken from [Neuromancer's BLS12-377 page](https://neuromancer.sk/std/bls/BLS12-377).
+    ///
+    /// # Safety
+    ///
+    /// - The generator point `(x, y, 1)` is predefined and is **known to be a valid point** on the curve.
+    /// - `unwrap` is used because this point is **mathematically verified**.
+    /// - Do **not** modify this function unless a new generator has been **mathematically verified**.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - These values are mathematically verified and known to be valid points on BLS12-377.
+        // - `unwrap()` is safe because we **ensure** the input values satisfy the curve equation.
+        let point= Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::new_base("8848defe740a67c8fc6225bf87ff5485951e2caa9d41bb188282c8bd37cb5cd5481512ffcd394eeab9b16eb21be9ef"),
             FieldElement::<Self::BaseField>::new_base("1914a69c5102eff1f674f5d30afeec4bd7fb348ca3e52d96d182ad44fb82305c2fe3d3634a9591afd82de55559c8ea6"),
             FieldElement::one()
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
@@ -93,17 +105,30 @@ impl ShortWeierstrassProjectivePoint<BLS12377Curve> {
 }
 
 impl ShortWeierstrassProjectivePoint<BLS12377TwistCurve> {
-    /// 𝜓(P) = 𝜁 ∘ 𝜋ₚ ∘ 𝜁⁻¹, where 𝜁 is the isomorphism u:E'(𝔽ₚ₆) −> E(𝔽ₚ₁₂) from the twist to E,, 𝜋ₚ is the p-power frobenius endomorphism
+    /// Computes 𝜓(P) = 𝜁 ∘ 𝜋ₚ ∘ 𝜁⁻¹, where 𝜁 is the isomorphism u:E'(𝔽ₚ₆) −> E(𝔽ₚ₁₂) from the twist to E,, 𝜋ₚ is the p-power frobenius endomorphism
     /// and 𝜓 satisifies minmal equation 𝑋² + 𝑡𝑋 + 𝑞 = 𝑂
     /// https://eprint.iacr.org/2022/352.pdf 4.2 (7)
     /// ψ(P) = (ψ_x * conjugate(x), ψ_y * conjugate(y), conjugate(z))
+    ///
+    /// # Safety
+    ///
+    /// - This function assumes `self` is a valid point on the BLS12-377 **twist** curve.
+    /// - The conjugation operation preserves validity.
+    /// - `unwrap()` is used because `psi()` is defined to **always return a valid point**.
     fn psi(&self) -> Self {
         let [x, y, z] = self.coordinates();
-        Self::new([
+        // SAFETY:
+        // - `conjugate()` preserves the validity of the field element.
+        // - `ENDO_U` and `ENDO_V` are precomputed constants that ensure the
+        //   resulting point satisfies the curve equation.
+        // - `unwrap()` is safe because the transformation follows
+        //   **a known valid isomorphism** between the twist and E.
+        let point = Self::new([
             x.conjugate() * GAMMA_12,
             y.conjugate() * GAMMA_13,
             z.conjugate(),
-        ])
+        ]);
+        point.unwrap()
     }
 
     /// 𝜓(P) = 𝑢P, where 𝑢 = SEED of the curve
@@ -186,7 +211,7 @@ mod tests {
         let x = FpE::new_base("134e4cc122cb62a06767fb98e86f2d5f77e2a12fefe23bb0c4c31d1bd5348b88d6f5e5dee2b54db4a2146cc9f249eea") * FpE::from(2);
         let y = FpE::new_base("17949c29effee7a9f13f69b1c28eccd78c1ed12b47068836473481ff818856594fd9c1935e3d9e621901a2d500257a2") * FpE::from(2);
         let z = FpE::from(2);
-        let point_2 = ShortWeierstrassProjectivePoint::<BLS12377Curve>::new([x, y, z]);
+        let point_2 = ShortWeierstrassProjectivePoint::<BLS12377Curve>::new([x, y, z]).unwrap();
 
         let first_algorithm_result = point_2.operate_with(&point_1).to_affine();
         let second_algorithm_result = point_2.operate_with_affine(&point_1).to_affine();
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/pairing.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/pairing.rs
index ff61e4edd..c17d654c5 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/pairing.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/pairing.rs
@@ -160,13 +160,19 @@ fn line(p: &G1Point, t: &G2Point, q: &G2Point) -> (G2Point, Fp12E) {
         let y_r = g.square() - (e_square.double() + e_square);
         let z_r = b * &h;
 
+        debug_assert_eq!(
+            BLS12377TwistCurve::defining_equation_projective(&x_r, &y_r, &z_r),
+            Fp2E::zero()
+        );
+        // SAFETY: `unwrap()` is used here because we ensure that `x_r, y_r, z_r`
+        // satisfy the curve equation. The previous assertion checks that this is indeed the case.
         let r = G2Point::new([x_r, y_r, z_r]);
 
         let l = Fp12E::new([
             Fp6E::new([y_p * (-h), Fp2E::zero(), Fp2E::zero()]),
             Fp6E::new([x_p * (j.double() + &j), i, Fp2E::zero()]),
         ]);
-        (r, l)
+        (r.unwrap(), l)
     } else {
         let [x_q, y_q, _] = q.coordinates();
         let [x_t, y_t, z_t] = t.coordinates();
@@ -188,13 +194,19 @@ fn line(p: &G1Point, t: &G2Point, q: &G2Point) -> (G2Point, Fp12E) {
         let y_r = &theta * (g - h) - i;
         let z_r = z_t * e;
 
+        debug_assert_eq!(
+            BLS12377TwistCurve::defining_equation_projective(&x_r, &y_r, &z_r),
+            Fp2E::zero()
+        );
+        // SAFETY: The values `x_r, y_r, z_r` are computed correctly to be on the curve.
+        // The assertion above verifies that the resulting point is valid.
         let r = G2Point::new([x_r, y_r, z_r]);
 
         let l = Fp12E::new([
             Fp6E::new([y_p * lambda, Fp2E::zero(), Fp2E::zero()]),
             Fp6E::new([x_p * (-theta), j, Fp2E::zero()]),
         ]);
-        (r, l)
+        (r.unwrap(), l)
     }
 }
 
@@ -405,15 +417,19 @@ mod tests {
 
     #[test]
     fn ate_pairing_errors_when_one_element_is_not_in_subgroup() {
+        // p = (0, 1, 1) is in the curve but not in the subgroup.
+        // Recall that the BLS 12-377 curve equation is y^2 = x^3 + 1.
         let p = ShortWeierstrassProjectivePoint::new([
+            FieldElement::zero(),
             FieldElement::one(),
             FieldElement::one(),
-            FieldElement::one(),
-        ]);
+        ])
+        .unwrap();
         let q = ShortWeierstrassProjectivePoint::neutral_element();
         let result = BLS12377AtePairing::compute_batch(&[(&p.to_affine(), &q)]);
         assert!(result.is_err())
     }
+
     #[test]
     fn apply_12_times_frobenius_is_identity() {
         let f = Fp12E::from_coefficients(&[
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/twist.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/twist.rs
index e0aa87c59..ff660ec14 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/twist.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/twist.rs
@@ -20,7 +20,10 @@ impl IsEllipticCurve for BLS12377TwistCurve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.s
+        let point = Self::PointRepresentation::new([
             FieldElement::new([
                 FieldElement::new(GENERATOR_X_0),
                 FieldElement::new(GENERATOR_X_1),
@@ -30,7 +33,8 @@ impl IsEllipticCurve for BLS12377TwistCurve {
                 FieldElement::new(GENERATOR_Y_1),
             ]),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/curve.rs
index a080aa77b..67e83761b 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/curve.rs
@@ -26,12 +26,22 @@ impl IsEllipticCurve for BLS12381Curve {
     type BaseField = BLS12381PrimeField;
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
+    /// Returns the generator point of the BLS12-381 curve.
+    ///
+    /// # Safety
+    ///
+    /// - The generator point is mathematically verified to be a valid point on the curve.
+    /// - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - These values are mathematically verified and known to be valid points on BLS12-381.
+        // - `unwrap()` is safe because we **ensure** the input values satisfy the curve equation.
+        let point= Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::new_base("17f1d3a73197d7942695638c4fa9ac0fc3688c4f9774b905a14e3a3f171bac586c55e83ff97a1aeffb3af00adb22c6bb"),
             FieldElement::<Self::BaseField>::new_base("8b3f481e3aaa0f1a09e30ed741d8ae4fcf5e095d5d00af600db18cb2c04b3edd03cc744a2888ae40caa232946c5e7e1"),
             FieldElement::one()
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
@@ -89,16 +99,29 @@ impl ShortWeierstrassProjectivePoint<BLS12381Curve> {
 }
 
 impl ShortWeierstrassProjectivePoint<BLS12381TwistCurve> {
-    /// 𝜓(P) = 𝜁 ∘ 𝜋ₚ ∘ 𝜁⁻¹, where 𝜁 is the isomorphism u:E'(𝔽ₚ₆) −> E(𝔽ₚ₁₂) from the twist to E,, 𝜋ₚ is the p-power frobenius endomorphism
+    /// Computes 𝜓(P) 𝜓(P) = 𝜁 ∘ 𝜋ₚ ∘ 𝜁⁻¹, where 𝜁 is the isomorphism u:E'(𝔽ₚ₆) −> E(𝔽ₚ₁₂) from the twist to E,, 𝜋ₚ is the p-power frobenius endomorphism
     /// and 𝜓 satisifies minmal equation 𝑋² + 𝑡𝑋 + 𝑞 = 𝑂
     /// https://eprint.iacr.org/2022/352.pdf 4.2 (7)
+    ///
+    /// # Safety
+    ///
+    /// - This function assumes `self` is a valid point on the BLS12-381 **twist** curve.
+    /// - The conjugation operation preserves validity.
+    /// - `unwrap()` is used because `psi()` is defined to **always return a valid point**.
     fn psi(&self) -> Self {
         let [x, y, z] = self.coordinates();
-        Self::new([
+        // SAFETY:
+        // - `conjugate()` preserves the validity of the field element.
+        // - `ENDO_U` and `ENDO_V` are precomputed constants that ensure the
+        //   resulting point satisfies the curve equation.
+        // - `unwrap()` is safe because the transformation follows
+        //   **a known valid isomorphism** between the twist and E.
+        let point = Self::new([
             x.conjugate() * ENDO_U,
             y.conjugate() * ENDO_V,
             z.conjugate(),
-        ])
+        ]);
+        point.unwrap()
     }
 
     /// 𝜓(P) = 𝑢P, where 𝑢 = SEED of the curve
@@ -137,13 +160,23 @@ mod tests {
         FieldElement::from_hex_unchecked("0")
     ]);
 
-    // Cmoputes the psi^2() 'Untwist Frobenius Endomorphism'
+    /// Computes the psi^2() 'Untwist Frobenius Endomorphism'
+    ///  
+    /// # Safety
+    ///
+    /// - This function assumes `p` is a valid point on the BLS12-381 twist curve.
+    /// - The transformation involves multiplying the x and y coordinates by known constants.
+    /// - `unwrap()` is used because the resulting point remains valid under the curve equations.
     fn psi_square(
         p: &ShortWeierstrassProjectivePoint<BLS12381TwistCurve>,
     ) -> ShortWeierstrassProjectivePoint<BLS12381TwistCurve> {
         let [x, y, z] = p.coordinates();
         // Since power of frobenius map is 2 we apply once as applying twice is inverse
-        ShortWeierstrassProjectivePoint::new([x * ENDO_U_2, y * ENDO_V_2, z.clone()])
+
+        // SAFETY:
+        // - `ENDO_U_2` and `ENDO_V_2` are known valid constants.
+        // - `unwrap()` is safe because the transformation preserves curve validity.
+        ShortWeierstrassProjectivePoint::new([x * ENDO_U_2, y * ENDO_V_2, z.clone()]).unwrap()
     }
 
     #[allow(clippy::upper_case_acronyms)]
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/pairing.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/pairing.rs
index 0cf9ea067..2f6385d2d 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/pairing.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/pairing.rs
@@ -449,15 +449,19 @@ mod tests {
 
     #[test]
     fn ate_pairing_errors_when_one_element_is_not_in_subgroup() {
+        // p = (0, 2, 1) is in the curve but not in the subgroup.
+        // Recall that the BLS 12-381 curve equation is y^2 = x^3 + 4.
         let p = ShortWeierstrassProjectivePoint::new([
+            FieldElement::zero(),
+            FieldElement::from(2),
             FieldElement::one(),
-            FieldElement::one(),
-            FieldElement::one(),
-        ]);
+        ])
+        .unwrap();
         let q = ShortWeierstrassProjectivePoint::neutral_element();
         let result = BLS12381AtePairing::compute_batch(&[(&p.to_affine(), &q)]);
         assert!(result.is_err())
     }
+
     #[test]
     fn apply_12_times_frobenius_is_identity() {
         let f = Fp12E::from_coefficients(&[
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/twist.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/twist.rs
index 668e94359..ccce55866 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/twist.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/twist.rs
@@ -22,7 +22,10 @@ impl IsEllipticCurve for BLS12381TwistCurve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::new([
                 FieldElement::new(GENERATOR_X_0),
                 FieldElement::new(GENERATOR_X_1),
@@ -32,7 +35,8 @@ impl IsEllipticCurve for BLS12381TwistCurve {
                 FieldElement::new(GENERATOR_Y_1),
             ]),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bn_254/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/bn_254/curve.rs
index 02ede27c4..dd6f25652 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bn_254/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bn_254/curve.rs
@@ -20,12 +20,22 @@ impl IsEllipticCurve for BN254Curve {
     type BaseField = BN254PrimeField;
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
+    /// Returns the generator point of the BN254 curve.
+    ///
+    /// # Safety
+    ///
+    /// - The generator point is mathematically verified to be a valid point on the curve.
+    /// - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator coordinates `(1, 2, 1)` are **predefined** and belong to the BN254 curve.
+        // - `unwrap()` is safe because we **ensure** the input values satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::one(),
             FieldElement::<Self::BaseField>::from(2),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
@@ -50,13 +60,23 @@ impl ShortWeierstrassProjectivePoint<BN254TwistCurve> {
     /// We also use phi at the last lines of the Miller Loop of the pairing.
     /// phi(q) = (x^p, y^p, z^p), where (x, y, z) are the projective coordinates of q.
     /// See https://hackmd.io/@Wimet/ry7z1Xj-2#Subgroup-Checks.
+    ///
+    /// # Safety
+    ///
+    /// - The function assumes `self` is a valid point on the BN254 twist curve.
+    /// - The transformation follows a known isomorphism and preserves validity.
     pub fn phi(&self) -> Self {
         let [x, y, z] = self.coordinates();
-        Self::new([
+        // SAFETY:
+        // - `conjugate()` preserves the validity of the field element.
+        // - `unwrap()` is safe because the transformation follows
+        //   **a known valid isomorphism** between the twist and E.
+        let point = Self::new([
             x.conjugate() * GAMMA_12,
             y.conjugate() * GAMMA_13,
             z.conjugate(),
-        ])
+        ]);
+        point.unwrap()
     }
 
     // Checks if a G2 point is in the subgroup of the twisted curve.
@@ -279,7 +299,8 @@ mod tests {
                 )),
             ]),
             Fp2E::one(),
-        ]);
+        ])
+        .unwrap();
         assert!(!q.is_in_subgroup())
     }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bn_254/pairing.rs b/math/src/elliptic_curve/short_weierstrass/curves/bn_254/pairing.rs
index b550d5013..afdc37e7d 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bn_254/pairing.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bn_254/pairing.rs
@@ -304,13 +304,19 @@ fn line_optimized(p: &G1Point, t: &G2Point, q: &G2Point) -> (G2Point, Fp12E) {
         let y_r = g.square() - (e_square.double() + e_square);
         let z_r = b * &h;
 
+        debug_assert_eq!(
+            BN254TwistCurve::defining_equation_projective(&x_r, &y_r, &z_r),
+            Fp2E::zero(),
+        );
+        // SAFETY: `unwrap()` is used here because we ensure that `x_r, y_r, z_r`
+        // satisfy the curve equation. The previous assertion checks that this is indeed the case.
         let r = G2Point::new([x_r, y_r, z_r]);
 
         let l = Fp12E::new([
             Fp6E::new([y_p * (-h), Fp2E::zero(), Fp2E::zero()]),
             Fp6E::new([x_p * (j.double() + &j), i, Fp2E::zero()]),
         ]);
-        (r, l)
+        (r.unwrap(), l)
     } else {
         let [x_q, y_q, _] = q.coordinates();
         let [x_t, y_t, z_t] = t.coordinates();
@@ -332,13 +338,19 @@ fn line_optimized(p: &G1Point, t: &G2Point, q: &G2Point) -> (G2Point, Fp12E) {
         let y_r = &theta * (g - h) - i;
         let z_r = z_t * e;
 
+        debug_assert_eq!(
+            BN254TwistCurve::defining_equation_projective(&x_r, &y_r, &z_r),
+            Fp2E::zero()
+        );
+        // SAFETY: `unwrap()` is used here because we ensure that `x_r, y_r, z_r`
+        // satisfy the curve equation. The previous assertion checks that this is indeed the case.
         let r = G2Point::new([x_r, y_r, z_r]);
 
         let l = Fp12E::new([
             Fp6E::new([y_p * lambda, Fp2E::zero(), Fp2E::zero()]),
             Fp6E::new([x_p * (-theta), j, Fp2E::zero()]),
         ]);
-        (r, l)
+        (r.unwrap(), l)
     }
 }
 
@@ -658,7 +670,8 @@ mod tests {
                 )),
             ]),
             Fp2E::one(),
-        ]);
+        ])
+        .unwrap();
         let result = BN254AtePairing::compute_batch(&[(&p, &q)]);
         assert!(result.is_err())
     }
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bn_254/twist.rs b/math/src/elliptic_curve/short_weierstrass/curves/bn_254/twist.rs
index 3f5068908..d9df4dcab 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bn_254/twist.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bn_254/twist.rs
@@ -29,7 +29,10 @@ impl IsEllipticCurve for BN254TwistCurve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::new([
                 FieldElement::new(GENERATOR_X_0),
                 FieldElement::new(GENERATOR_X_1),
@@ -39,7 +42,8 @@ impl IsEllipticCurve for BN254TwistCurve {
                 FieldElement::new(GENERATOR_Y_1),
             ]),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/grumpkin/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/grumpkin/curve.rs
index b8d9fff5f..e0d1e5e3d 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/grumpkin/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/grumpkin/curve.rs
@@ -23,13 +23,17 @@ impl IsEllipticCurve for GrumpkinCurve {
 
     // G = (1, sprt(-16)) = (1, 17631683881184975370165255887551781615748388533673675138860) = (0x1, 0x2cf135e7506a45d632d270d45f1181294833fc48d823f272c)
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::one(),
             FieldElement::<Self::BaseField>::from_hex_unchecked(
                 "0x2cf135e7506a45d632d270d45f1181294833fc48d823f272c",
             ),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/pallas/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/pallas/curve.rs
index d932d7522..62d5a5604 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/pallas/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/pallas/curve.rs
@@ -13,11 +13,15 @@ impl IsEllipticCurve for PallasCurve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             -FieldElement::<Self::BaseField>::one(),
             FieldElement::<Self::BaseField>::from(2),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/secp256k1/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/secp256k1/curve.rs
index 99ad0f610..b885ecb9f 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/secp256k1/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/secp256k1/curve.rs
@@ -14,7 +14,10 @@ impl IsEllipticCurve for Secp256k1Curve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::from_hex_unchecked(
                 "79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798",
             ),
@@ -22,7 +25,8 @@ impl IsEllipticCurve for Secp256k1Curve {
                 "483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8",
             ),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/secp256r1/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/secp256r1/curve.rs
index 1fcc358f5..0d643f706 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/secp256r1/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/secp256r1/curve.rs
@@ -14,7 +14,10 @@ impl IsEllipticCurve for Secp256r1Curve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::from_hex_unchecked(
                 "6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296",
             ),
@@ -22,7 +25,8 @@ impl IsEllipticCurve for Secp256r1Curve {
                 "4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5",
             ),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/secq256k1/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/secq256k1/curve.rs
index e52997294..8101aadc7 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/secq256k1/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/secq256k1/curve.rs
@@ -13,7 +13,10 @@ impl IsEllipticCurve for Secq256k1Curve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::from_hex_unchecked(
                 "76C39F5585CB160EB6B06C87A2CE32E23134E45A097781A6A24288E37702EDA6",
             ),
@@ -21,7 +24,8 @@ impl IsEllipticCurve for Secq256k1Curve {
                 "3FFC646C7B2918B5DC2D265A8E82A7F7D18983D26E8DC055A4120DDAD952677F",
             ),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
@@ -68,7 +72,7 @@ mod tests {
         let z = FE::from_hex_unchecked(
             "bb26eae3d2b9603d98dff86d87175f442e539c07bbe4ef5712e47c4d72c89734",
         );
-        ShortWeierstrassProjectivePoint::<Secq256k1Curve>::new([x, y, z])
+        ShortWeierstrassProjectivePoint::<Secq256k1Curve>::new([x, y, z]).unwrap()
     }
 
     #[test]
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/stark_curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/stark_curve.rs
index c006e3682..a6adea101 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/stark_curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/stark_curve.rs
@@ -1,5 +1,6 @@
 use crate::{
     elliptic_curve::{
+        point::ProjectivePoint,
         short_weierstrass::{point::ShortWeierstrassProjectivePoint, traits::IsShortWeierstrass},
         traits::IsEllipticCurve,
     },
@@ -16,11 +17,11 @@ impl StarkCurve {
         x_hex: &str,
         y_hex: &str,
     ) -> ShortWeierstrassProjectivePoint<Self> {
-        ShortWeierstrassProjectivePoint::new([
+        ShortWeierstrassProjectivePoint(ProjectivePoint::new([
             FieldElement::<Stark252PrimeField>::from_hex_unchecked(x_hex),
             FieldElement::<Stark252PrimeField>::from_hex_unchecked(y_hex),
             FieldElement::<Stark252PrimeField>::from_hex_unchecked("1"),
-        ])
+        ]))
     }
 }
 
@@ -29,7 +30,10 @@ impl IsEllipticCurve for StarkCurve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::<Self::BaseField>::from_hex_unchecked(
                 "1EF15C18599971B7BECED415A40F0C7DEACFD9B0D1819E03D723D8BC943CFCA",
             ),
@@ -37,7 +41,8 @@ impl IsEllipticCurve for StarkCurve {
                 "5668060AA49730B7BE4801DF46EC62DE53ECD11ABE43A32873000C36E8DC1F",
             ),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/test_curve_1.rs b/math/src/elliptic_curve/short_weierstrass/curves/test_curve_1.rs
index c8601ce23..2612fea28 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/test_curve_1.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/test_curve_1.rs
@@ -41,11 +41,15 @@ impl IsEllipticCurve for TestCurve1 {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::from(35),
             FieldElement::from(31),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/test_curve_2.rs b/math/src/elliptic_curve/short_weierstrass/curves/test_curve_2.rs
index ca71c4a9d..0357ab41c 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/test_curve_2.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/test_curve_2.rs
@@ -49,7 +49,10 @@ impl IsEllipticCurve for TestCurve2 {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             FieldElement::new([
                 FieldElement::new(U384::from_hex_unchecked(
                     "21acedb641ca6d0f8b60148123a999801",
@@ -67,7 +70,8 @@ impl IsEllipticCurve for TestCurve2 {
                 )),
             ]),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/vesta/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/vesta/curve.rs
index ac871a3d6..c76c079e2 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/vesta/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/vesta/curve.rs
@@ -13,11 +13,15 @@ impl IsEllipticCurve for VestaCurve {
     type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
 
     fn generator() -> Self::PointRepresentation {
-        Self::PointRepresentation::new([
+        // SAFETY:
+        // - The generator point is mathematically verified to be a valid point on the curve.
+        // - `unwrap()` is safe because the provided coordinates satisfy the curve equation.
+        let point = Self::PointRepresentation::new([
             -FieldElement::<Self::BaseField>::one(),
             FieldElement::<Self::BaseField>::from(2),
             FieldElement::one(),
-        ])
+        ]);
+        point.unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/point.rs b/math/src/elliptic_curve/short_weierstrass/point.rs
index 679edeaea..2a82e4c1e 100644
--- a/math/src/elliptic_curve/short_weierstrass/point.rs
+++ b/math/src/elliptic_curve/short_weierstrass/point.rs
@@ -15,13 +15,32 @@ use super::traits::IsShortWeierstrass;
 use crate::traits::AsBytes;
 #[cfg(feature = "alloc")]
 use alloc::vec::Vec;
+
 #[derive(Clone, Debug)]
 pub struct ShortWeierstrassProjectivePoint<E: IsEllipticCurve>(pub ProjectivePoint<E>);
 
 impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
     /// Creates an elliptic curve point giving the projective [x: y: z] coordinates.
-    pub const fn new(value: [FieldElement<E::BaseField>; 3]) -> Self {
-        Self(ProjectivePoint::new(value))
+    pub fn new(value: [FieldElement<E::BaseField>; 3]) -> Result<Self, EllipticCurveError> {
+        let (x, y, z) = (&value[0], &value[1], &value[2]);
+
+        if z != &FieldElement::<E::BaseField>::zero()
+            && E::defining_equation_projective(x, y, z) == FieldElement::<E::BaseField>::zero()
+        {
+            Ok(Self(ProjectivePoint::new(value)))
+        // The point at infinity is (0, 1, 0)
+        // We convert every (0, _, 0) into the infinity.
+        } else if x == &FieldElement::<E::BaseField>::zero()
+            && z == &FieldElement::<E::BaseField>::zero()
+        {
+            Ok(Self(ProjectivePoint::new([
+                FieldElement::<E::BaseField>::zero(),
+                FieldElement::<E::BaseField>::one(),
+                FieldElement::<E::BaseField>::zero(),
+            ])))
+        } else {
+            Err(EllipticCurveError::InvalidPoint)
+        }
     }
 
     /// Returns the `x` coordinate of the point.
@@ -44,9 +63,6 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
         self.0.coordinates()
     }
 
-    /// Creates the same point in affine coordinates. That is,
-    /// returns [x / z: y / z: 1] where `self` is [x: y: z].
-    /// Panics if `self` is the point at infinity.
     pub fn to_affine(&self) -> Self {
         Self(self.0.to_affine())
     }
@@ -85,7 +101,15 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
         let xp = &hs + &hs;
         let yp = w * (four_b - &h) - eight_pys_square;
         let zp = eight_s_cube;
-        Self::new([xp, yp, zp])
+
+        debug_assert_eq!(
+            E::defining_equation_projective(&xp, &yp, &zp),
+            FieldElement::<E::BaseField>::zero()
+        );
+        // SAFETY: The values `x_p, y_p, z_p` are computed correctly to be on the curve.
+        // The assertion above verifies that the resulting point is valid.
+        let point = Self::new([xp, yp, zp]);
+        point.unwrap()
     }
     // https://hyperelliptic.org/EFD/g1p/data/shortw/projective/addition/madd-1998-cmo
     pub fn operate_with_affine(&self, other: &Self) -> Self {
@@ -103,11 +127,13 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
 
         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();
             }
@@ -125,7 +151,14 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
         let y = &u * (&r - &a) - &vvv * py;
         let z = &vvv * pz;
 
-        Self::new([x, y, z])
+        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()
     }
 }
 
@@ -141,24 +174,24 @@ impl<E: IsShortWeierstrass> FromAffine<E::BaseField> for ShortWeierstrassProject
     fn from_affine(
         x: FieldElement<E::BaseField>,
         y: FieldElement<E::BaseField>,
-    ) -> Result<Self, crate::elliptic_curve::traits::EllipticCurveError> {
-        if E::defining_equation(&x, &y) != FieldElement::zero() {
-            Err(EllipticCurveError::InvalidPoint)
-        } else {
-            let coordinates = [x, y, FieldElement::one()];
-            Ok(ShortWeierstrassProjectivePoint::new(coordinates))
-        }
+    ) -> Result<Self, EllipticCurveError> {
+        let coordinates = [x, y, FieldElement::one()];
+        ShortWeierstrassProjectivePoint::new(coordinates)
     }
 }
 
 impl<E: IsShortWeierstrass> IsGroup for ShortWeierstrassProjectivePoint<E> {
     /// The point at infinity.
     fn neutral_element() -> Self {
+        // SAFETY:
+        // - `(0, 1, 0)` is **mathematically valid** as the neutral element.
+        // - `unwrap()` is safe because this is **a known valid point**.
         Self::new([
             FieldElement::zero(),
             FieldElement::one(),
             FieldElement::zero(),
         ])
+        .unwrap()
     }
 
     fn is_neutral_element(&self) -> bool {
@@ -201,7 +234,14 @@ impl<E: IsShortWeierstrass> IsGroup for ShortWeierstrassProjectivePoint<E> {
                 let xp = &v * &a;
                 let yp = u * (&v_square_v2 - a) - &v_cube * u2;
                 let zp = &v_cube * w;
-                Self::new([xp, yp, zp])
+
+                debug_assert_eq!(
+                    E::defining_equation_projective(&xp, &yp, &zp),
+                    FieldElement::<E::BaseField>::zero()
+                );
+                // SAFETY: The values `x_p, y_p, z_p` are computed correctly to be on the curve.
+                // The assertion above verifies that the resulting point is valid.
+                Self::new([xp, yp, zp]).unwrap()
             }
         }
     }
@@ -209,7 +249,10 @@ impl<E: IsShortWeierstrass> IsGroup for ShortWeierstrassProjectivePoint<E> {
     /// Returns the additive inverse of the projective point `p`
     fn neg(&self) -> Self {
         let [px, py, pz] = self.coordinates();
-        Self::new([px.clone(), -py, pz.clone()])
+        // SAFETY:
+        // - Negating `y` maintains the curve structure.
+        // - `unwraps()` is safe because negation **is always valid**.
+        Self::new([px.clone(), -py, pz.clone()]).unwrap()
     }
 }
 
@@ -304,18 +347,9 @@ where
                     z = ByteConversion::from_bytes_le(&bytes[len * 2..])?;
                 }
 
-                if z == FieldElement::zero() {
-                    let point = Self::new([x, y, z]);
-                    if point.is_neutral_element() {
-                        Ok(point)
-                    } else {
-                        Err(DeserializationError::FieldFromBytesError)
-                    }
-                } else if E::defining_equation(&(&x / &z), &(&y / &z)) == FieldElement::zero() {
-                    Ok(Self::new([x, y, z]))
-                } else {
-                    Err(DeserializationError::FieldFromBytesError)
-                }
+                let point =
+                    Self::new([x, y, z]).map_err(|_| DeserializationError::FieldFromBytesError)?;
+                Ok(point)
             }
             PointFormat::Uncompressed => {
                 if bytes.len() % 2 != 0 {
@@ -334,11 +368,10 @@ where
                     y = ByteConversion::from_bytes_le(&bytes[len..])?;
                 }
 
-                if E::defining_equation(&x, &y) == FieldElement::zero() {
-                    Ok(Self::new([x, y, FieldElement::one()]))
-                } else {
-                    Err(DeserializationError::FieldFromBytesError)
-                }
+                let z = FieldElement::<E::BaseField>::one();
+                let point =
+                    Self::new([x, y, z]).map_err(|_| DeserializationError::FieldFromBytesError)?;
+                Ok(point)
             }
         }
     }
@@ -383,9 +416,25 @@ where
 pub struct ShortWeierstrassJacobianPoint<E: IsEllipticCurve>(pub JacobianPoint<E>);
 
 impl<E: IsShortWeierstrass> ShortWeierstrassJacobianPoint<E> {
-    /// Creates an elliptic curve point giving the projective [x: y: z] coordinates.
-    pub const fn new(value: [FieldElement<E::BaseField>; 3]) -> Self {
-        Self(JacobianPoint::new(value))
+    /// Creates an elliptic curve point giving the jacobian [x: y: z] coordinates.
+    pub fn new(value: [FieldElement<E::BaseField>; 3]) -> Result<Self, EllipticCurveError> {
+        let (x, y, z) = (&value[0], &value[1], &value[2]);
+
+        if z != &FieldElement::<E::BaseField>::zero()
+            && E::defining_equation_jacobian(x, y, z) == FieldElement::<E::BaseField>::zero()
+        {
+            Ok(Self(JacobianPoint::new(value)))
+        // The point at infinity is (1, 1, 0)
+        // We convert every (x, x, 0) into the infinity.
+        } else if z == &FieldElement::<E::BaseField>::zero() && x == y {
+            Ok(Self(JacobianPoint::new([
+                FieldElement::<E::BaseField>::one(),
+                FieldElement::<E::BaseField>::one(),
+                FieldElement::<E::BaseField>::zero(),
+            ])))
+        } else {
+            Err(EllipticCurveError::InvalidPoint)
+        }
     }
 
     /// Returns the `x` coordinate of the point.
@@ -435,7 +484,13 @@ impl<E: IsShortWeierstrass> ShortWeierstrassJacobianPoint<E> {
             let y3 = &e * (&d - &x3) - &c.double().double().double(); // Y3 = E * (D - X3) - 8 * C
             let z3 = (y1 * z1).double(); // Z3 = 2 * Y1 * Z1
 
-            Self::new([x3, y3, z3])
+            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()
         } else {
             // http://www.hyperelliptic.org/EFD/g1p/data/shortw/jacobian-0/doubling/dbl-2009-alnr
             // http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
@@ -449,7 +504,13 @@ impl<E: IsShortWeierstrass> ShortWeierstrassJacobianPoint<E> {
             let y3 = m * (&s - &x3) - &yyyy.double().double().double(); // Y3 = M * (S - X3) - 8 * YYYY
             let z3 = (y1 + z1).square() - &yy - &zz; // Z3 = (Y1 + Z1)^2 - YY - ZZ
 
-            Self::new([x3, y3, z3])
+            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()
         }
     }
 
@@ -486,7 +547,13 @@ impl<E: IsShortWeierstrass> ShortWeierstrassJacobianPoint<E> {
             let y3 = r * (&v - &x3) - s1 * &hhh;
             let z3 = z1 * &h;
 
-            Self::new([x3, y3, z3])
+            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()
         }
     }
 }
@@ -503,24 +570,25 @@ impl<E: IsShortWeierstrass> FromAffine<E::BaseField> for ShortWeierstrassJacobia
     fn from_affine(
         x: FieldElement<E::BaseField>,
         y: FieldElement<E::BaseField>,
-    ) -> Result<Self, crate::elliptic_curve::traits::EllipticCurveError> {
-        if E::defining_equation(&x, &y) != FieldElement::zero() {
-            Err(EllipticCurveError::InvalidPoint)
-        } else {
-            let coordinates = [x, y, FieldElement::one()];
-            Ok(ShortWeierstrassJacobianPoint::new(coordinates))
-        }
+    ) -> Result<Self, EllipticCurveError> {
+        let coordinates = [x, y, FieldElement::one()];
+        ShortWeierstrassJacobianPoint::new(coordinates)
     }
 }
 
 impl<E: IsShortWeierstrass> IsGroup for ShortWeierstrassJacobianPoint<E> {
     /// The point at infinity.
     fn neutral_element() -> Self {
+        // SAFETY:
+        // - `(1, 1, 0)` is **mathematically valid** as the neutral element.
+        // - `unwrap()` is safe because this is **a known valid point**.
+
         Self::new([
             FieldElement::one(),
             FieldElement::one(),
             FieldElement::zero(),
         ])
+        .unwrap()
     }
 
     fn is_neutral_element(&self) -> bool {
@@ -584,13 +652,22 @@ impl<E: IsShortWeierstrass> IsGroup for ShortWeierstrassJacobianPoint<E> {
         let z3 = z1 * z2;
         let z3 = z3.double() * h;
 
-        Self::new([x3, y3, z3])
+        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()
     }
 
     /// Returns the additive inverse of the jacobian point `p`
     fn neg(&self) -> Self {
         let [x, y, z] = self.coordinates();
-        Self::new([x.clone(), -y, z.clone()])
+        // SAFETY:
+        // - The negation formula for Short Weierstrass curves is well-defined.
+        // - The result remains a valid curve point.
+        Self::new([x.clone(), -y, z.clone()]).unwrap()
     }
 }
 
diff --git a/math/src/elliptic_curve/short_weierstrass/traits.rs b/math/src/elliptic_curve/short_weierstrass/traits.rs
index be456a3bd..a21f2edb3 100644
--- a/math/src/elliptic_curve/short_weierstrass/traits.rs
+++ b/math/src/elliptic_curve/short_weierstrass/traits.rs
@@ -18,6 +18,28 @@ pub trait IsShortWeierstrass: IsEllipticCurve + Clone + Debug {
     ) -> FieldElement<Self::BaseField> {
         y.square() - ((x.square() + Self::a()) * x + Self::b())
     }
+
+    // Evaluates the projective equation:
+    // y^2 * z = x^3 + a * x * z^2 + b * z^3
+    fn defining_equation_projective(
+        x: &FieldElement<Self::BaseField>,
+        y: &FieldElement<Self::BaseField>,
+        z: &FieldElement<Self::BaseField>,
+    ) -> FieldElement<Self::BaseField> {
+        y.square() * z - ((x.square() + Self::a() * z.square()) * x + Self::b() * z.square() * z)
+    }
+
+    // Evaluates the jacobian equation:
+    // y^2  = x^3 + a * x * z^4 + b * z^6
+    fn defining_equation_jacobian(
+        x: &FieldElement<Self::BaseField>,
+        y: &FieldElement<Self::BaseField>,
+        z: &FieldElement<Self::BaseField>,
+    ) -> FieldElement<Self::BaseField> {
+        y.square()
+            - ((x.square() + Self::a() * z.square().square()) * x
+                + Self::b() * z.square().square() * z.square())
+    }
 }
 
 pub trait Compress {