diff --git a/README.md b/README.md
index befe19dba..1edb80a8f 100644
--- a/README.md
+++ b/README.md
@@ -15,6 +15,16 @@ This library provides efficient implementation of cryptographic primitives used
 
 </div>
 
+## Examples - mini apps
+
+Below is a list of examples to understand lambdaworks and learn what you can build with the tools provided. 
+
+- [Merkle Tree CLI](./examples/merkle-tree-cli/)
+- [Proving Miden](./examples/prove-miden/)
+- [Shamir's secret sharing](./examples/shamir_secret_sharing/)
+- [BabySNARK](./examples/baby-snark/)
+- [Pinocchio](./examples/pinocchio/)
+
 ## Why we built lambdaworks
 
 Zero-Knowledge and Validity Proofs have gained a lot of attention over the last few years. We strongly believe in this potential and that is why we decided to start working in this challenging ecosystem, where math, cryptography and distributed systems meet. The main barrier in the beginning was not the cryptography or math but the lack of good libraries which are performant and developer friendly. There are some exceptions, though, like gnark or halo2. Some have nice APIs and are easy to work with, but they are not written in Rust, and some are written in Rust but have poor programming and engineering practices. Most of them don't have support for CUDA, Metal and WebGPU or distributed FFT calculation using schedulers like Dask.
@@ -41,14 +51,6 @@ Most of math and crypto crates supports no-std without allocation with `no-defau
 
 Both Math and Crypto support wasm with target `wasm32-unknown-unknown`. To see an example of how to use this to deploy a verifier in a browser, check the Cairo Prover wasm-pack verifier.
 
-## Examples - mini apps
-
-- [Merkle Tree CLI](https://github.com/lambdaclass/lambdaworks/tree/main/examples/merkle-tree-cli)
-- [Proving Miden](https://github.com/lambdaclass/lambdaworks/tree/main/examples/prove-miden)
-- [Shamir's secret sharing](https://github.com/lambdaclass/lambdaworks/tree/main/examples/shamir_secret_sharing)
-- [BabySNARK](https://github.com/lambdaclass/lambdaworks/tree/main/examples/baby-snark)
-- [Pinocchio](https://github.com/lambdaclass/lambdaworks/tree/main/examples/pinocchio)
-
 ## Exercises and Challenges
 
 - [lambdaworks exercises and challenges](https://github.com/lambdaclass/lambdaworks_exercises/tree/main)
diff --git a/math/benches/criterion_elliptic_curve.rs b/math/benches/criterion_elliptic_curve.rs
index 3812d74c3..f9bffc781 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_377_elliptic_curve_benchmarks,bls12_381_elliptic_curve_benchmarks
+    targets = bn_254_elliptic_curve_benchmarks,bls12_381_elliptic_curve_benchmarks,bls12_377_elliptic_curve_benchmarks
 );
 criterion_main!(elliptic_curve_benches);
diff --git a/math/benches/criterion_field.rs b/math/benches/criterion_field.rs
index 1c21822de..6e41cbb4b 100644
--- a/math/benches/criterion_field.rs
+++ b/math/benches/criterion_field.rs
@@ -2,7 +2,7 @@ use criterion::{criterion_group, criterion_main, Criterion};
 use pprof::criterion::{Output, PProfProfiler};
 
 mod fields;
-use fields::mersenne31::mersenne31_ops_benchmarks;
+use fields::mersenne31::{mersenne31_extension_ops_benchmarks, mersenne31_ops_benchmarks};
 use fields::mersenne31_montgomery::mersenne31_mont_ops_benchmarks;
 use fields::{
     stark252::starkfield_ops_benchmarks, u64_goldilocks::u64_goldilocks_ops_benchmarks,
@@ -12,6 +12,6 @@ use fields::{
 criterion_group!(
     name = field_benches;
     config = Criterion::default().with_profiler(PProfProfiler::new(100, Output::Flamegraph(None)));
-    targets = starkfield_ops_benchmarks, mersenne31_ops_benchmarks, mersenne31_mont_ops_benchmarks, u64_goldilocks_ops_benchmarks, u64_goldilocks_montgomery_ops_benchmarks
+    targets = mersenne31_ops_benchmarks, mersenne31_extension_ops_benchmarks, mersenne31_mont_ops_benchmarks, starkfield_ops_benchmarks, u64_goldilocks_ops_benchmarks, u64_goldilocks_montgomery_ops_benchmarks
 );
 criterion_main!(field_benches);
diff --git a/math/benches/elliptic_curves/bls12_381.rs b/math/benches/elliptic_curves/bls12_381.rs
index d33e8d2ca..452480697 100644
--- a/math/benches/elliptic_curves/bls12_381.rs
+++ b/math/benches/elliptic_curves/bls12_381.rs
@@ -4,7 +4,9 @@ use lambdaworks_math::{
     elliptic_curve::{
         short_weierstrass::{
             curves::bls12_381::{
-                curve::BLS12381Curve, pairing::BLS12381AtePairing, twist::BLS12381TwistCurve,
+                curve::BLS12381Curve,
+                pairing::{final_exponentiation, miller, BLS12381AtePairing},
+                twist::BLS12381TwistCurve,
             },
             traits::Compress,
         },
@@ -24,6 +26,8 @@ pub fn bls12_381_elliptic_curve_benchmarks(c: &mut Criterion) {
     let a_g2 = BLS12381TwistCurve::generator();
     let b_g2 = BLS12381TwistCurve::generator();
 
+    let miller_loop_output = miller(&a_g2, &a_g1);
+
     let mut group = c.benchmark_group("BLS12-381 Ops");
     group.significance_level(0.1).sample_size(10000);
     group.throughput(criterion::Throughput::Elements(1));
@@ -93,4 +97,14 @@ pub fn bls12_381_elliptic_curve_benchmarks(c: &mut Criterion) {
             ))
         });
     });
+
+    // Miller
+    group.bench_function("Miller", |bencher| {
+        bencher.iter(|| black_box(miller(black_box(&a_g2), black_box(&a_g1))))
+    });
+
+    // Final Exponentiation Optimized
+    group.bench_function("Final Exponentiation", |bencher| {
+        bencher.iter(|| black_box(final_exponentiation(black_box(&miller_loop_output))))
+    });
 }
diff --git a/math/benches/fields/mersenne31.rs b/math/benches/fields/mersenne31.rs
index 99e3921a5..e8d99d1c2 100644
--- a/math/benches/fields/mersenne31.rs
+++ b/math/benches/fields/mersenne31.rs
@@ -1,10 +1,18 @@
 use std::hint::black_box;
 
 use criterion::Criterion;
-use lambdaworks_math::field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field};
+use lambdaworks_math::field::{
+    element::FieldElement,
+    fields::mersenne31::{
+        extensions::{Degree2ExtensionField, Degree4ExtensionField},
+        field::Mersenne31Field,
+    },
+};
 use rand::random;
 
 pub type F = FieldElement<Mersenne31Field>;
+pub type Fp2E = FieldElement<Degree2ExtensionField>;
+pub type Fp4E = FieldElement<Degree4ExtensionField>;
 
 #[inline(never)]
 #[no_mangle]
@@ -17,6 +25,60 @@ pub fn rand_field_elements(num: usize) -> Vec<(F, F)> {
     result
 }
 
+//TODO: Check if this is the correct way to bench.
+pub fn rand_fp4e(num: usize) -> Vec<(Fp4E, Fp4E)> {
+    let mut result = Vec::with_capacity(num);
+    for _ in 0..result.capacity() {
+        result.push((
+            Fp4E::new([
+                Fp2E::new([F::new(random()), F::new(random())]),
+                Fp2E::new([F::new(random()), F::new(random())]),
+            ]),
+            Fp4E::new([
+                Fp2E::new([F::new(random()), F::new(random())]),
+                Fp2E::new([F::new(random()), F::new(random())]),
+            ]),
+        ));
+    }
+    result
+}
+
+pub fn mersenne31_extension_ops_benchmarks(c: &mut Criterion) {
+    let input: Vec<Vec<(Fp4E, Fp4E)>> = [1000000].into_iter().map(rand_fp4e).collect::<Vec<_>>();
+
+    let mut group = c.benchmark_group("Mersenne31 Fp4 operations");
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("Mul of Fp4 {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, y) in i {
+                    black_box(black_box(x) * black_box(y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("Square of Fp4 {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).square());
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("Inv of Fp4 {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).inv().unwrap());
+                }
+            });
+        });
+    }
+}
+
 pub fn mersenne31_ops_benchmarks(c: &mut Criterion) {
     let input: Vec<Vec<(F, F)>> = [1, 10, 100, 1000, 10000, 100000, 1000000]
         .into_iter()
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 c324b35de..efa52bdbc 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
@@ -1,10 +1,25 @@
-use super::field_extension::BLS12377PrimeField;
+use super::{
+    field_extension::{BLS12377PrimeField, Degree2ExtensionField},
+    twist::BLS12377TwistCurve,
+};
+use crate::cyclic_group::IsGroup;
 use crate::elliptic_curve::short_weierstrass::point::ShortWeierstrassProjectivePoint;
 use crate::elliptic_curve::traits::IsEllipticCurve;
+use crate::unsigned_integer::element::U256;
+
 use crate::{
     elliptic_curve::short_weierstrass::traits::IsShortWeierstrass, field::element::FieldElement,
 };
 
+pub const SUBGROUP_ORDER: U256 =
+    U256::from_hex_unchecked("12ab655e9a2ca55660b44d1e5c37b00159aa76fed00000010a11800000000001");
+
+pub const CURVE_COFACTOR: U256 =
+    U256::from_hex_unchecked("0x30631250834960419227450344600217059328");
+
+pub type BLS12377FieldElement = FieldElement<BLS12377PrimeField>;
+pub type BLS12377TwistCurveFieldElement = FieldElement<Degree2ExtensionField>;
+
 /// The description of the curve.
 #[derive(Clone, Debug)]
 pub struct BLS12377Curve;
@@ -32,6 +47,71 @@ impl IsShortWeierstrass for BLS12377Curve {
     }
 }
 
+/// This is equal to the frobenius trace of the BLS12 377 curve minus one or seed value z.
+pub const MILLER_LOOP_CONSTANT: u64 = 0x8508c00000000001;
+
+/// 𝛽 : primitive cube root of unity of 𝐹ₚ that §satisfies the minimal equation
+/// 𝛽² + 𝛽 + 1 = 0 mod 𝑝
+pub const CUBE_ROOT_OF_UNITY_G1: BLS12377FieldElement = FieldElement::from_hex_unchecked(
+    "0x1ae3a4617c510eabc8756ba8f8c524eb8882a75cc9bc8e359064ee822fb5bffd1e945779fffffffffffffffffffffff",
+);
+
+/// x-coordinate of 𝜁 ∘ 𝜋_q ∘ 𝜁⁻¹, where 𝜁 is the isomorphism u:E'(𝔽ₚ₆) −> E(𝔽ₚ₁₂) from the twist to E
+pub const ENDO_U: BLS12377TwistCurveFieldElement =
+    BLS12377TwistCurveFieldElement::const_from_raw([
+        FieldElement::from_hex_unchecked(
+            "9B3AF05DD14F6EC619AAF7D34594AABC5ED1347970DEC00452217CC900000008508C00000000002",
+        ),
+        FieldElement::from_hex_unchecked("0"),
+    ]);
+
+/// y-coordinate of 𝜁 ∘ 𝜋_q ∘ 𝜁⁻¹, where 𝜁 is the isomorphism u:E'(𝔽ₚ₆) −> E(𝔽ₚ₁₂) from the twist to E
+pub const ENDO_V: BLS12377TwistCurveFieldElement =
+    BLS12377TwistCurveFieldElement::const_from_raw([
+        FieldElement::from_hex_unchecked("1680A40796537CAC0C534DB1A79BEB1400398F50AD1DEC1BCE649CF436B0F6299588459BFF27D8E6E76D5ECF1391C63"),
+        FieldElement::from_hex_unchecked("0"),
+    ]);
+
+impl ShortWeierstrassProjectivePoint<BLS12377Curve> {
+    /// Returns 𝜙(P) = (𝑥, 𝑦) ⇒ (𝛽𝑥, 𝑦), where 𝛽 is the Cube Root of Unity in the base prime field
+    /// https://eprint.iacr.org/2022/352.pdf 2 Preliminaries
+    fn phi(&self) -> Self {
+        let mut a = self.clone();
+        a.0.value[0] = a.x() * CUBE_ROOT_OF_UNITY_G1;
+        a
+    }
+
+    /// 𝜙(P) = −𝑢²P
+    /// https://eprint.iacr.org/2022/352.pdf 4.3 Prop. 4
+    pub fn is_in_subgroup(&self) -> bool {
+        self.operate_with_self(MILLER_LOOP_CONSTANT)
+            .operate_with_self(MILLER_LOOP_CONSTANT)
+            .neg()
+            == self.phi()
+    }
+}
+
+impl ShortWeierstrassProjectivePoint<BLS12377TwistCurve> {
+    /// 𝜓(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))
+    fn psi(&self) -> Self {
+        let [x, y, z] = self.coordinates();
+        Self::new([
+            x.conjugate() * ENDO_U,
+            y.conjugate() * ENDO_V,
+            z.conjugate(),
+        ])
+    }
+
+    /// 𝜓(P) = 𝑢P, where 𝑢 = SEED of the curve
+    /// https://eprint.iacr.org/2022/352.pdf 4.2
+    pub fn is_in_subgroup(&self) -> bool {
+        self.psi() == self.operate_with_self(MILLER_LOOP_CONSTANT)
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -43,17 +123,18 @@ mod tests {
     use super::BLS12377Curve;
 
     #[allow(clippy::upper_case_acronyms)]
-    type FEE = FieldElement<BLS12377PrimeField>;
+    type FpE = FieldElement<BLS12377PrimeField>;
+    type Fp2 = FieldElement<Degree2ExtensionField>;
 
     fn point_1() -> ShortWeierstrassProjectivePoint<BLS12377Curve> {
-        let x = FEE::new_base("134e4cc122cb62a06767fb98e86f2d5f77e2a12fefe23bb0c4c31d1bd5348b88d6f5e5dee2b54db4a2146cc9f249eea");
-        let y = FEE::new_base("17949c29effee7a9f13f69b1c28eccd78c1ed12b47068836473481ff818856594fd9c1935e3d9e621901a2d500257a2");
+        let x = FpE::new_base("134e4cc122cb62a06767fb98e86f2d5f77e2a12fefe23bb0c4c31d1bd5348b88d6f5e5dee2b54db4a2146cc9f249eea");
+        let y = FpE::new_base("17949c29effee7a9f13f69b1c28eccd78c1ed12b47068836473481ff818856594fd9c1935e3d9e621901a2d500257a2");
         BLS12377Curve::create_point_from_affine(x, y).unwrap()
     }
 
     fn point_1_times_5() -> ShortWeierstrassProjectivePoint<BLS12377Curve> {
-        let x = FEE::new_base("3c852d5aab73fbb51e57fbf5a0a8b5d6513ec922b2611b7547bfed74cba0dcdfc3ad2eac2733a4f55d198ec82b9964");
-        let y = FEE::new_base("a71425e68e55299c64d7eada9ae9c3fb87a9626b941d17128b64685fc07d0e635f3c3a512903b4e0a43e464045967b");
+        let x = FpE::new_base("3c852d5aab73fbb51e57fbf5a0a8b5d6513ec922b2611b7547bfed74cba0dcdfc3ad2eac2733a4f55d198ec82b9964");
+        let y = FpE::new_base("a71425e68e55299c64d7eada9ae9c3fb87a9626b941d17128b64685fc07d0e635f3c3a512903b4e0a43e464045967b");
         BLS12377Curve::create_point_from_affine(x, y).unwrap()
     }
 
@@ -101,9 +182,9 @@ mod tests {
         let point_1 = point_1().to_affine();
 
         // Create point 2
-        let x = FEE::new_base("134e4cc122cb62a06767fb98e86f2d5f77e2a12fefe23bb0c4c31d1bd5348b88d6f5e5dee2b54db4a2146cc9f249eea") * FEE::from(2);
-        let y = FEE::new_base("17949c29effee7a9f13f69b1c28eccd78c1ed12b47068836473481ff818856594fd9c1935e3d9e621901a2d500257a2") * FEE::from(2);
-        let z = FEE::from(2);
+        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 first_algorithm_result = point_2.operate_with(&point_1).to_affine();
@@ -115,15 +196,15 @@ mod tests {
     #[test]
     fn create_valid_point_works() {
         let p = point_1();
-        assert_eq!(*p.x(), FEE::new_base("134e4cc122cb62a06767fb98e86f2d5f77e2a12fefe23bb0c4c31d1bd5348b88d6f5e5dee2b54db4a2146cc9f249eea"));
-        assert_eq!(*p.y(), FEE::new_base("17949c29effee7a9f13f69b1c28eccd78c1ed12b47068836473481ff818856594fd9c1935e3d9e621901a2d500257a2"));
-        assert_eq!(*p.z(), FEE::new_base("1"));
+        assert_eq!(*p.x(), FpE::new_base("134e4cc122cb62a06767fb98e86f2d5f77e2a12fefe23bb0c4c31d1bd5348b88d6f5e5dee2b54db4a2146cc9f249eea"));
+        assert_eq!(*p.y(), FpE::new_base("17949c29effee7a9f13f69b1c28eccd78c1ed12b47068836473481ff818856594fd9c1935e3d9e621901a2d500257a2"));
+        assert_eq!(*p.z(), FpE::new_base("1"));
     }
 
     #[test]
     fn create_invalid_points_panics() {
         assert_eq!(
-            BLS12377Curve::create_point_from_affine(FEE::from(1), FEE::from(1)).unwrap_err(),
+            BLS12377Curve::create_point_from_affine(FpE::from(1), FpE::from(1)).unwrap_err(),
             EllipticCurveError::InvalidPoint
         )
     }
@@ -144,4 +225,36 @@ mod tests {
             g.operate_with_self(3_u16)
         );
     }
+
+    #[test]
+    fn generator_g1_is_in_subgroup() {
+        let g = BLS12377Curve::generator();
+        assert!(g.is_in_subgroup())
+    }
+
+    #[test]
+    fn point1_is_in_subgroup() {
+        let p = point_1();
+        assert!(p.is_in_subgroup())
+    }
+
+    #[test]
+    fn arbitrary_g1_point_is_in_subgroup() {
+        let g = BLS12377Curve::generator().operate_with_self(32u64);
+        assert!(g.is_in_subgroup())
+    }
+    #[test]
+    fn generator_g2_is_in_subgroup() {
+        let g = BLS12377TwistCurve::generator();
+        assert!(g.is_in_subgroup())
+    }
+
+    #[test]
+    fn g2_conjugate_works() {
+        let a = Fp2::zero();
+        let mut expected = a.conjugate();
+        expected = expected.conjugate();
+
+        assert_eq!(a, expected);
+    }
 }
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/field_extension.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/field_extension.rs
index 32f9d465d..96c227ddb 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/field_extension.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/field_extension.rs
@@ -1,10 +1,18 @@
 use crate::field::{
     element::FieldElement,
+    errors::FieldError,
+    extensions::{
+        cubic::{CubicExtensionField, HasCubicNonResidue},
+        quadratic::{HasQuadraticNonResidue, QuadraticExtensionField},
+    },
     fields::montgomery_backed_prime_fields::{IsModulus, MontgomeryBackendPrimeField},
+    traits::{IsField, IsSubFieldOf},
 };
+use crate::traits::ByteConversion;
 use crate::unsigned_integer::element::U384;
 
 pub const BLS12377_PRIME_FIELD_ORDER: U384 = U384::from_hex_unchecked("1ae3a4617c510eac63b05c06ca1493b1a22d9f300f5138f1ef3622fba094800170b5d44300000008508c00000000001");
+pub const FP2_RESIDUE: FieldElement<BLS12377PrimeField> =FieldElement::from_hex_unchecked("0x1ae3a4617c510eac63b05c06ca1493b1a22d9f300f5138f1ef3622fba094800170b5d44300000008508bffffffffffc");
 
 // FPBLS12377
 #[derive(Clone, Debug)]
@@ -15,8 +23,340 @@ impl IsModulus<U384> for BLS12377FieldModulus {
 
 pub type BLS12377PrimeField = MontgomeryBackendPrimeField<BLS12377FieldModulus, 6>;
 
+#[derive(Clone, Debug)]
+pub struct Degree2ExtensionField;
+
+impl IsField for Degree2ExtensionField {
+    type BaseType = [FieldElement<BLS12377PrimeField>; 2];
+    /// Returns the component wise addition of `a` and `b`
+    /// Returns the component wise addition of `a` and `b`
+    fn add(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        [&a[0] + &b[0], &a[1] + &b[1]]
+    }
+
+    /// Returns the multiplication of `a` and `b` using the following
+    /// equation:
+    /// (a0 + a1 * t) * (b0 + b1 * t) = a0 * b0 + a1 * b1 * Self::residue() + (a0 * b1 + a1 * b0) * t
+    /// where `t.pow(2)` equals `Q::residue()`.
+    fn mul(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        let a0b0 = &a[0] * &b[0];
+        let a1b1 = &a[1] * &b[1];
+        let z = (&a[0] + &a[1]) * (&b[0] + &b[1]);
+        [&a0b0 + &a1b1 * FP2_RESIDUE, z - a0b0 - a1b1]
+    }
+
+    fn square(a: &Self::BaseType) -> Self::BaseType {
+        let [a0, a1] = a;
+        let v0 = a0 * a1;
+        let c0 = (a0 + a1) * (a0 + &FP2_RESIDUE * a1) - &v0 - FP2_RESIDUE * &v0;
+        let c1 = &v0 + &v0;
+        [c0, c1]
+    }
+    /// Returns the component wise subtraction of `a` and `b`
+    fn sub(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        [&a[0] - &b[0], &a[1] - &b[1]]
+    }
+
+    /// Returns the component wise negation of `a`
+    fn neg(a: &Self::BaseType) -> Self::BaseType {
+        [-&a[0], -&a[1]]
+    }
+
+    /// Returns the multiplicative inverse of `a`
+    /// This uses the equality `(a0 + a1 * t) * (a0 - a1 * t) = a0.pow(2) - a1.pow(2) * Q::residue()`
+    fn inv(a: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
+        let q: FieldElement<BLS12377PrimeField> = -FieldElement::from(5);
+        let inv_norm = (a[0].pow(2_u64) - q * a[1].pow(2_u64)).inv()?;
+        Ok([&a[0] * &inv_norm, -&a[1] * inv_norm])
+    }
+
+    /// Returns the division of `a` and `b`
+    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    }
+
+    /// Returns a boolean indicating whether `a` and `b` are equal component wise.
+    fn eq(a: &Self::BaseType, b: &Self::BaseType) -> bool {
+        a[0] == b[0] && a[1] == b[1]
+    }
+
+    /// Returns the additive neutral element of the field extension.
+    fn zero() -> Self::BaseType {
+        [FieldElement::zero(), FieldElement::zero()]
+    }
+
+    /// Returns the multiplicative neutral element of the field extension.
+    fn one() -> Self::BaseType {
+        [FieldElement::one(), FieldElement::zero()]
+    }
+
+    /// Returns the element `x * 1` where 1 is the multiplicative neutral element.
+    fn from_u64(x: u64) -> Self::BaseType {
+        [FieldElement::from(x), FieldElement::zero()]
+    }
+
+    /// Takes as input an element of BaseType and returns the internal representation
+    /// of that element in the field.
+    /// Note: for this case this is simply the identity, because the components
+    /// already have correct representations.
+    fn from_base_type(x: Self::BaseType) -> Self::BaseType {
+        x
+    }
+}
+
+impl IsSubFieldOf<Degree2ExtensionField> for BLS12377PrimeField {
+    fn mul(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        let c0 = FieldElement::from_raw(<Self as IsField>::mul(a, b[0].value()));
+        let c1 = FieldElement::from_raw(<Self as IsField>::mul(a, b[1].value()));
+        [c0, c1]
+    }
+
+    fn add(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        let c0 = FieldElement::from_raw(<Self as IsField>::add(a, b[0].value()));
+        let c1 = FieldElement::from_raw(*b[1].value());
+        [c0, c1]
+    }
+
+    fn div(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        let b_inv = Degree2ExtensionField::inv(b).unwrap();
+        <Self as IsSubFieldOf<Degree2ExtensionField>>::mul(a, &b_inv)
+    }
+
+    fn sub(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        let c0 = FieldElement::from_raw(<Self as IsField>::sub(a, b[0].value()));
+        let c1 = FieldElement::from_raw(<Self as IsField>::neg(b[1].value()));
+        [c0, c1]
+    }
+
+    fn embed(a: Self::BaseType) -> <Degree2ExtensionField as IsField>::BaseType {
+        [FieldElement::from_raw(a), FieldElement::zero()]
+    }
+
+    #[cfg(feature = "alloc")]
+    fn to_subfield_vec(
+        b: <Degree2ExtensionField as IsField>::BaseType,
+    ) -> alloc::vec::Vec<Self::BaseType> {
+        b.into_iter().map(|x| x.to_raw()).collect()
+    }
+}
+
+impl ByteConversion for FieldElement<Degree2ExtensionField> {
+    #[cfg(feature = "alloc")]
+    fn to_bytes_be(&self) -> alloc::vec::Vec<u8> {
+        let mut byte_slice = ByteConversion::to_bytes_be(&self.value()[0]);
+        byte_slice.extend(ByteConversion::to_bytes_be(&self.value()[1]));
+        byte_slice
+    }
+
+    #[cfg(feature = "alloc")]
+    fn to_bytes_le(&self) -> alloc::vec::Vec<u8> {
+        let mut byte_slice = ByteConversion::to_bytes_le(&self.value()[0]);
+        byte_slice.extend(ByteConversion::to_bytes_le(&self.value()[1]));
+        byte_slice
+    }
+
+    fn from_bytes_be(bytes: &[u8]) -> Result<Self, crate::errors::ByteConversionError>
+    where
+        Self: core::marker::Sized,
+    {
+        const BYTES_PER_FIELD: usize = 48;
+        let x0 = FieldElement::from_bytes_be(&bytes[0..BYTES_PER_FIELD])?;
+        let x1 = FieldElement::from_bytes_be(&bytes[BYTES_PER_FIELD..BYTES_PER_FIELD * 2])?;
+        Ok(Self::new([x0, x1]))
+    }
+
+    fn from_bytes_le(bytes: &[u8]) -> Result<Self, crate::errors::ByteConversionError>
+    where
+        Self: core::marker::Sized,
+    {
+        const BYTES_PER_FIELD: usize = 48;
+        let x0 = FieldElement::from_bytes_le(&bytes[0..BYTES_PER_FIELD])?;
+        let x1 = FieldElement::from_bytes_le(&bytes[BYTES_PER_FIELD..BYTES_PER_FIELD * 2])?;
+        Ok(Self::new([x0, x1]))
+    }
+}
+
+impl FieldElement<Degree2ExtensionField> {
+    pub fn new_base(a_hex: &str) -> Self {
+        Self::new([FieldElement::new(U384::from(a_hex)), FieldElement::zero()])
+    }
+
+    pub fn conjugate(&self) -> Self {
+        let [a, b] = self.value();
+        Self::new([a.clone(), -b])
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct BLS12377Residue;
+impl HasQuadraticNonResidue<BLS12377PrimeField> for BLS12377Residue {
+    fn residue() -> FieldElement<BLS12377PrimeField> {
+        FP2_RESIDUE
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct LevelTwoResidue;
+
+impl HasCubicNonResidue<Degree2ExtensionField> for LevelTwoResidue {
+    fn residue() -> FieldElement<Degree2ExtensionField> {
+        FieldElement::new([
+            FieldElement::new(U384::from("0")),
+            -FieldElement::new(U384::from("1")),
+        ])
+    }
+}
+
+pub type Degree6ExtensionField = CubicExtensionField<Degree2ExtensionField, LevelTwoResidue>;
+
+#[derive(Debug, Clone)]
+pub struct LevelThreeResidue;
+impl HasQuadraticNonResidue<Degree6ExtensionField> for LevelThreeResidue {
+    fn residue() -> FieldElement<Degree6ExtensionField> {
+        FieldElement::new([
+            FieldElement::zero(),
+            FieldElement::one(),
+            FieldElement::zero(),
+        ])
+    }
+}
+
+/// We define Fp12 = Fp6 [w] / (w^2 - v)
+pub type Degree12ExtensionField = QuadraticExtensionField<Degree6ExtensionField, LevelThreeResidue>;
+
+impl FieldElement<Degree6ExtensionField> {
+    pub fn new_base(a_hex: &str) -> Self {
+        Self::new([
+            FieldElement::new([FieldElement::new(U384::from(a_hex)), FieldElement::zero()]),
+            FieldElement::zero(),
+            FieldElement::zero(),
+        ])
+    }
+}
+
 impl FieldElement<BLS12377PrimeField> {
     pub fn new_base(a_hex: &str) -> Self {
-        Self::new(U384::from_hex_unchecked(a_hex))
+        Self::new(U384::from(a_hex))
+    }
+}
+impl FieldElement<Degree12ExtensionField> {
+    pub fn new_base(a_hex: &str) -> Self {
+        Self::new([
+            FieldElement::<Degree6ExtensionField>::new_base(a_hex),
+            FieldElement::zero(),
+        ])
+    }
+
+    pub fn from_coefficients(coefficients: &[&str; 12]) -> Self {
+        FieldElement::<Degree12ExtensionField>::new([
+            FieldElement::new([
+                FieldElement::new([
+                    FieldElement::new(U384::from(coefficients[0])),
+                    FieldElement::new(U384::from(coefficients[1])),
+                ]),
+                FieldElement::new([
+                    FieldElement::new(U384::from(coefficients[2])),
+                    FieldElement::new(U384::from(coefficients[3])),
+                ]),
+                FieldElement::new([
+                    FieldElement::new(U384::from(coefficients[4])),
+                    FieldElement::new(U384::from(coefficients[5])),
+                ]),
+            ]),
+            FieldElement::new([
+                FieldElement::new([
+                    FieldElement::new(U384::from(coefficients[6])),
+                    FieldElement::new(U384::from(coefficients[7])),
+                ]),
+                FieldElement::new([
+                    FieldElement::new(U384::from(coefficients[8])),
+                    FieldElement::new(U384::from(coefficients[9])),
+                ]),
+                FieldElement::new([
+                    FieldElement::new(U384::from(coefficients[10])),
+                    FieldElement::new(U384::from(coefficients[11])),
+                ]),
+            ]),
+        ])
+    }
+}
+
+#[cfg(test)]
+mod tests {
+
+    use super::*;
+    type FpE = FieldElement<BLS12377PrimeField>;
+    type Fp2E = FieldElement<Degree2ExtensionField>;
+    type Fp6E = FieldElement<Degree6ExtensionField>;
+    type Fp12E = FieldElement<Degree12ExtensionField>;
+
+    #[test]
+    fn embed_base_field_with_degree_2_extension() {
+        let a = FpE::from(3);
+        let a_extension = Fp2E::from(3);
+        assert_eq!(a.to_extension::<Degree2ExtensionField>(), a_extension);
+    }
+    #[test]
+    fn add_base_field_with_degree_2_extension() {
+        let a = FpE::from(3);
+        let a_extension = Fp2E::from(3);
+        let b = Fp2E::from(2);
+        assert_eq!(a + &b, a_extension + b);
+    }
+    #[test]
+    fn mul_degree_2_with_degree_6_extension() {
+        let a = Fp2E::new([FpE::from(3), FpE::from(4)]);
+        let a_extension = a.clone().to_extension::<Degree2ExtensionField>();
+        let b = Fp6E::from(2);
+        assert_eq!(a * &b, a_extension * b);
+    }
+    #[test]
+    fn div_degree_6_degree_12_extension() {
+        let a = Fp6E::from(3);
+        let a_extension = Fp12E::from(3);
+        let b = Fp12E::from(2);
+        assert_eq!(a / &b, a_extension / b);
+    }
+
+    #[test]
+    fn double_equals_sum_two_times() {
+        let a = FpE::from(3);
+        assert_eq!(a.double(), a.clone() + a);
+    }
+    #[test]
+    fn base_field_sum_is_asociative() {
+        let a = FpE::from(3);
+        let b = FpE::from(2);
+        let c = &a + &b;
+        assert_eq!(a.double() + b, a + c);
+    }
+    #[test]
+    fn degree_2_extension_mul_is_conmutative() {
+        let a = Fp2E::from(3);
+        let b = Fp2E::new([FpE::from(2), FpE::from(4)]);
+        assert_eq!(&a * &b, b * a);
+    }
+
+    #[test]
+    fn base_field_pow_p_is_identity() {
+        let a = FpE::from(3);
+        assert_eq!(a.pow(BLS12377_PRIME_FIELD_ORDER), a);
+    }
+    #[test]
+    fn square_equals_mul_two_times() {
+        let a = FpE::from(3);
+        assert_eq!(a.square(), a.clone() * a);
     }
 }
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/mod.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/mod.rs
index 8fc713a88..5fe6ce8bd 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/mod.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/mod.rs
@@ -1,2 +1,3 @@
 pub mod curve;
 pub mod field_extension;
+pub mod twist;
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
new file mode 100644
index 000000000..e0aa87c59
--- /dev/null
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/twist.rs
@@ -0,0 +1,172 @@
+use super::field_extension::Degree2ExtensionField;
+use crate::elliptic_curve::short_weierstrass::point::ShortWeierstrassProjectivePoint;
+use crate::elliptic_curve::traits::IsEllipticCurve;
+use crate::unsigned_integer::element::U384;
+use crate::{
+    elliptic_curve::short_weierstrass::traits::IsShortWeierstrass, field::element::FieldElement,
+};
+
+const GENERATOR_X_0: U384 = U384::from_hex_unchecked("0x018480be71c785fec89630a2a3841d01c565f071203e50317ea501f557db6b9b71889f52bb53540274e3e48f7c005196");
+const GENERATOR_X_1: U384 = U384::from_hex_unchecked("0x00ea6040e700403170dc5a51b1b140d5532777ee6651cecbe7223ece0799c9de5cf89984bff76fe6b26bfefa6ea16afe");
+const GENERATOR_Y_0: U384 = U384::from_hex_unchecked("0x00690d665d446f7bd960736bcbb2efb4de03ed7274b49a58e458c282f832d204f2cf88886d8c7c2ef094094409fd4ddf");
+const GENERATOR_Y_1: U384 = U384::from_hex_unchecked("0x00f8169fd28355189e549da3151a70aa61ef11ac3d591bf12463b01acee304c24279b83f5e52270bd9a1cdd185eb8f93");
+
+/// The description of the curve.
+#[derive(Clone, Debug)]
+pub struct BLS12377TwistCurve;
+
+impl IsEllipticCurve for BLS12377TwistCurve {
+    type BaseField = Degree2ExtensionField;
+    type PointRepresentation = ShortWeierstrassProjectivePoint<Self>;
+
+    fn generator() -> Self::PointRepresentation {
+        Self::PointRepresentation::new([
+            FieldElement::new([
+                FieldElement::new(GENERATOR_X_0),
+                FieldElement::new(GENERATOR_X_1),
+            ]),
+            FieldElement::new([
+                FieldElement::new(GENERATOR_Y_0),
+                FieldElement::new(GENERATOR_Y_1),
+            ]),
+            FieldElement::one(),
+        ])
+    }
+}
+
+impl IsShortWeierstrass for BLS12377TwistCurve {
+    fn a() -> FieldElement<Self::BaseField> {
+        FieldElement::zero()
+    }
+
+    fn b() -> FieldElement<Self::BaseField> {
+        FieldElement::new([FieldElement::zero(),
+        FieldElement::from_hex_unchecked("0x10222f6db0fd6f343bd03737460c589dc7b4f91cd5fd889129207b63c6bf8000dd39e5c1ccccccd1c9ed9999999999a")])
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{
+        cyclic_group::IsGroup,
+        elliptic_curve::{
+            short_weierstrass::{
+                curves::bls12_377::field_extension::{BLS12377PrimeField, Degree2ExtensionField},
+                traits::IsShortWeierstrass,
+            },
+            traits::IsEllipticCurve,
+        },
+        field::element::FieldElement,
+        unsigned_integer::element::U384,
+    };
+
+    use super::BLS12377TwistCurve;
+    type Level0FE = FieldElement<BLS12377PrimeField>;
+    type Level1FE = FieldElement<Degree2ExtensionField>;
+
+    #[test]
+    fn create_generator() {
+        let g = BLS12377TwistCurve::generator();
+        let [x, y, _] = g.coordinates();
+        assert_eq!(
+            BLS12377TwistCurve::defining_equation(x, y),
+            Level1FE::zero()
+        );
+    }
+    #[test]
+    // Values checked in sage
+    fn add_point_three_times_equals_operate_with_self_3() {
+        let px = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x11a87eda97b96e733c2eb833ae35531b87878b416d57b370c7cd13b5f3c413387633b0ca6dfead19305318501376087")),
+            Level0FE::new(U384::from_hex_unchecked("0xa4a6d842722f2636937acf0e889ab343e121a599b8a3a9bd3be766da7d84f8e060be00f06bb2d29df963bf2d847598"))
+        ]);
+        let py = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x75589c0925d6cf45e715460ea7cb3388d4e21d9b79aa8411567d8de85ba5561bcab80a5c0b363e31817d458e5b2a2a")),
+            Level0FE::new(U384::from_hex_unchecked("0xcb4e1e4b160cc6c92e7b3dd0b2f4053bc7178d201e7788796a6035c59ccd586635796e97003b1f76eca273576f01ac"))
+        ]);
+
+        let expectedx = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x1aba32e70b88834d0cd5fb3a27eda8211eb94b6a191cd287f798145c09dc64dabda377836c31d31cc728635425ccc61")),
+            Level0FE::new(U384::from_hex_unchecked("0x146b578a9c3d92f64baafa082d27c3446fb197659bbd4ac45e290f5f49ae308599448dc288140a4049bac3c6e3e1aac"))
+        ]);
+        let expectedy = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x1417486a114a5a9074151a3a2c710dc105cce81de69a91067758355cda7049a2ad8077a2343679bba473484ad0cddd6")),
+            Level0FE::new(U384::from_hex_unchecked("0x195c939995782a07b26e3b44b49d58eb0951d452d0e8928f218a8e63f74f74860cb56265e437da80df67b6254b27cd"))
+        ]);
+        let p = BLS12377TwistCurve::create_point_from_affine(px, py).unwrap();
+        let expected = BLS12377TwistCurve::create_point_from_affine(expectedx, expectedy).unwrap();
+        assert_eq!(p.operate_with_self(3_u16), expected);
+    }
+
+    #[test]
+    // Values checked in sage
+    fn operate_with_self_test() {
+        let px = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x2233db786c8cb6a3b6846aebad4ce3f5346961c8bade4c129d920170d1ceeb02d84fd4e12b592f0cba64e083d75167")),
+            Level0FE::new(U384::from_hex_unchecked("0x12d1c7ac03cb1991993cdb9d38c50588b67c18ed9b9db5f84ac5b59c201f6493a42f690169b96b7a9f12fae4718b6cb"))
+        ]);
+
+        let py = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x160cc59bb3929b1073996aa370c880e002b81f3e4f7275636caf55754bbfcfe1d43c5d91ee7f3cb49254be2366d5d0")),
+            Level0FE::new(U384::from_hex_unchecked("0xe66460970bf2fc79831983744d7c83fad910266fd56f08b4a8f737e7609d88930503091e06228a184627b57c02352"))
+        ]);
+
+        let qx = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x1124128cf7166bb67207327079f8a3e75d876e3b6dd54cc05c766951c5d832aa202c57ed2308d78283e4c859be8fee3")),
+            Level0FE::new(U384::from_hex_unchecked("0x1889e19d4f67d120d367c15f7bc23529fe4e335627e0eb16ec2bafe6199f0e7d393c5413fc7157ec03d5d56e1eb333"))
+        ]);
+
+        let qy = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x16db05c371bf6f28e92faa77d3abf5c66c4feea84d334807f46ee69227a2144a827b00f8bcf4179b97922a85ebd5776")),
+            Level0FE::new(U384::from_hex_unchecked("0x105def68752ac5cb73f3b692b3c877f2febe6cd8a679584f4b1c64f9886f12b15dd7909251bc1e90d559fadd6b1f8f5"))
+        ]);
+
+        let scalar = U384::from_hex_unchecked(
+            "0x3061aa3679b1865fa09dac7339a87511147f015aa8997fae59b475d751bc16f",
+        );
+
+        let p = BLS12377TwistCurve::create_point_from_affine(px, py).unwrap();
+        let q = BLS12377TwistCurve::create_point_from_affine(qx, qy).unwrap();
+
+        assert_eq!(p.operate_with_self(scalar), q);
+    }
+
+    #[test]
+    // Values checked in sage
+    fn add_two_points() {
+        let px = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x11a87eda97b96e733c2eb833ae35531b87878b416d57b370c7cd13b5f3c413387633b0ca6dfead19305318501376087")),
+            Level0FE::new(U384::from_hex_unchecked("0xa4a6d842722f2636937acf0e889ab343e121a599b8a3a9bd3be766da7d84f8e060be00f06bb2d29df963bf2d847598"))
+        ]);
+
+        let py = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x75589c0925d6cf45e715460ea7cb3388d4e21d9b79aa8411567d8de85ba5561bcab80a5c0b363e31817d458e5b2a2a")),
+            Level0FE::new(U384::from_hex_unchecked("0xcb4e1e4b160cc6c92e7b3dd0b2f4053bc7178d201e7788796a6035c59ccd586635796e97003b1f76eca273576f01ac"))
+        ]);
+
+        let qx = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0xa27279ea48d8e812223956c84ae89c5de09e67c7052c056d936da4578dcf0b30799f2212af0017fe50e146c5c2ce05")),
+            Level0FE::new(U384::from_hex_unchecked("0xaed1579f3de386644e0ffac471e98f8f1d97a6be4d88a0c516bc9fd4d7780649424c3d51bb85d074d66560e2c2bb07"))
+        ]);
+
+        let qy = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x6e7e249e0d7c81d0fec62a3f96d6fefd7d6c87019559b03664a015a2ed536d98e6432b93ec219426f9729f119c7982")),
+            Level0FE::new(U384::from_hex_unchecked("0x104c4e4adb48273511a0473d742724f48042b967591ab9b86731bb902debfd44d89571af704340614f6618863fc5841"))
+        ]);
+
+        let expectedx = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x19d74f8769a27aae05c5b50b2d051849f6feab554d0fd4f51c7fe43d918fac1d966c97dd6e61b2f7fc13694495dc259")),
+            Level0FE::new(U384::from_hex_unchecked("0x11f4194886ad75e3baff9853734e9ef4fe5f20d333951ba126bca7dd54ebc4998b17d4d315d29147dc22124c5464a64"))
+        ]);
+        let expectedy = Level1FE::new([
+            Level0FE::new(U384::from_hex_unchecked("0x14faa941c0cdf5567c294bc9e73e83c4602f331a7400222c8475fba4c6a688b12ce8f7cc20e54eaac1af6046aec58bd")),
+            Level0FE::new(U384::from_hex_unchecked("0x7cee0c4760e2bf76047cda31141e6242c5893ba5c2fba5f23c0048545912e309dc647f4cfe2db17b23aa2f57c3d8c3"))
+        ]);
+
+        let p = BLS12377TwistCurve::create_point_from_affine(px, py).unwrap();
+        let q = BLS12377TwistCurve::create_point_from_affine(qx, qy).unwrap();
+        let expected = BLS12377TwistCurve::create_point_from_affine(expectedx, expectedy).unwrap();
+
+        assert_eq!(p.operate_with(&q), expected);
+    }
+}
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/field_extension.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/field_extension.rs
index ecf87e2d9..cb41e5aeb 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/field_extension.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_381/field_extension.rs
@@ -21,6 +21,7 @@ impl IsModulus<U384> for BLS12381FieldModulus {
 }
 
 pub type BLS12381PrimeField = MontgomeryBackendPrimeField<BLS12381FieldModulus, 6>;
+type Fp2E = FieldElement<Degree2ExtensionField>;
 
 //////////////////
 #[derive(Clone, Debug)]
@@ -199,6 +200,16 @@ impl HasCubicNonResidue<Degree2ExtensionField> for LevelTwoResidue {
     }
 }
 
+impl HasQuadraticNonResidue<Degree2ExtensionField> for LevelTwoResidue {
+    fn residue() -> FieldElement<Degree2ExtensionField> {
+        FieldElement::new([
+            FieldElement::new(U384::from("1")),
+            FieldElement::new(U384::from("1")),
+        ])
+    }
+}
+pub type Degree4ExtensionField = QuadraticExtensionField<Degree2ExtensionField, LevelTwoResidue>;
+
 pub type Degree6ExtensionField = CubicExtensionField<Degree2ExtensionField, LevelTwoResidue>;
 
 #[derive(Debug, Clone)]
@@ -284,6 +295,14 @@ impl FieldElement<Degree12ExtensionField> {
     }
 }
 
+/// Computes the multiplication of an element of fp2 by the level two non-residue 9+u.
+pub fn mul_fp2_by_nonresidue(a: &Fp2E) -> Fp2E {
+    // (c0 + c1 * u) * (1 + u) = (c0 - c1) + (c1 + c0) * u
+    let c0 = &a.value()[0] - &a.value()[1]; // c0 - c1
+    let c1 = &a.value()[0] + &a.value()[1]; // c1 + c0
+
+    Fp2E::new([c0, c1])
+}
 #[cfg(test)]
 mod tests {
     use crate::elliptic_curve::{
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 b0c0c6cd2..8107f69ec 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
@@ -1,11 +1,12 @@
-use super::curve::MILLER_LOOP_CONSTANT;
 use super::{
     curve::BLS12381Curve,
-    field_extension::{BLS12381PrimeField, Degree12ExtensionField, Degree2ExtensionField},
+    field_extension::{
+        mul_fp2_by_nonresidue, BLS12381PrimeField, Degree12ExtensionField, Degree2ExtensionField,
+        Degree4ExtensionField,
+    },
     twist::BLS12381TwistCurve,
 };
 use crate::{cyclic_group::IsGroup, elliptic_curve::traits::IsPairing, errors::PairingError};
-
 use crate::{
     elliptic_curve::short_weierstrass::{
         curves::bls12_381::field_extension::{Degree6ExtensionField, LevelTwoResidue},
@@ -13,11 +14,76 @@ use crate::{
         traits::IsShortWeierstrass,
     },
     field::{element::FieldElement, extensions::cubic::HasCubicNonResidue},
-    unsigned_integer::element::{UnsignedInteger, U256},
+    unsigned_integer::element::U256,
 };
 
+type FpE = FieldElement<BLS12381PrimeField>;
+type Fp2E = FieldElement<Degree2ExtensionField>;
+type Fp4E = FieldElement<Degree4ExtensionField>;
+type Fp6E = FieldElement<Degree6ExtensionField>;
+type Fp12E = FieldElement<Degree12ExtensionField>;
+
 pub const SUBGROUP_ORDER: U256 =
     U256::from_hex_unchecked("73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001");
+// value of x in binary
+
+// We use |x|, then if needed we apply the minus sign in the final exponentiation by applying the inverse (in this case the conjugate because the element is in the cyclotomic subgroup)
+pub const X: u64 = 0xd201000000010000;
+
+// X = 1101001000000001000000000000000000000000000000010000000000000000
+pub const X_BINARY: &[bool] = &[
+    true, true, false, true, false, false, true, false, false, false, false, false, false, false,
+    false, true, false, false, false, false, false, false, false, false, false, false, false,
+    false, false, false, false, false, false, false, false, false, false, false, false, false,
+    false, false, false, false, false, false, false, true, false, false, false, false, false,
+    false, false, false, false, false, false, false, false, false, false, false,
+];
+
+// GAMMA constants used to compute the Frobenius morphisms
+/// We took these constants from https://github.com/hecmas/zkNotebook/blob/main/src/BLS12381/constants.ts
+
+pub const GAMMA_11: Fp2E = Fp2E::const_from_raw([
+    FpE::from_hex_unchecked("1904D3BF02BB0667C231BEB4202C0D1F0FD603FD3CBD5F4F7B2443D784BAB9C4F67EA53D63E7813D8D0775ED92235FB8"),
+    FpE::from_hex_unchecked("FC3E2B36C4E03288E9E902231F9FB854A14787B6C7B36FEC0C8EC971F63C5F282D5AC14D6C7EC22CF78A126DDC4AF3"),
+]);
+
+pub const GAMMA_12: Fp2E = Fp2E::const_from_raw([
+    FpE::from_hex_unchecked("0"),
+    FpE::from_hex_unchecked("1A0111EA397FE699EC02408663D4DE85AA0D857D89759AD4897D29650FB85F9B409427EB4F49FFFD8BFD00000000AAAC"),
+]);
+
+pub const GAMMA_13: Fp2E = Fp2E::const_from_raw([
+    FpE::from_hex_unchecked("6AF0E0437FF400B6831E36D6BD17FFE48395DABC2D3435E77F76E17009241C5EE67992F72EC05F4C81084FBEDE3CC09"),
+    FpE::from_hex_unchecked("6AF0E0437FF400B6831E36D6BD17FFE48395DABC2D3435E77F76E17009241C5EE67992F72EC05F4C81084FBEDE3CC09"),
+]);
+
+pub const GAMMA_14: Fp2E = Fp2E::const_from_raw([
+    FpE::from_hex_unchecked("1A0111EA397FE699EC02408663D4DE85AA0D857D89759AD4897D29650FB85F9B409427EB4F49FFFD8BFD00000000AAAD"),
+    FpE::from_hex_unchecked("0"),
+]);
+
+pub const GAMMA_15: Fp2E = Fp2E::const_from_raw([
+    FpE::from_hex_unchecked("5B2CFD9013A5FD8DF47FA6B48B1E045F39816240C0B8FEE8BEADF4D8E9C0566C63A3E6E257F87329B18FAE980078116"),
+    FpE::from_hex_unchecked("144E4211384586C16BD3AD4AFA99CC9170DF3560E77982D0DB45F3536814F0BD5871C1908BD478CD1EE605167FF82995"),
+]);
+
+/// GAMMA_2i = GAMMA_1i * GAMMA_1i.conjugate()
+pub const GAMMA_21: FpE = FpE::from_hex_unchecked(
+    "5F19672FDF76CE51BA69C6076A0F77EADDB3A93BE6F89688DE17D813620A00022E01FFFFFFFEFFFF",
+);
+
+pub const GAMMA_22: FpE = FpE::from_hex_unchecked(
+    "5F19672FDF76CE51BA69C6076A0F77EADDB3A93BE6F89688DE17D813620A00022E01FFFFFFFEFFFE",
+);
+
+pub const GAMMA_23: FpE =
+    FpE::from_hex_unchecked("1A0111EA397FE69A4B1BA7B6434BACD764774B84F38512BF6730D2A0F6B0F6241EABFFFEB153FFFFB9FEFFFFFFFFAAAA");
+
+pub const GAMMA_24: FpE =
+    FpE::from_hex_unchecked("1A0111EA397FE699EC02408663D4DE85AA0D857D89759AD4897D29650FB85F9B409427EB4F49FFFD8BFD00000000AAAC");
+
+pub const GAMMA_25: FpE =
+    FpE::from_hex_unchecked("1A0111EA397FE699EC02408663D4DE85AA0D857D89759AD4897D29650FB85F9B409427EB4F49FFFD8BFD00000000AAAD");
 
 #[derive(Clone)]
 pub struct BLS12381AtePairing;
@@ -46,6 +112,26 @@ impl IsPairing for BLS12381AtePairing {
     }
 }
 
+/// Implements the miller loop for the ate pairing of the BLS12 381 curve.
+/// Based on algorithm 9.2, page 212 of the book
+/// "Topics in computational number theory" by W. Bons and K. Lenstra
+#[allow(unused)]
+pub fn miller(
+    q: &ShortWeierstrassProjectivePoint<BLS12381TwistCurve>,
+    p: &ShortWeierstrassProjectivePoint<BLS12381Curve>,
+) -> FieldElement<Degree12ExtensionField> {
+    let mut r = q.clone();
+    let mut f = FieldElement::<Degree12ExtensionField>::one();
+    X_BINARY.iter().skip(1).for_each(|bit| {
+        double_accumulate_line(&mut r, p, &mut f);
+        if *bit {
+            add_accumulate_line(&mut r, q, p, &mut f);
+        }
+    });
+
+    f.conjugate()
+}
+
 fn double_accumulate_line(
     t: &mut ShortWeierstrassProjectivePoint<BLS12381TwistCurve>,
     p: &ShortWeierstrassProjectivePoint<BLS12381Curve>,
@@ -152,69 +238,148 @@ fn add_accumulate_line(
         ]),
     ]);
 }
-/// Implements the miller loop for the ate pairing of the BLS12 381 curve.
-/// Based on algorithm 9.2, page 212 of the book
-/// "Topics in computational number theory" by W. Bons and K. Lenstra
-#[allow(unused)]
-fn miller(
-    q: &ShortWeierstrassProjectivePoint<BLS12381TwistCurve>,
-    p: &ShortWeierstrassProjectivePoint<BLS12381Curve>,
-) -> FieldElement<Degree12ExtensionField> {
-    let mut r = q.clone();
-    let mut f = FieldElement::<Degree12ExtensionField>::one();
-    let mut miller_loop_constant = MILLER_LOOP_CONSTANT;
-    let mut miller_loop_constant_bits: alloc::vec::Vec<bool> = alloc::vec![];
 
-    while miller_loop_constant > 0 {
-        miller_loop_constant_bits.insert(0, (miller_loop_constant & 1) == 1);
-        miller_loop_constant >>= 1;
-    }
+// To understand more about how to reduce the final exponentiation
+// read "Efficient Final Exponentiation via Cyclotomic Structure for
+// Pairings over Families of Elliptic Curves" (https://eprint.iacr.org/2020/875.pdf)
+// Hard part from https://eprint.iacr.org/2020/875.pdf p14
+// Same implementation as in Constantine https://github.com/mratsim/constantine/blob/master/constantine/math/pairings/pairings_bls12.nim
+pub fn final_exponentiation(f: &Fp12E) -> Fp12E {
+    let f_easy_aux = f.conjugate() * f.inv().unwrap();
+    let mut f_easy = frobenius_square(&f_easy_aux) * &f_easy_aux;
+
+    let mut v2 = cyclotomic_square(&f_easy); // v2 = f²
+    let mut v0 = cyclotomic_pow_x(&f_easy).conjugate(); //  v0 = f^x
+    let mut v1 = f_easy.conjugate(); // v1 = f^-1
+
+    //  (x−1)²
+    v0 *= v1; // v0 = f^(x-1)
+    v1 = cyclotomic_pow_x(&v0).conjugate(); // v1 = (f^(x-1))^(x)
+
+    v0 = v0.conjugate(); // v0 = (f^(x-1))^(-1)
+    v0 *= &v1; // v0 = (f^(x-1))^(-1) * (f^(x-1))^x = (f^(x-1))^(x-1) =  f^((x-1)²)
+
+    // (x+p)
+    v1 = cyclotomic_pow_x(&v0).conjugate(); // v1 = f^((x-1)².x)
+    v0 = frobenius(&v0); // f^((x-1)².p)
+    v0 *= &v1; // f^((x-1)².p + (x-1)².x) = f^((x-1)².(x+p))
+
+    // + 3
+    f_easy *= v2; // f^3
+
+    // (x²+p²−1)
+    v2 = cyclotomic_pow_x(&v0).conjugate();
+    v1 = cyclotomic_pow_x(&v2).conjugate(); // v1 = f^((x-1)².(x+p).x²)
+    v2 = frobenius_square(&v0); // v2 = f^((x-1)².(x+p).p²)
+    v0 = v0.conjugate(); // v0 = f^((x-1)².(x+p).-1)
+    v0 *= &v1; // v0 = f^((x-1)².(x+p).(x²-1))
+    v0 *= &v2; // v0 = f^((x-1)².(x+p).(x²+p²-1))
+
+    f_easy *= &v0;
+    f_easy
+}
 
-    for bit in miller_loop_constant_bits[1..].iter() {
-        double_accumulate_line(&mut r, p, &mut f);
-        if *bit {
-            add_accumulate_line(&mut r, q, p, &mut f);
-        }
-    }
-    f.inv().unwrap()
+////////////////// FROBENIUS MORPHISIMS //////////////////
+pub fn frobenius(f: &Fp12E) -> Fp12E {
+    let [a, b] = f.value(); // f = a + bw, where a and b in Fp6.
+    let [a0, a1, a2] = a.value(); // a = a0 + a1 * v + a2 * v^2, where a0, a1 and a2 in Fp2.
+    let [b0, b1, b2] = b.value(); // b = b0 + b1 * v + b2 * v^2, where b0, b1 and b2 in Fp2.
+
+    // c1 = a0.conjugate() + a1.conjugate() * GAMMA_12 * v + a2.conjugate() * GAMMA_14 * v^2
+    let c1 = Fp6E::new([
+        a0.conjugate(),
+        a1.conjugate() * GAMMA_12,
+        a2.conjugate() * GAMMA_14,
+    ]);
+
+    let c2 = Fp6E::new([
+        b0.conjugate() * GAMMA_11,
+        b1.conjugate() * GAMMA_13,
+        b2.conjugate() * GAMMA_15,
+    ]);
+
+    Fp12E::new([c1, c2]) //c1 + c2 * w
 }
 
-/// Auxiliary function for the final exponentiation of the ate pairing.
 fn frobenius_square(
     f: &FieldElement<Degree12ExtensionField>,
 ) -> FieldElement<Degree12ExtensionField> {
     let [a, b] = f.value();
-    let w_raised_to_p_squared_minus_one = FieldElement::<Degree6ExtensionField>::new_base("1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad");
-    let omega_3 = FieldElement::<Degree2ExtensionField>::new_base("1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac");
-    let omega_3_squared = FieldElement::<Degree2ExtensionField>::new_base(
-        "5f19672fdf76ce51ba69c6076a0f77eaddb3a93be6f89688de17d813620a00022e01fffffffefffe",
-    );
-
     let [a0, a1, a2] = a.value();
     let [b0, b1, b2] = b.value();
 
-    let f0 = FieldElement::new([a0.clone(), a1 * &omega_3, a2 * &omega_3_squared]);
-    let f1 = FieldElement::new([b0.clone(), b1 * omega_3, b2 * omega_3_squared]);
+    let c1 = Fp6E::new([a0.clone(), GAMMA_22 * a1, GAMMA_24 * a2]);
+    let c2 = Fp6E::new([GAMMA_21 * b0, GAMMA_23 * b1, GAMMA_25 * b2]);
 
-    FieldElement::new([f0, w_raised_to_p_squared_minus_one * f1])
+    Fp12E::new([c1, c2])
 }
 
-// To understand more about how to reduce the final exponentiation
-// read "Efficient Final Exponentiation via Cyclotomic Structure for
-// Pairings over Families of Elliptic Curves" (https://eprint.iacr.org/2020/875.pdf)
-//
-// TODO: implement optimizations for the hard part of the final exponentiation.
-#[allow(unused)]
-fn final_exponentiation(
-    base: &FieldElement<Degree12ExtensionField>,
-) -> FieldElement<Degree12ExtensionField> {
-    const PHI_DIVIDED_BY_R: UnsignedInteger<20> = UnsignedInteger::from_hex_unchecked("f686b3d807d01c0bd38c3195c899ed3cde88eeb996ca394506632528d6a9a2f230063cf081517f68f7764c28b6f8ae5a72bce8d63cb9f827eca0ba621315b2076995003fc77a17988f8761bdc51dc2378b9039096d1b767f17fcbde783765915c97f36c6f18212ed0b283ed237db421d160aeb6a1e79983774940996754c8c71a2629b0dea236905ce937335d5b68fa9912aae208ccf1e516c3f438e3ba79");
-
-    let f1 = base.conjugate() * base.inv().unwrap();
-    let f2 = frobenius_square(&f1) * f1;
-    f2.pow(PHI_DIVIDED_BY_R)
+////////////////// CYCLOTOMIC SUBGROUP OPERATIONS //////////////////
+/// Since the result of the Easy Part of the Final Exponentiation belongs to the cyclotomic
+/// subgroup of Fp12, we can optimize the square and pow operations used in the Hard Part.
+
+/// Computes the square of an element of a cyclotomic subgroup of Fp12.
+/// Algorithm from Constantine's cyclotomic_square_quad_over_cube
+/// https://github.com/mratsim/constantine/blob/master/constantine/math/pairings/cyclotomic_subgroups.nim#L354
+pub fn cyclotomic_square(a: &Fp12E) -> Fp12E {
+    // a = g + h * w
+    let [g, h] = a.value();
+    let [b0, b1, b2] = g.value();
+    let [b3, b4, b5] = h.value();
+
+    let v0 = Fp4E::new([b0.clone(), b4.clone()]).square();
+    let v1 = Fp4E::new([b3.clone(), b2.clone()]).square();
+    let v2 = Fp4E::new([b1.clone(), b5.clone()]).square();
+
+    // r = r0 + r1 * w
+    // r0 = r00 + r01 * v + r02 * v^2
+    // r1 = r10 + r11 * v + r12 * v^2
+
+    // r00 = 3v00 - 2b0
+    let mut r00 = &v0.value()[0] - b0;
+    r00 = r00.double();
+    r00 += v0.value()[0].clone();
+
+    // r01 = 3v10 -2b1
+    let mut r01 = &v1.value()[0] - b1;
+    r01 = r01.double();
+    r01 += v1.value()[0].clone();
+
+    // r11 = 3v01 - 2b4
+    let mut r11 = &v0.value()[1] + b4;
+    r11 = r11.double();
+    r11 += v0.value()[1].clone();
+
+    // r12 = 3v11 - 2b5
+    let mut r12 = &v1.value()[1] + b5;
+    r12 = r12.double();
+    r12 += v1.value()[1].clone();
+
+    // 3 * (9 + u) * v21 + 2b3
+    let v21 = mul_fp2_by_nonresidue(&v2.value()[1]);
+    let mut r10 = &v21 + b3;
+    r10 = r10.double();
+    r10 += v21;
+
+    // 3 * (9 + u) * v20 - 2b3
+    let mut r02 = &v2.value()[0] - b2;
+    r02 = r02.double();
+    r02 += v2.value()[0].clone();
+
+    Fp12E::new([Fp6E::new([r00, r01, r02]), Fp6E::new([r10, r11, r12])])
 }
 
+#[allow(clippy::needless_range_loop)]
+pub fn cyclotomic_pow_x(f: &Fp12E) -> Fp12E {
+    let mut result = Fp12E::one();
+    X_BINARY.iter().for_each(|&bit| {
+        result = cyclotomic_square(&result);
+        if bit {
+            result = &result * f;
+        }
+    });
+    result
+}
 #[cfg(test)]
 mod tests {
     use crate::{
@@ -294,4 +459,58 @@ mod tests {
         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(&[
+            "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
+        ]);
+        let mut result = frobenius(&f);
+        for _ in 1..12 {
+            result = frobenius(&result);
+        }
+        assert_eq!(f, result)
+    }
+
+    #[test]
+    fn apply_6_times_frobenius_square_is_identity() {
+        let f = Fp12E::from_coefficients(&[
+            "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
+        ]);
+        let mut result = frobenius_square(&f);
+        for _ in 1..6 {
+            result = frobenius_square(&result);
+        }
+        assert_eq!(f, result)
+    }
+
+    #[test]
+    fn cyclotomic_square_equals_square() {
+        let p = BLS12381Curve::generator();
+        let q = BLS12381TwistCurve::generator();
+        let f = miller(&q, &p);
+        let f_easy_aux = f.conjugate() * f.inv().unwrap(); // f ^ (p^6 - 1) because f^(p^6) = f.conjugate().
+        let f_easy = &frobenius_square(&f_easy_aux) * f_easy_aux; // (f^{p^6 - 1})^(p^2) * (f^{p^6 - 1}).
+        assert_eq!(cyclotomic_square(&f_easy), f_easy.square());
+    }
+
+    #[test]
+    fn test_double_accumulate_line_doubles_point_correctl_2() {
+        let g1 = BLS12381Curve::generator();
+        let g2 = BLS12381TwistCurve::generator();
+        let mut r = g2.clone();
+        let mut f = FieldElement::one();
+        double_accumulate_line(&mut r, &g1, &mut f);
+        let expected_r = g2.operate_with(&g2);
+        assert_eq!(r.to_affine(), expected_r.to_affine());
+    }
+
+    #[test]
+    fn cyclotomic_pow_x_equals_pow() {
+        let p = BLS12381Curve::generator();
+        let q = BLS12381TwistCurve::generator();
+        let f = miller(&q, &p);
+        let f_easy_aux = f.conjugate() * f.inv().unwrap(); // f ^ (p^6 - 1) because f^(p^6) = f.conjugate().
+        let f_easy = &frobenius_square(&f_easy_aux) * f_easy_aux; // (f^{p^6 - 1})^(p^2) * (f^{p^6 - 1}).
+        assert_eq!(cyclotomic_pow_x(&f_easy), f_easy.pow(X));
+    }
 }
diff --git a/math/src/elliptic_curve/short_weierstrass/point.rs b/math/src/elliptic_curve/short_weierstrass/point.rs
index 888d3f50a..679edeaea 100644
--- a/math/src/elliptic_curve/short_weierstrass/point.rs
+++ b/math/src/elliptic_curve/short_weierstrass/point.rs
@@ -52,6 +52,9 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
     }
 
     pub fn double(&self) -> Self {
+        if self.is_neutral_element() {
+            return self.clone();
+        }
         let [px, py, pz] = self.coordinates();
 
         let px_square = px * px;
@@ -86,11 +89,6 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
     }
     // https://hyperelliptic.org/EFD/g1p/data/shortw/projective/addition/madd-1998-cmo
     pub fn operate_with_affine(&self, other: &Self) -> Self {
-        let [px, py, pz] = self.coordinates();
-        let [qx, qy, _qz] = other.coordinates();
-        let u = qy * pz;
-        let v = qx * pz;
-
         if self.is_neutral_element() {
             return other.clone();
         }
@@ -98,6 +96,11 @@ impl<E: IsShortWeierstrass> ShortWeierstrassProjectivePoint<E> {
             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() {
                 return Self::new([
diff --git a/math/src/field/fields/mersenne31/extension.rs b/math/src/field/fields/mersenne31/extension.rs
deleted file mode 100644
index 3c89a2147..000000000
--- a/math/src/field/fields/mersenne31/extension.rs
+++ /dev/null
@@ -1,304 +0,0 @@
-use crate::field::{
-    element::FieldElement,
-    errors::FieldError,
-    extensions::{
-        cubic::{CubicExtensionField, HasCubicNonResidue},
-        quadratic::{HasQuadraticNonResidue, QuadraticExtensionField},
-    },
-    traits::IsField,
-};
-
-use super::field::Mersenne31Field;
-
-//Note: The inverse calculation in mersenne31/plonky3 differs from the default quadratic extension so I implemented the complex extension.
-//////////////////
-#[derive(Clone, Debug)]
-pub struct Mersenne31Complex;
-
-impl IsField for Mersenne31Complex {
-    //Elements represents a[0] = real, a[1] = imaginary
-    type BaseType = [FieldElement<Mersenne31Field>; 2];
-
-    /// Returns the component wise addition of `a` and `b`
-    fn add(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        [a[0] + b[0], a[1] + b[1]]
-    }
-
-    //NOTE: THIS uses Gauss algorithm. Bench this against plonky 3 implementation to see what is faster.
-    /// Returns the multiplication of `a` and `b` using the following
-    /// equation:
-    /// (a0 + a1 * t) * (b0 + b1 * t) = a0 * b0 + a1 * b1 * Self::residue() + (a0 * b1 + a1 * b0) * t
-    /// where `t.pow(2)` equals `Q::residue()`.
-    fn mul(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        let a0b0 = a[0] * b[0];
-        let a1b1 = a[1] * b[1];
-        let z = (a[0] + a[1]) * (b[0] + b[1]);
-        [a0b0 - a1b1, z - a0b0 - a1b1]
-    }
-
-    fn square(a: &Self::BaseType) -> Self::BaseType {
-        let [a0, a1] = a;
-        let v0 = a0 * a1;
-        let c0 = (a0 + a1) * (a0 - a1);
-        let c1 = v0 + v0;
-        [c0, c1]
-    }
-    /// Returns the component wise subtraction of `a` and `b`
-    fn sub(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        [a[0] - b[0], a[1] - b[1]]
-    }
-
-    /// Returns the component wise negation of `a`
-    fn neg(a: &Self::BaseType) -> Self::BaseType {
-        [-a[0], -a[1]]
-    }
-
-    /// Returns the multiplicative inverse of `a`
-    fn inv(a: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
-        let inv_norm = (a[0].pow(2_u64) + a[1].pow(2_u64)).inv()?;
-        Ok([a[0] * inv_norm, -a[1] * inv_norm])
-    }
-
-    /// Returns the division of `a` and `b`
-    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        Self::mul(a, &Self::inv(b).unwrap())
-    }
-
-    /// Returns a boolean indicating whether `a` and `b` are equal component wise.
-    fn eq(a: &Self::BaseType, b: &Self::BaseType) -> bool {
-        a[0] == b[0] && a[1] == b[1]
-    }
-
-    /// Returns the additive neutral element of the field extension.
-    fn zero() -> Self::BaseType {
-        [FieldElement::zero(), FieldElement::zero()]
-    }
-
-    /// Returns the multiplicative neutral element of the field extension.
-    fn one() -> Self::BaseType {
-        [FieldElement::one(), FieldElement::zero()]
-    }
-
-    /// Returns the element `x * 1` where 1 is the multiplicative neutral element.
-    fn from_u64(x: u64) -> Self::BaseType {
-        [FieldElement::from(x), FieldElement::zero()]
-    }
-
-    /// Takes as input an element of BaseType and returns the internal representation
-    /// of that element in the field.
-    /// Note: for this case this is simply the identity, because the components
-    /// already have correct representations.
-    fn from_base_type(x: Self::BaseType) -> Self::BaseType {
-        x
-    }
-}
-
-pub type Mersenne31ComplexQuadraticExtensionField =
-    QuadraticExtensionField<Mersenne31Field, Mersenne31Complex>;
-
-//TODO: Check this should be for complex and not base field
-impl HasQuadraticNonResidue<Mersenne31Complex> for Mersenne31Complex {
-    // Verifiable in Sage with
-    // ```sage
-    // p = 2**31 - 1  # Mersenne31
-    // F = GF(p)  # The base field GF(p)
-    // R.<x> = F[]  # The polynomial ring over F
-    // K.<i> = F.extension(x^2 + 1)  # The complex extension field
-    // R2.<y> = K[]
-    // f2 = y^2 - i - 2
-    // assert f2.is_irreducible()
-    // ```
-    fn residue() -> FieldElement<Mersenne31Complex> {
-        FieldElement::from(&Mersenne31Complex::from_base_type([
-            FieldElement::<Mersenne31Field>::from(2),
-            FieldElement::<Mersenne31Field>::one(),
-        ]))
-    }
-}
-
-pub type Mersenne31ComplexCubicExtensionField =
-    CubicExtensionField<Mersenne31Field, Mersenne31Complex>;
-
-impl HasCubicNonResidue<Mersenne31Complex> for Mersenne31Complex {
-    // Verifiable in Sage with
-    // ```sage
-    // p = 2**31 - 1  # Mersenne31
-    // F = GF(p)  # The base field GF(p)
-    // R.<x> = F[]  # The polynomial ring over F
-    // K.<i> = F.extension(x^2 + 1)  # The complex extension field
-    // R2.<y> = K[]
-    // f2 = y^3 - 5*i
-    // assert f2.is_irreducible()
-    // ```
-    fn residue() -> FieldElement<Mersenne31Complex> {
-        FieldElement::from(&Mersenne31Complex::from_base_type([
-            FieldElement::<Mersenne31Field>::zero(),
-            FieldElement::<Mersenne31Field>::from(5),
-        ]))
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use crate::field::fields::mersenne31::field::MERSENNE_31_PRIME_FIELD_ORDER;
-
-    use super::*;
-
-    type Fi = Mersenne31Complex;
-    type F = FieldElement<Mersenne31Field>;
-
-    //NOTE: from_u64 reflects from_real
-    //NOTE: for imag use from_base_type
-
-    #[test]
-    fn add_real_one_plus_one_is_two() {
-        assert_eq!(Fi::add(&Fi::one(), &Fi::one()), Fi::from_u64(2))
-    }
-
-    #[test]
-    fn add_real_neg_one_plus_one_is_zero() {
-        assert_eq!(Fi::add(&Fi::neg(&Fi::one()), &Fi::one()), Fi::zero())
-    }
-
-    #[test]
-    fn add_real_neg_one_plus_two_is_one() {
-        assert_eq!(Fi::add(&Fi::neg(&Fi::one()), &Fi::from_u64(2)), Fi::one())
-    }
-
-    #[test]
-    fn add_real_neg_one_plus_neg_one_is_order_sub_two() {
-        assert_eq!(
-            Fi::add(&Fi::neg(&Fi::one()), &Fi::neg(&Fi::one())),
-            Fi::from_u64((MERSENNE_31_PRIME_FIELD_ORDER - 2).into())
-        )
-    }
-
-    #[test]
-    fn add_complex_one_plus_one_two() {
-        //Manually declare the complex part to one
-        let one = Fi::from_base_type([F::zero(), F::one()]);
-        let two = Fi::from_base_type([F::zero(), F::from(2)]);
-        assert_eq!(Fi::add(&one, &one), two)
-    }
-
-    #[test]
-    fn add_complex_neg_one_plus_one_is_zero() {
-        //Manually declare the complex part to one
-        let neg_one = Fi::from_base_type([F::zero(), -F::one()]);
-        let one = Fi::from_base_type([F::zero(), F::one()]);
-        assert_eq!(Fi::add(&neg_one, &one), Fi::zero())
-    }
-
-    #[test]
-    fn add_complex_neg_one_plus_two_is_one() {
-        let neg_one = Fi::from_base_type([F::zero(), -F::one()]);
-        let two = Fi::from_base_type([F::zero(), F::from(2)]);
-        let one = Fi::from_base_type([F::zero(), F::one()]);
-        assert_eq!(Fi::add(&neg_one, &two), one)
-    }
-
-    #[test]
-    fn add_complex_neg_one_plus_neg_one_imag_is_order_sub_two() {
-        let neg_one = Fi::from_base_type([F::zero(), -F::one()]);
-        assert_eq!(
-            Fi::add(&neg_one, &neg_one)[1],
-            F::new(MERSENNE_31_PRIME_FIELD_ORDER - 2)
-        )
-    }
-
-    #[test]
-    fn add_order() {
-        let a = Fi::from_base_type([-F::one(), F::one()]);
-        let b = Fi::from_base_type([F::from(2), F::new(MERSENNE_31_PRIME_FIELD_ORDER - 2)]);
-        let c = Fi::from_base_type([F::one(), -F::one()]);
-        assert_eq!(Fi::add(&a, &b), c)
-    }
-
-    #[test]
-    fn add_equal_zero() {
-        let a = Fi::from_base_type([-F::one(), -F::one()]);
-        let b = Fi::from_base_type([F::one(), F::one()]);
-        assert_eq!(Fi::add(&a, &b), Fi::zero())
-    }
-
-    #[test]
-    fn add_plus_one() {
-        let a = Fi::from_base_type([F::one(), F::from(2)]);
-        let b = Fi::from_base_type([F::one(), F::one()]);
-        let c = Fi::from_base_type([F::from(2), F::from(3)]);
-        assert_eq!(Fi::add(&a, &b), c)
-    }
-
-    #[test]
-    fn sub_real_one_sub_one_is_zero() {
-        assert_eq!(Fi::sub(&Fi::one(), &Fi::one()), Fi::zero())
-    }
-
-    #[test]
-    fn sub_real_two_sub_two_is_zero() {
-        assert_eq!(
-            Fi::sub(&Fi::from_u64(2u64), &Fi::from_u64(2u64)),
-            Fi::zero()
-        )
-    }
-
-    #[test]
-    fn sub_real_neg_one_sub_neg_one_is_zero() {
-        assert_eq!(
-            Fi::sub(&Fi::neg(&Fi::one()), &Fi::neg(&Fi::one())),
-            Fi::zero()
-        )
-    }
-
-    #[test]
-    fn sub_real_two_sub_one_is_one() {
-        assert_eq!(Fi::sub(&Fi::from_u64(2), &Fi::one()), Fi::one())
-    }
-
-    #[test]
-    fn sub_real_neg_one_sub_zero_is_neg_one() {
-        assert_eq!(
-            Fi::sub(&Fi::neg(&Fi::one()), &Fi::zero()),
-            Fi::neg(&Fi::one())
-        )
-    }
-
-    #[test]
-    fn sub_complex_one_sub_one_is_zero() {
-        let one = Fi::from_base_type([F::zero(), F::one()]);
-        assert_eq!(Fi::sub(&one, &one), Fi::zero())
-    }
-
-    #[test]
-    fn sub_complex_two_sub_two_is_zero() {
-        let two = Fi::from_base_type([F::zero(), F::from(2)]);
-        assert_eq!(Fi::sub(&two, &two), Fi::zero())
-    }
-
-    #[test]
-    fn sub_complex_neg_one_sub_neg_one_is_zero() {
-        let neg_one = Fi::from_base_type([F::zero(), -F::one()]);
-        assert_eq!(Fi::sub(&neg_one, &neg_one), Fi::zero())
-    }
-
-    #[test]
-    fn sub_complex_two_sub_one_is_one() {
-        let two = Fi::from_base_type([F::zero(), F::from(2)]);
-        let one = Fi::from_base_type([F::zero(), F::one()]);
-        assert_eq!(Fi::sub(&two, &one), one)
-    }
-
-    #[test]
-    fn sub_complex_neg_one_sub_zero_is_neg_one() {
-        let neg_one = Fi::from_base_type([F::zero(), -F::one()]);
-        assert_eq!(Fi::sub(&neg_one, &Fi::zero()), neg_one)
-    }
-
-    #[test]
-    fn mul() {
-        let a = Fi::from_base_type([F::from(2), F::from(2)]);
-        let b = Fi::from_base_type([F::from(4), F::from(5)]);
-        let c = Fi::from_base_type([-F::from(2), F::from(18)]);
-        assert_eq!(Fi::mul(&a, &b), c)
-    }
-}
diff --git a/math/src/field/fields/mersenne31/extensions.rs b/math/src/field/fields/mersenne31/extensions.rs
new file mode 100644
index 000000000..27c2ab118
--- /dev/null
+++ b/math/src/field/fields/mersenne31/extensions.rs
@@ -0,0 +1,627 @@
+use super::field::Mersenne31Field;
+use crate::field::{
+    element::FieldElement,
+    errors::FieldError,
+    traits::{IsField, IsSubFieldOf},
+};
+#[cfg(feature = "alloc")]
+use alloc::vec::Vec;
+
+type FpE = FieldElement<Mersenne31Field>;
+
+#[derive(Clone, Debug)]
+pub struct Degree2ExtensionField;
+
+impl Degree2ExtensionField {
+    pub fn mul_fp2_by_nonresidue(a: &Fp2E) -> Fp2E {
+        Fp2E::new([
+            a.value()[0].double() - a.value()[1],
+            a.value()[1].double() + a.value()[0],
+        ])
+    }
+}
+
+impl IsField for Degree2ExtensionField {
+    //Element representation: a[0] = real part, a[1] = imaginary part
+    type BaseType = [FpE; 2];
+
+    /// Returns the component wise addition of `a` and `b`
+    fn add(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        [a[0] + b[0], a[1] + b[1]]
+    }
+
+    /// Returns the multiplication of `a` and `b`.
+    fn mul(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        let a0b0 = a[0] * b[0];
+        let a1b1 = a[1] * b[1];
+        let z = (a[0] + a[1]) * (b[0] + b[1]);
+        [a0b0 - a1b1, z - a0b0 - a1b1]
+    }
+
+    fn square(a: &Self::BaseType) -> Self::BaseType {
+        let [a0, a1] = a;
+        let v0 = a0 * a1;
+        let c0 = (a0 + a1) * (a0 - a1);
+        let c1 = v0.double();
+        [c0, c1]
+    }
+    /// Returns the component wise subtraction of `a` and `b`
+    fn sub(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        [a[0] - b[0], a[1] - b[1]]
+    }
+
+    /// Returns the component wise negation of `a`
+    fn neg(a: &Self::BaseType) -> Self::BaseType {
+        [-a[0], -a[1]]
+    }
+
+    /// Returns the multiplicative inverse of `a`
+    fn inv(a: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
+        let inv_norm = (a[0].square() + a[1].square()).inv()?;
+        Ok([a[0] * inv_norm, -a[1] * inv_norm])
+    }
+
+    /// Returns the division of `a` and `b`
+    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    }
+
+    /// Returns a boolean indicating whether `a` and `b` are equal component wise.
+    fn eq(a: &Self::BaseType, b: &Self::BaseType) -> bool {
+        a[0] == b[0] && a[1] == b[1]
+    }
+
+    /// Returns the multiplicative neutral element of the field extension.
+    fn one() -> Self::BaseType {
+        [FpE::one(), FpE::zero()]
+    }
+
+    /// Returns the element `x * 1` where 1 is the multiplicative neutral element.
+    fn from_u64(x: u64) -> Self::BaseType {
+        [FpE::from(x), FpE::zero()]
+    }
+
+    /// Takes as input an element of BaseType and returns the internal representation
+    /// of that element in the field.
+    /// Note: for this case this is simply the identity, because the components
+    /// already have correct representations.
+    fn from_base_type(x: Self::BaseType) -> Self::BaseType {
+        x
+    }
+}
+
+impl IsSubFieldOf<Degree2ExtensionField> for Mersenne31Field {
+    fn add(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        [FpE::from(a) + b[0], b[1]]
+    }
+
+    fn sub(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        [FpE::from(a) - b[0], -b[1]]
+    }
+
+    fn mul(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        [FpE::from(a) * b[0], FpE::from(a) * b[1]]
+    }
+
+    fn div(
+        a: &Self::BaseType,
+        b: &<Degree2ExtensionField as IsField>::BaseType,
+    ) -> <Degree2ExtensionField as IsField>::BaseType {
+        let b_inv = Degree2ExtensionField::inv(b).unwrap();
+        <Self as IsSubFieldOf<Degree2ExtensionField>>::mul(a, &b_inv)
+    }
+
+    fn embed(a: Self::BaseType) -> <Degree2ExtensionField as IsField>::BaseType {
+        [FieldElement::from_raw(a), FieldElement::zero()]
+    }
+
+    #[cfg(feature = "alloc")]
+    fn to_subfield_vec(
+        b: <Degree2ExtensionField as IsField>::BaseType,
+    ) -> alloc::vec::Vec<Self::BaseType> {
+        b.into_iter().map(|x| x.to_raw()).collect()
+    }
+}
+
+type Fp2E = FieldElement<Degree2ExtensionField>;
+
+#[derive(Clone, Debug)]
+pub struct Degree4ExtensionField;
+
+impl IsField for Degree4ExtensionField {
+    type BaseType = [Fp2E; 2];
+
+    fn add(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        [&a[0] + &b[0], &a[1] + &b[1]]
+    }
+
+    fn sub(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        [&a[0] - &b[0], &a[1] - &b[1]]
+    }
+
+    fn neg(a: &Self::BaseType) -> Self::BaseType {
+        [-&a[0], -&a[1]]
+    }
+
+    fn mul(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        // Algorithm from: https://github.com/ingonyama-zk/papers/blob/main/Mersenne31_polynomial_arithmetic.pdf (page 5):
+        let a0b0 = &a[0] * &b[0];
+        let a1b1 = &a[1] * &b[1];
+        [
+            &a0b0 + Degree2ExtensionField::mul_fp2_by_nonresidue(&a1b1),
+            (&a[0] + &a[1]) * (&b[0] + &b[1]) - a0b0 - a1b1,
+        ]
+    }
+
+    fn square(a: &Self::BaseType) -> Self::BaseType {
+        let a0_square = &a[0].square();
+        let a1_square = &a[1].square();
+        [
+            a0_square + Degree2ExtensionField::mul_fp2_by_nonresidue(a1_square),
+            (&a[0] + &a[1]).square() - a0_square - a1_square,
+        ]
+    }
+
+    fn inv(a: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
+        let inv_norm =
+            (a[0].square() - Degree2ExtensionField::mul_fp2_by_nonresidue(&a[1].square())).inv()?;
+        Ok([&a[0] * &inv_norm, -&a[1] * &inv_norm])
+    }
+
+    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
+        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    }
+
+    fn eq(a: &Self::BaseType, b: &Self::BaseType) -> bool {
+        a[0] == b[0] && a[1] == b[1]
+    }
+
+    fn zero() -> Self::BaseType {
+        [Fp2E::zero(), Fp2E::zero()]
+    }
+
+    fn one() -> Self::BaseType {
+        [Fp2E::one(), Fp2E::zero()]
+    }
+
+    fn from_u64(x: u64) -> Self::BaseType {
+        [Fp2E::from(x), Fp2E::zero()]
+    }
+
+    fn from_base_type(x: Self::BaseType) -> Self::BaseType {
+        x
+    }
+}
+
+impl IsSubFieldOf<Degree4ExtensionField> for Mersenne31Field {
+    fn add(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        [FpE::from(a) + &b[0], b[1].clone()]
+    }
+
+    fn sub(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        [FpE::from(a) - &b[0], -&b[1]]
+    }
+
+    fn mul(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        let c0 = FpE::from(a) * &b[0];
+        let c1 = FpE::from(a) * &b[1];
+        [c0, c1]
+    }
+
+    fn div(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        let b_inv = Degree4ExtensionField::inv(b).unwrap();
+        <Self as IsSubFieldOf<Degree4ExtensionField>>::mul(a, &b_inv)
+    }
+
+    fn embed(a: Self::BaseType) -> <Degree4ExtensionField as IsField>::BaseType {
+        [
+            Fp2E::from_raw(<Self as IsSubFieldOf<Degree2ExtensionField>>::embed(a)),
+            Fp2E::zero(),
+        ]
+    }
+
+    #[cfg(feature = "alloc")]
+    fn to_subfield_vec(
+        b: <Degree4ExtensionField as IsField>::BaseType,
+    ) -> alloc::vec::Vec<Self::BaseType> {
+        // TODO: Repace this for with a map similarly to this:
+        // b.into_iter().map(|x| x.to_raw()).collect()
+        let mut result = Vec::new();
+        for fp2e in b {
+            result.push(fp2e.value()[0].to_raw());
+            result.push(fp2e.value()[1].to_raw());
+        }
+        result
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use core::ops::Neg;
+
+    use crate::field::fields::mersenne31::field::MERSENNE_31_PRIME_FIELD_ORDER;
+
+    use super::*;
+
+    type FpE = FieldElement<Mersenne31Field>;
+    type Fp2E = FieldElement<Degree2ExtensionField>;
+    type Fp4E = FieldElement<Degree4ExtensionField>;
+
+    #[test]
+    fn add_real_one_plus_one_is_two() {
+        assert_eq!(Fp2E::one() + Fp2E::one(), Fp2E::from(2))
+    }
+
+    #[test]
+    fn add_real_neg_one_plus_one_is_zero() {
+        assert_eq!(Fp2E::one() + Fp2E::one().neg(), Fp2E::zero())
+    }
+
+    #[test]
+    fn add_real_neg_one_plus_two_is_one() {
+        assert_eq!(Fp2E::one().neg() + Fp2E::from(2), Fp2E::one())
+    }
+
+    #[test]
+    fn add_real_neg_one_plus_neg_one_is_order_sub_two() {
+        assert_eq!(
+            Fp2E::one().neg() + Fp2E::one().neg(),
+            Fp2E::new([FpE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 2)), FpE::zero()])
+        )
+    }
+
+    #[test]
+    fn add_complex_one_plus_one_two() {
+        let one_i = Fp2E::new([FpE::zero(), FpE::one()]);
+        let two_i = Fp2E::new([FpE::zero(), FpE::from(2)]);
+        assert_eq!(&one_i + &one_i, two_i)
+    }
+
+    #[test]
+    fn add_complex_neg_one_plus_one_is_zero() {
+        //Manually declare the complex part to one
+        let neg_one_i = Fp2E::new([FpE::zero(), -FpE::one()]);
+        let one_i = Fp2E::new([FpE::zero(), FpE::one()]);
+        assert_eq!(neg_one_i + one_i, Fp2E::zero())
+    }
+
+    #[test]
+    fn add_complex_neg_one_plus_two_is_one() {
+        let neg_one_i = Fp2E::new([FpE::zero(), -FpE::one()]);
+        let two_i = Fp2E::new([FpE::zero(), FpE::from(2)]);
+        let one_i = Fp2E::new([FpE::zero(), FpE::one()]);
+        assert_eq!(&neg_one_i + &two_i, one_i)
+    }
+
+    #[test]
+    fn add_complex_neg_one_plus_neg_one_imag_is_order_sub_two() {
+        let neg_one_i = Fp2E::new([FpE::zero(), -FpE::one()]);
+        assert_eq!(
+            (&neg_one_i + &neg_one_i).value()[1],
+            FpE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 2))
+        )
+    }
+
+    #[test]
+    fn add_order() {
+        let a = Fp2E::new([-FpE::one(), FpE::one()]);
+        let b = Fp2E::new([
+            FpE::from(2),
+            FpE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 2)),
+        ]);
+        let c = Fp2E::new([FpE::one(), -FpE::one()]);
+        assert_eq!(&a + &b, c)
+    }
+
+    #[test]
+    fn add_equal_zero() {
+        let a = Fp2E::new([-FpE::one(), -FpE::one()]);
+        let b = Fp2E::new([FpE::one(), FpE::one()]);
+        assert_eq!(&a + &b, Fp2E::zero())
+    }
+
+    #[test]
+    fn add_plus_one() {
+        let a = Fp2E::new([FpE::one(), FpE::from(2)]);
+        let b = Fp2E::new([FpE::one(), FpE::one()]);
+        let c = Fp2E::new([FpE::from(2), FpE::from(3)]);
+        assert_eq!(&a + &b, c)
+    }
+
+    #[test]
+    fn sub_real_one_sub_one_is_zero() {
+        assert_eq!(&Fp2E::one() - &Fp2E::one(), Fp2E::zero())
+    }
+
+    #[test]
+    fn sub_real_two_sub_two_is_zero() {
+        assert_eq!(&Fp2E::from(2) - &Fp2E::from(2), Fp2E::zero())
+    }
+
+    #[test]
+    fn sub_real_neg_one_sub_neg_one_is_zero() {
+        assert_eq!(Fp2E::one().neg() - Fp2E::one().neg(), Fp2E::zero())
+    }
+
+    #[test]
+    fn sub_real_two_sub_one_is_one() {
+        assert_eq!(Fp2E::from(2) - Fp2E::one(), Fp2E::one())
+    }
+
+    #[test]
+    fn sub_real_neg_one_sub_zero_is_neg_one() {
+        assert_eq!(Fp2E::one().neg() - Fp2E::zero(), Fp2E::one().neg())
+    }
+
+    #[test]
+    fn sub_complex_one_sub_one_is_zero() {
+        let one = Fp2E::new([FpE::zero(), FpE::one()]);
+        assert_eq!(&one - &one, Fp2E::zero())
+    }
+
+    #[test]
+    fn sub_complex_two_sub_two_is_zero() {
+        let two = Fp2E::new([FpE::zero(), FpE::from(2)]);
+        assert_eq!(&two - &two, Fp2E::zero())
+    }
+
+    #[test]
+    fn sub_complex_neg_one_sub_neg_one_is_zero() {
+        let neg_one = Fp2E::new([FpE::zero(), -FpE::one()]);
+        assert_eq!(&neg_one - &neg_one, Fp2E::zero())
+    }
+
+    #[test]
+    fn sub_complex_two_sub_one_is_one() {
+        let two = Fp2E::new([FpE::zero(), FpE::from(2)]);
+        let one = Fp2E::new([FpE::zero(), FpE::one()]);
+        assert_eq!(&two - &one, one)
+    }
+
+    #[test]
+    fn sub_complex_neg_one_sub_zero_is_neg_one() {
+        let neg_one = Fp2E::new([FpE::zero(), -FpE::one()]);
+        assert_eq!(&neg_one - &Fp2E::zero(), neg_one)
+    }
+
+    #[test]
+    fn mul_fp2_is_correct() {
+        let a = Fp2E::new([FpE::from(2), FpE::from(2)]);
+        let b = Fp2E::new([FpE::from(4), FpE::from(5)]);
+        let c = Fp2E::new([-FpE::from(2), FpE::from(18)]);
+        assert_eq!(&a * &b, c)
+    }
+
+    #[test]
+    fn square_equals_mul_by_itself() {
+        let a = Fp2E::new([FpE::from(2), FpE::from(3)]);
+        assert_eq!(a.square(), &a * &a)
+    }
+
+    #[test]
+    fn test_fp2_add() {
+        let a = Fp2E::new([FpE::from(0), FpE::from(3)]);
+        let b = Fp2E::new([-FpE::from(2), FpE::from(8)]);
+        let expected_result = Fp2E::new([FpE::from(0) - FpE::from(2), FpE::from(3) + FpE::from(8)]);
+        assert_eq!(a + b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_add_2() {
+        let a = Fp2E::new([FpE::from(2), FpE::from(4)]);
+        let b = Fp2E::new([-FpE::from(2), -FpE::from(4)]);
+        let expected_result = Fp2E::new([FpE::from(2) - FpE::from(2), FpE::from(4) - FpE::from(4)]);
+        assert_eq!(a + b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_add_3() {
+        let a = Fp2E::new([FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER), FpE::from(1)]);
+        let b = Fp2E::new([FpE::from(1), FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER)]);
+        let expected_result = Fp2E::new([FpE::from(1), FpE::from(1)]);
+        assert_eq!(a + b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_sub() {
+        let a = Fp2E::new([FpE::from(0), FpE::from(3)]);
+        let b = Fp2E::new([-FpE::from(2), FpE::from(8)]);
+        let expected_result = Fp2E::new([FpE::from(0) + FpE::from(2), FpE::from(3) - FpE::from(8)]);
+        assert_eq!(a - b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_sub_2() {
+        let a = Fp2E::new([FpE::zero(), FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER)]);
+        let b = Fp2E::new([FpE::one(), -FpE::one()]);
+        let expected_result =
+            Fp2E::new([FpE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 1)), FpE::one()]);
+        assert_eq!(a - b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_sub_3() {
+        let a = Fp2E::new([FpE::from(5), FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER)]);
+        let b = Fp2E::new([FpE::from(5), FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER)]);
+        let expected_result = Fp2E::new([FpE::zero(), FpE::zero()]);
+        assert_eq!(a - b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_mul() {
+        let a = Fp2E::new([FpE::from(12), FpE::from(5)]);
+        let b = Fp2E::new([-FpE::from(4), FpE::from(2)]);
+        let expected_result = Fp2E::new([-FpE::from(58), FpE::new(4)]);
+        assert_eq!(a * b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_mul_2() {
+        let a = Fp2E::new([FpE::one(), FpE::zero()]);
+        let b = Fp2E::new([FpE::from(12), -FpE::from(8)]);
+        let expected_result = Fp2E::new([FpE::from(12), -FpE::new(8)]);
+        assert_eq!(a * b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_mul_3() {
+        let a = Fp2E::new([FpE::zero(), FpE::zero()]);
+        let b = Fp2E::new([FpE::from(2), FpE::from(7)]);
+        let expected_result = Fp2E::new([FpE::zero(), FpE::zero()]);
+        assert_eq!(a * b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_mul_4() {
+        let a = Fp2E::new([FpE::from(2), FpE::from(7)]);
+        let b = Fp2E::new([FpE::zero(), FpE::zero()]);
+        let expected_result = Fp2E::new([FpE::zero(), FpE::zero()]);
+        assert_eq!(a * b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_mul_5() {
+        let a = Fp2E::new([FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER), FpE::one()]);
+        let b = Fp2E::new([FpE::from(2), FpE::from(&MERSENNE_31_PRIME_FIELD_ORDER)]);
+        let expected_result = Fp2E::new([FpE::zero(), FpE::from(2)]);
+        assert_eq!(a * b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_inv() {
+        let a = Fp2E::new([FpE::one(), FpE::zero()]);
+        let expected_result = Fp2E::new([FpE::one(), FpE::zero()]);
+        assert_eq!(a.inv().unwrap(), expected_result);
+    }
+
+    #[test]
+    fn test_fp2_inv_2() {
+        let a = Fp2E::new([FpE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 1)), FpE::one()]);
+        let expected_result = Fp2E::new([FpE::from(1073741823), FpE::from(1073741823)]);
+        assert_eq!(a.inv().unwrap(), expected_result);
+    }
+
+    #[test]
+    fn test_fp2_inv_3() {
+        let a = Fp2E::new([FpE::from(2063384121), FpE::from(1232183486)]);
+        let expected_result = Fp2E::new([FpE::from(1244288232), FpE::from(1321511038)]);
+        assert_eq!(a.inv().unwrap(), expected_result);
+    }
+
+    #[test]
+    fn test_fp2_mul_inv() {
+        let a = Fp2E::new([FpE::from(12), FpE::from(5)]);
+        let b = a.inv().unwrap();
+        let expected_result = Fp2E::new([FpE::one(), FpE::zero()]);
+        assert_eq!(a * b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_div() {
+        let a = Fp2E::new([FpE::from(12), FpE::from(5)]);
+        let b = Fp2E::new([FpE::from(4), FpE::from(2)]);
+        let expected_result = Fp2E::new([FpE::from(644245097), FpE::from(1288490188)]);
+        assert_eq!(a / b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_div_2() {
+        let a = Fp2E::new([FpE::from(4), FpE::from(7)]);
+        let b = Fp2E::new([FpE::one(), FpE::zero()]);
+        let expected_result = Fp2E::new([FpE::from(4), FpE::from(7)]);
+        assert_eq!(a / b, expected_result);
+    }
+
+    #[test]
+    fn test_fp2_div_3() {
+        let a = Fp2E::new([FpE::zero(), FpE::zero()]);
+        let b = Fp2E::new([FpE::from(3), FpE::from(12)]);
+        let expected_result = Fp2E::new([FpE::zero(), FpE::zero()]);
+        assert_eq!(a / b, expected_result);
+    }
+
+    #[test]
+    fn mul_fp4_by_zero_is_zero() {
+        let a = Fp4E::new([
+            Fp2E::new([FpE::from(2), FpE::from(3)]),
+            Fp2E::new([FpE::from(4), FpE::from(5)]),
+        ]);
+        assert_eq!(Fp4E::zero(), a * Fp4E::zero())
+    }
+
+    #[test]
+    fn mul_fp4_by_one_is_identity() {
+        let a = Fp4E::new([
+            Fp2E::new([FpE::from(2), FpE::from(3)]),
+            Fp2E::new([FpE::from(4), FpE::from(5)]),
+        ]);
+        assert_eq!(a, a.clone() * Fp4E::one())
+    }
+
+    #[test]
+    fn square_fp4_equals_mul_two_times() {
+        let a = Fp4E::new([
+            Fp2E::new([FpE::from(3), FpE::from(4)]),
+            Fp2E::new([FpE::from(5), FpE::from(6)]),
+        ]);
+
+        assert_eq!(a.square(), &a * &a)
+    }
+
+    #[test]
+    fn fp4_mul_by_inv_is_one() {
+        let a = Fp4E::new([
+            Fp2E::new([FpE::from(2147483647), FpE::from(2147483648)]),
+            Fp2E::new([FpE::from(2147483649), FpE::from(2147483650)]),
+        ]);
+
+        assert_eq!(&a * a.inv().unwrap(), Fp4E::one())
+    }
+
+    #[test]
+    fn embed_fp_with_fp4() {
+        let a = FpE::from(3);
+        let a_extension = Fp4E::from(3);
+        assert_eq!(a.to_extension::<Degree4ExtensionField>(), a_extension);
+    }
+
+    #[test]
+    fn add_fp_and_fp4() {
+        let a = FpE::from(3);
+        let a_extension = Fp4E::from(3);
+        let b = Fp4E::from(2);
+        assert_eq!(a + &b, a_extension + b);
+    }
+
+    #[test]
+    fn mul_fp_by_fp4() {
+        let a = FpE::from(30000000000);
+        let a_extension = a.to_extension::<Degree4ExtensionField>();
+        let b = Fp4E::new([
+            Fp2E::new([FpE::from(1), FpE::from(2)]),
+            Fp2E::new([FpE::from(3), FpE::from(4)]),
+        ]);
+        assert_eq!(a * &b, a_extension * b);
+    }
+}
diff --git a/math/src/field/fields/mersenne31/field.rs b/math/src/field/fields/mersenne31/field.rs
index e4abfab0f..1c8b2dc58 100644
--- a/math/src/field/fields/mersenne31/field.rs
+++ b/math/src/field/fields/mersenne31/field.rs
@@ -42,6 +42,29 @@ impl Mersenne31Field {
         // Delayed reduction
         Self::from_u64(iter.map(|x| (x as u64)).sum::<u64>())
     }
+
+    /// Computes a * 2^k, with 0 < k < 31
+    pub fn mul_power_two(a: u32, k: u32) -> u32 {
+        let msb = (a & (u32::MAX << (31 - k))) >> (31 - k); // The k + 1 msf shifted right .
+        let lsb = (a & (u32::MAX >> (k + 1))) << k; // The 31 - k lsb shifted left.
+        Self::weak_reduce(msb + lsb)
+    }
+
+    pub fn pow_2(a: &u32, order: u32) -> u32 {
+        let mut res = *a;
+        (0..order).for_each(|_| res = Self::square(&res));
+        res
+    }
+
+    /// TODO: See if we can optimize this function.
+    /// Computes 2a^2 - 1
+    pub fn two_square_minus_one(a: &u32) -> u32 {
+        if *a == 0 {
+            MERSENNE_31_PRIME_FIELD_ORDER - 1
+        } else {
+            Self::from_u64(((u64::from(*a) * u64::from(*a)) << 1) - 1)
+        }
+    }
 }
 
 pub const MERSENNE_31_PRIME_FIELD_ORDER: u32 = (1 << 31) - 1;
@@ -54,18 +77,9 @@ impl IsField for Mersenne31Field {
 
     /// Returns the sum of `a` and `b`.
     fn add(a: &u32, b: &u32) -> u32 {
-        // Avoids conditional https://github.com/Plonky3/Plonky3/blob/6049a30c3b1f5351c3eb0f7c994dc97e8f68d10d/mersenne-31/src/lib.rs#L249
-        // Working with i32 means we get a flag which informs us if overflow happens
-        let (sum_i32, over) = (*a as i32).overflowing_add(*b as i32);
-        let sum_u32 = sum_i32 as u32;
-        let sum_corr = sum_u32.wrapping_sub(MERSENNE_31_PRIME_FIELD_ORDER);
-
-        //assert 31 bit clear
-        // If self + rhs did not overflow, return it.
-        // If self + rhs overflowed, sum_corr = self + rhs - (2**31 - 1).
-        let sum = if over { sum_corr } else { sum_u32 };
-        debug_assert!((sum >> 31) == 0);
-        Self::as_representative(&sum)
+        // We are using that if a and b are field elements of Mersenne31, then
+        // a + b has at most 32 bits, so we can use the weak_reduce function to take mudulus p.
+        Self::weak_reduce(a + b)
     }
 
     /// Returns the multiplication of `a` and `b`.
@@ -75,13 +89,7 @@ impl IsField for Mersenne31Field {
     }
 
     fn sub(a: &u32, b: &u32) -> u32 {
-        let (mut sub, over) = a.overflowing_sub(*b);
-
-        // If we didn't overflow we have the correct value.
-        // Otherwise we have added 2**32 = 2**31 + 1 mod 2**31 - 1.
-        // Hence we need to remove the most significant bit and subtract 1.
-        sub -= over as u32;
-        sub & MERSENNE_31_PRIME_FIELD_ORDER
+        Self::weak_reduce(a + MERSENNE_31_PRIME_FIELD_ORDER - b)
     }
 
     /// Returns the additive inverse of `a`.
@@ -91,20 +99,20 @@ impl IsField for Mersenne31Field {
     }
 
     /// Returns the multiplicative inverse of `a`.
-    fn inv(a: &u32) -> Result<u32, FieldError> {
-        if *a == Self::zero() || *a == MERSENNE_31_PRIME_FIELD_ORDER {
+    fn inv(x: &u32) -> Result<u32, FieldError> {
+        if *x == Self::zero() || *x == MERSENNE_31_PRIME_FIELD_ORDER {
             return Err(FieldError::InvZeroError);
         }
-        let p101 = Self::mul(&Self::pow(a, 4u32), a);
+        let p101 = Self::mul(&Self::pow_2(x, 2), x);
         let p1111 = Self::mul(&Self::square(&p101), &p101);
-        let p11111111 = Self::mul(&Self::pow(&p1111, 16u32), &p1111);
-        let p111111110000 = Self::pow(&p11111111, 16u32);
+        let p11111111 = Self::mul(&Self::pow_2(&p1111, 4u32), &p1111);
+        let p111111110000 = Self::pow_2(&p11111111, 4u32);
         let p111111111111 = Self::mul(&p111111110000, &p1111);
-        let p1111111111111111 = Self::mul(&Self::pow(&p111111110000, 16u32), &p11111111);
+        let p1111111111111111 = Self::mul(&Self::pow_2(&p111111110000, 4u32), &p11111111);
         let p1111111111111111111111111111 =
-            Self::mul(&Self::pow(&p1111111111111111, 4096u32), &p111111111111);
+            Self::mul(&Self::pow_2(&p1111111111111111, 12u32), &p111111111111);
         let p1111111111111111111111111111101 =
-            Self::mul(&Self::pow(&p1111111111111111111111111111, 8u32), &p101);
+            Self::mul(&Self::pow_2(&p1111111111111111111111111111, 3u32), &p101);
         Ok(p1111111111111111111111111111101)
     }
 
@@ -120,7 +128,7 @@ impl IsField for Mersenne31Field {
     }
 
     /// Returns the additive neutral element.
-    fn zero() -> Self::BaseType {
+    fn zero() -> u32 {
         0u32
     }
 
@@ -131,16 +139,7 @@ impl IsField for Mersenne31Field {
 
     /// Returns the element `x * 1` where 1 is the multiplicative neutral element.
     fn from_u64(x: u64) -> u32 {
-        let (lo, hi) = (x as u32 as u64, x >> 32);
-        // 2^32 = 2 (mod Mersenne 31 bit prime)
-        // t <= (2^32 - 1) + 2 * (2^32 - 1) = 3 * 2^32 - 3 = 6 * 2^31 - 3
-        let t = lo + 2 * hi;
-
-        const MASK: u64 = (1 << 31) - 1;
-        let (lo, hi) = ((t & MASK) as u32, (t >> 31) as u32);
-        // 2^31 = 1 mod Mersenne31
-        // lo < 2^31, hi < 6, so lo + hi < 2^32.
-        Self::weak_reduce(lo + hi)
+        (((((x >> 31) + x + 1) >> 31) + x) & (MERSENNE_31_PRIME_FIELD_ORDER as u64)) as u32
     }
 
     /// Takes as input an element of BaseType and returns the internal representation
@@ -148,6 +147,9 @@ impl IsField for Mersenne31Field {
     fn from_base_type(x: u32) -> u32 {
         Self::weak_reduce(x)
     }
+    fn double(a: &u32) -> u32 {
+        Self::weak_reduce(a << 1)
+    }
 }
 
 impl IsPrimeField for Mersenne31Field {
@@ -205,12 +207,45 @@ impl Display for FieldElement<Mersenne31Field> {
 mod tests {
     use super::*;
     type F = Mersenne31Field;
+    type FE = FieldElement<F>;
+
+    #[test]
+    fn mul_power_two_is_correct() {
+        let a = 3u32;
+        let k = 2;
+        let expected_result = FE::from(&a) * FE::from(2).pow(k as u16);
+        let result = F::mul_power_two(a, k);
+        assert_eq!(FE::from(&result), expected_result)
+    }
+
+    #[test]
+    fn mul_power_two_is_correct_2() {
+        let a = 229287u32;
+        let k = 4;
+        let expected_result = FE::from(&a) * FE::from(2).pow(k as u16);
+        let result = F::mul_power_two(a, k);
+        assert_eq!(FE::from(&result), expected_result)
+    }
+
+    #[test]
+    fn pow_2_is_correct() {
+        let a = 3u32;
+        let order = 12;
+        let result = F::pow_2(&a, order);
+        let expected_result = FE::pow(&FE::from(&a), 4096u32);
+        assert_eq!(FE::from(&result), expected_result)
+    }
 
     #[test]
     fn from_hex_for_b_is_11() {
         assert_eq!(F::from_hex("B").unwrap(), 11);
     }
 
+    #[test]
+    fn from_hex_for_b_is_11_v2() {
+        assert_eq!(FE::from_hex("B").unwrap(), FE::from(11));
+    }
+
     #[test]
     fn sum_delayed_reduction() {
         let up_to = u32::pow(2, 16);
@@ -236,190 +271,195 @@ mod tests {
 
     #[test]
     fn one_plus_1_is_2() {
-        let a = F::one();
-        let b = F::one();
-        let c = F::add(&a, &b);
-        assert_eq!(c, 2u32);
+        assert_eq!(FE::one() + FE::one(), FE::from(&2u32));
     }
 
     #[test]
     fn neg_1_plus_1_is_0() {
-        let a = F::neg(&F::one());
-        let b = F::one();
-        let c = F::add(&a, &b);
-        assert_eq!(c, F::zero());
+        assert_eq!(-FE::one() + FE::one(), FE::zero());
     }
 
     #[test]
     fn neg_1_plus_2_is_1() {
-        let a = F::neg(&F::one());
-        let b = F::from_base_type(2u32);
-        let c = F::add(&a, &b);
-        assert_eq!(c, F::one());
+        assert_eq!(-FE::one() + FE::from(&2u32), FE::one());
     }
 
     #[test]
     fn max_order_plus_1_is_0() {
-        let a = F::from_base_type(MERSENNE_31_PRIME_FIELD_ORDER - 1);
-        let b = F::one();
-        let c = F::add(&a, &b);
-        assert_eq!(c, F::zero());
+        assert_eq!(
+            FE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 1)) + FE::from(1),
+            FE::from(0)
+        );
     }
 
     #[test]
     fn comparing_13_and_13_are_equal() {
-        let a = F::from_base_type(13);
-        let b = F::from_base_type(13);
-        assert_eq!(a, b);
+        assert_eq!(FE::from(&13u32), FE::from(13));
     }
 
     #[test]
     fn comparing_13_and_8_they_are_not_equal() {
-        let a = F::from_base_type(13);
-        let b = F::from_base_type(8);
-        assert_ne!(a, b);
+        assert_ne!(FE::from(&13u32), FE::from(8));
     }
 
     #[test]
     fn one_sub_1_is_0() {
-        let a = F::one();
-        let b = F::one();
-        let c = F::sub(&a, &b);
-        assert_eq!(c, F::zero());
+        assert_eq!(FE::one() - FE::one(), FE::zero());
     }
 
     #[test]
     fn zero_sub_1_is_order_minus_1() {
-        let a = F::zero();
-        let b = F::one();
-        let c = F::sub(&a, &b);
-        assert_eq!(c, MERSENNE_31_PRIME_FIELD_ORDER - 1);
+        assert_eq!(
+            FE::zero() - FE::one(),
+            FE::from(&(MERSENNE_31_PRIME_FIELD_ORDER - 1))
+        );
     }
 
     #[test]
     fn neg_1_sub_neg_1_is_0() {
-        let a = F::neg(&F::one());
-        let b = F::neg(&F::one());
-        let c = F::sub(&a, &b);
-        assert_eq!(c, F::zero());
+        assert_eq!(-FE::one() - (-FE::one()), FE::zero());
     }
 
     #[test]
-    fn neg_1_sub_1_is_neg_1() {
-        let a = F::neg(&F::one());
-        let b = F::zero();
-        let c = F::sub(&a, &b);
-        assert_eq!(c, F::neg(&F::one()));
+    fn neg_1_sub_0_is_neg_1() {
+        assert_eq!(-FE::one() - FE::zero(), -FE::one());
     }
 
     #[test]
     fn mul_neutral_element() {
-        let a = F::from_base_type(1);
-        let b = F::from_base_type(2);
-        let c = F::mul(&a, &b);
-        assert_eq!(c, F::from_base_type(2));
+        assert_eq!(FE::one() * FE::from(&2u32), FE::from(&2u32));
     }
 
     #[test]
     fn mul_2_3_is_6() {
-        let a = F::from_base_type(2);
-        let b = F::from_base_type(3);
-        assert_eq!(a * b, F::from_base_type(6));
+        assert_eq!(FE::from(&2u32) * FE::from(&3u32), FE::from(&6u32));
     }
 
     #[test]
     fn mul_order_neg_1() {
-        let a = F::from_base_type(MERSENNE_31_PRIME_FIELD_ORDER - 1);
-        let b = F::from_base_type(MERSENNE_31_PRIME_FIELD_ORDER - 1);
-        let c = F::mul(&a, &b);
-        assert_eq!(c, F::from_base_type(1));
+        assert_eq!(
+            FE::from(MERSENNE_31_PRIME_FIELD_ORDER as u64 - 1)
+                * FE::from(MERSENNE_31_PRIME_FIELD_ORDER as u64 - 1),
+            FE::one()
+        );
     }
 
     #[test]
     fn pow_p_neg_1() {
         assert_eq!(
-            F::pow(&F::from_base_type(2), MERSENNE_31_PRIME_FIELD_ORDER - 1),
-            F::one()
+            FE::pow(&FE::from(&2u32), MERSENNE_31_PRIME_FIELD_ORDER - 1),
+            FE::one()
         )
     }
 
     #[test]
     fn inv_0_error() {
-        let result = F::inv(&F::zero());
+        let result = FE::inv(&FE::zero());
         assert!(matches!(result, Err(FieldError::InvZeroError)));
     }
 
     #[test]
     fn inv_2() {
-        let result = F::inv(&F::from_base_type(2u32)).unwrap();
+        let result = FE::inv(&FE::from(&2u32)).unwrap();
         // sage: 1 / F(2) = 1073741824
-        assert_eq!(result, 1073741824);
+        assert_eq!(result, FE::from(1073741824));
     }
 
     #[test]
     fn pow_2_3() {
-        assert_eq!(F::pow(&F::from_base_type(2), 3_u64), 8)
+        assert_eq!(FE::pow(&FE::from(&2u32), 3u64), FE::from(8));
     }
 
     #[test]
     fn div_1() {
-        assert_eq!(F::div(&F::from_base_type(2), &F::from_base_type(1)), 2)
+        assert_eq!(FE::from(&2u32) / FE::from(&1u32), FE::from(&2u32));
     }
 
     #[test]
     fn div_4_2() {
-        assert_eq!(F::div(&F::from_base_type(4), &F::from_base_type(2)), 2)
+        assert_eq!(FE::from(&4u32) / FE::from(&2u32), FE::from(&2u32));
     }
 
-    // 1431655766
     #[test]
     fn div_4_3() {
         // sage: F(4) / F(3) = 1431655766
-        assert_eq!(
-            F::div(&F::from_base_type(4), &F::from_base_type(3)),
-            1431655766
-        )
+        assert_eq!(FE::from(&4u32) / FE::from(&3u32), FE::from(1431655766));
     }
 
     #[test]
     fn two_plus_its_additive_inv_is_0() {
-        let two = F::from_base_type(2);
-
-        assert_eq!(F::add(&two, &F::neg(&two)), F::zero())
+        assert_eq!(FE::from(&2u32) + (-FE::from(&2u32)), FE::zero());
     }
 
     #[test]
     fn from_u64_test() {
-        let num = F::from_u64(1u64);
-        assert_eq!(num, F::one());
+        assert_eq!(FE::from(1u64), FE::one());
     }
 
     #[test]
     fn creating_a_field_element_from_its_representative_returns_the_same_element_1() {
-        let change = 1;
-        let f1 = F::from_base_type(MERSENNE_31_PRIME_FIELD_ORDER + change);
-        let f2 = F::from_base_type(Mersenne31Field::representative(&f1));
+        let change: u32 = MERSENNE_31_PRIME_FIELD_ORDER + 1;
+        let f1 = FE::from(&change);
+        let f2 = FE::from(&FE::representative(&f1));
         assert_eq!(f1, f2);
     }
 
     #[test]
     fn creating_a_field_element_from_its_representative_returns_the_same_element_2() {
-        let change = 8;
-        let f1 = F::from_base_type(MERSENNE_31_PRIME_FIELD_ORDER + change);
-        let f2 = F::from_base_type(Mersenne31Field::representative(&f1));
+        let change: u32 = MERSENNE_31_PRIME_FIELD_ORDER + 8;
+        let f1 = FE::from(&change);
+        let f2 = FE::from(&FE::representative(&f1));
         assert_eq!(f1, f2);
     }
 
     #[test]
     fn from_base_type_test() {
-        let b = F::from_base_type(1u32);
-        assert_eq!(b, F::one());
+        assert_eq!(FE::from(&1u32), FE::one());
     }
 
     #[cfg(feature = "std")]
     #[test]
     fn to_hex_test() {
-        let num = F::from_hex("B").unwrap();
-        assert_eq!(F::to_hex(&num), "B");
+        let num = FE::from_hex("B").unwrap();
+        assert_eq!(FE::to_hex(&num), "B");
+    }
+
+    #[test]
+    fn double_equals_add_itself() {
+        let a = FE::from(1234);
+        assert_eq!(a + a, a.double())
+    }
+
+    #[test]
+    fn two_square_minus_one_is_correct() {
+        let a = FE::from(2147483650);
+        assert_eq!(
+            FE::from(&F::two_square_minus_one(a.value())),
+            a.square().double() - FE::one()
+        )
+    }
+
+    #[test]
+    fn two_square_zero_minus_one_is_minus_one() {
+        let a = FE::from(0);
+        assert_eq!(
+            FE::from(&F::two_square_minus_one(a.value())),
+            a.square().double() - FE::one()
+        )
+    }
+
+    #[test]
+    fn two_square_p_minus_one_is_minus_one() {
+        let a = FE::from(&MERSENNE_31_PRIME_FIELD_ORDER);
+        assert_eq!(
+            FE::from(&F::two_square_minus_one(a.value())),
+            a.square().double() - FE::one()
+        )
+    }
+
+    #[test]
+    fn mul_by_inv() {
+        let x = 3476715743_u32;
+        assert_eq!(FE::from(&x).inv().unwrap() * FE::from(&x), FE::one());
     }
 }
diff --git a/math/src/field/fields/mersenne31/mod.rs b/math/src/field/fields/mersenne31/mod.rs
index 4bfd3daf1..2272e7d5e 100644
--- a/math/src/field/fields/mersenne31/mod.rs
+++ b/math/src/field/fields/mersenne31/mod.rs
@@ -1,2 +1,2 @@
-pub mod extension;
+pub mod extensions;
 pub mod field;
diff --git a/math/src/field/fields/p448_goldilocks_prime_field.rs b/math/src/field/fields/p448_goldilocks_prime_field.rs
index 9f67fd1d2..eadd61541 100644
--- a/math/src/field/fields/p448_goldilocks_prime_field.rs
+++ b/math/src/field/fields/p448_goldilocks_prime_field.rs
@@ -16,7 +16,7 @@ pub const P448_GOLDILOCKS_PRIME_FIELD_ORDER: U448 =
 /// 448-bit unsigned integer represented as
 /// a size 8 `u64` array `limbs` of 56-bit words.
 /// The least significant word is in the left most position.
-#[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord)]
+#[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Default)]
 pub struct U56x8 {
     limbs: [u64; 8],
 }
diff --git a/math/src/field/traits.rs b/math/src/field/traits.rs
index a9a39d2dd..320531f42 100644
--- a/math/src/field/traits.rs
+++ b/math/src/field/traits.rs
@@ -99,9 +99,9 @@ pub trait IsField: Debug + Clone {
     /// The underlying base type for representing elements from the field.
     // TODO: Relax Unpin for non cuda usage
     #[cfg(feature = "lambdaworks-serde-binary")]
-    type BaseType: Clone + Debug + Unpin + ByteConversion;
+    type BaseType: Clone + Debug + Unpin + ByteConversion + Default;
     #[cfg(not(feature = "lambdaworks-serde-binary"))]
-    type BaseType: Clone + Debug + Unpin;
+    type BaseType: Clone + Debug + Unpin + Default;
 
     /// Returns the sum of `a` and `b`.
     fn add(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType;
@@ -173,7 +173,9 @@ pub trait IsField: Debug + Clone {
     fn eq(a: &Self::BaseType, b: &Self::BaseType) -> bool;
 
     /// Returns the additive neutral element.
-    fn zero() -> Self::BaseType;
+    fn zero() -> Self::BaseType {
+        Self::BaseType::default()
+    }
 
     /// Returns the multiplicative neutral element.
     fn one() -> Self::BaseType;
diff --git a/math/src/unsigned_integer/element.rs b/math/src/unsigned_integer/element.rs
index 021afd9a3..1613e192b 100644
--- a/math/src/unsigned_integer/element.rs
+++ b/math/src/unsigned_integer/element.rs
@@ -36,6 +36,14 @@ pub struct UnsignedInteger<const NUM_LIMBS: usize> {
     pub limbs: [u64; NUM_LIMBS],
 }
 
+impl<const NUM_LIMBS: usize> Default for UnsignedInteger<NUM_LIMBS> {
+    fn default() -> Self {
+        Self {
+            limbs: [0; NUM_LIMBS],
+        }
+    }
+}
+
 // NOTE: manually implementing `PartialOrd` may seem unorthodox, but the
 // derived implementation had terrible performance.
 #[allow(clippy::non_canonical_partial_ord_impl)]