From 9617e52e88d5422787192244cd9648727b9cd4ab Mon Sep 17 00:00:00 2001
From: erhant <erhany96@gmail.com>
Date: Tue, 1 Oct 2024 17:16:14 +0300
Subject: [PATCH 1/3] add prime field for stark101 tutorial (#875)

Co-authored-by: Diego K <43053772+diegokingston@users.noreply.github.com>
---
 math/src/field/fields/fft_friendly/mod.rs     |   2 +
 .../fft_friendly/stark_101_prime_field.rs     | 113 ++++++++++++++++++
 2 files changed, 115 insertions(+)
 create mode 100644 math/src/field/fields/fft_friendly/stark_101_prime_field.rs

diff --git a/math/src/field/fields/fft_friendly/mod.rs b/math/src/field/fields/fft_friendly/mod.rs
index a9b2f0f2e..535b94ecc 100644
--- a/math/src/field/fields/fft_friendly/mod.rs
+++ b/math/src/field/fields/fft_friendly/mod.rs
@@ -2,6 +2,8 @@
 pub mod babybear;
 /// Implemenation of the quadratic extension of the babybear field
 pub mod quadratic_babybear;
+/// Implementation of the prime field used in [Stark101](https://starkware.co/stark-101/) tutorial, p = 3 * 2^30 + 1
+pub mod stark_101_prime_field;
 /// Implementation of two-adic prime field over 256 bit unsigned integers.
 pub mod stark_252_prime_field;
 /// Implemenation of the Goldilocks Prime Field p = 2^64 - 2^32 + 1
diff --git a/math/src/field/fields/fft_friendly/stark_101_prime_field.rs b/math/src/field/fields/fft_friendly/stark_101_prime_field.rs
new file mode 100644
index 000000000..d6232a414
--- /dev/null
+++ b/math/src/field/fields/fft_friendly/stark_101_prime_field.rs
@@ -0,0 +1,113 @@
+use crate::{
+    field::{
+        element::FieldElement,
+        fields::montgomery_backed_prime_fields::{IsModulus, U64PrimeField},
+        traits::IsFFTField,
+    },
+    unsigned_integer::element::{UnsignedInteger, U64},
+};
+
+#[derive(Clone, Debug, Hash, Copy)]
+pub struct MontgomeryConfigStark101PrimeField;
+impl IsModulus<U64> for MontgomeryConfigStark101PrimeField {
+    /// 3 * 2^30 + 1
+    const MODULUS: U64 = U64::from_hex_unchecked("c0000001");
+}
+
+pub type Stark101PrimeField = U64PrimeField<MontgomeryConfigStark101PrimeField>;
+
+impl IsFFTField for Stark101PrimeField {
+    const TWO_ADICITY: u64 = 30;
+
+    const TWO_ADIC_PRIMITVE_ROOT_OF_UNITY: U64 = UnsignedInteger::from_hex_unchecked("bb6e79d");
+
+    fn field_name() -> &'static str {
+        "stark101"
+    }
+}
+
+impl FieldElement<Stark101PrimeField> {
+    pub fn to_bytes_le(&self) -> [u8; 8] {
+        let limbs = self.representative().limbs;
+        limbs[0].to_le_bytes()
+    }
+
+    pub fn to_bytes_be(&self) -> [u8; 8] {
+        let limbs = self.representative().limbs;
+        limbs[0].to_be_bytes()
+    }
+}
+
+#[allow(clippy::non_canonical_partial_ord_impl)]
+impl PartialOrd for FieldElement<Stark101PrimeField> {
+    fn partial_cmp(&self, other: &Self) -> Option<core::cmp::Ordering> {
+        self.representative().partial_cmp(&other.representative())
+    }
+}
+
+impl Ord for FieldElement<Stark101PrimeField> {
+    fn cmp(&self, other: &Self) -> core::cmp::Ordering {
+        self.representative().cmp(&other.representative())
+    }
+}
+
+#[cfg(test)]
+mod test_stark101_prime_field {
+    use super::Stark101PrimeField;
+    use crate::{
+        field::{element::FieldElement, traits::IsFFTField},
+        traits::ByteConversion,
+    };
+
+    #[test]
+    fn two_adic_order() {
+        let w = FieldElement::<Stark101PrimeField>::from(
+            Stark101PrimeField::TWO_ADIC_PRIMITVE_ROOT_OF_UNITY.limbs[0],
+        );
+
+        assert_eq!(
+            w.pow(1u64 << Stark101PrimeField::TWO_ADICITY),
+            FieldElement::one()
+        );
+        assert_ne!(
+            w.pow(1u64 << (Stark101PrimeField::TWO_ADICITY >> 1)),
+            FieldElement::one()
+        );
+    }
+
+    #[test]
+    #[cfg(feature = "alloc")]
+    fn byte_serialization_for_a_number_matches_with_byte_conversion_implementation_le() {
+        let element = FieldElement::<Stark101PrimeField>::from_hex_unchecked("0123456701234567");
+        let bytes = element.to_bytes_le();
+        let expected_bytes: [u8; 8] = ByteConversion::to_bytes_le(&element).try_into().unwrap();
+        assert_eq!(bytes, expected_bytes);
+    }
+
+    #[test]
+    #[cfg(feature = "alloc")]
+    fn byte_serialization_for_a_number_matches_with_byte_conversion_implementation_be() {
+        let element = FieldElement::<Stark101PrimeField>::from_hex_unchecked("0123456701234567");
+        let bytes = element.to_bytes_be();
+        let expected_bytes: [u8; 8] = ByteConversion::to_bytes_be(&element).try_into().unwrap();
+        assert_eq!(bytes, expected_bytes);
+    }
+
+    #[test]
+
+    fn byte_serialization_and_deserialization_works_le() {
+        let element = FieldElement::<Stark101PrimeField>::from_hex_unchecked("7654321076543210");
+        let bytes = element.to_bytes_le();
+        let from_bytes = FieldElement::<Stark101PrimeField>::from_bytes_le(&bytes).unwrap();
+        assert_eq!(element, from_bytes);
+    }
+
+    #[test]
+
+    fn byte_serialization_and_deserialization_works_be() {
+        let element = FieldElement::<Stark101PrimeField>::from_hex_unchecked("7654321076543210");
+        let bytes = element.to_bytes_be();
+        let from_bytes = FieldElement::<Stark101PrimeField>::from_bytes_be(&bytes).unwrap();
+        assert_eq!(element, from_bytes);
+    }
+}

From 7a427aacfd9241276c17f38ecbfd5e6dfae4e032 Mon Sep 17 00:00:00 2001
From: Nicole Graus <nicole.graus@lambdaclass.com>
Date: Tue, 1 Oct 2024 16:34:56 -0300
Subject: [PATCH 2/3] Optimize Mersenne31 Field (#921)

* optimize add

* save changes. Add, sub and mul checked

* fix tests

* add new inv

* add mult by powers of two

* replace inverse

* test new inv

* modify old algorithm for inv

* fix tests extension

* add mul for degree 4 extension

* add fp4 isField and isSubField operations and benchmarks

* new version for fp4 mul based on the paper

* add mul of a fp2e by non-residue

* change inv using mul_fp2_by_non_resiude

* save work

* wip fp2 test

* add fp2 tests

* add 2 * a^2 - 1 function

* use karatsuba in fp4 mul version 1

* clean up

* fix Fp as subfield of Fp2. Tests Fp plus Fp4 is now correct

* fix inv

* fix comments

* fix comments

* fixes

* fix clippy

* fix cargo check no-std

* fix clippy

* change zero function of isField to rust default

* fix two_square_minus_one function and optimize inv function

* fix clippy

---------

Co-authored-by: Nicole <nicole@Nicoles-MacBook-Air.local>
Co-authored-by: Joaquin Carletti <joaquin.carletti@lambdaclass.com>
Co-authored-by: diegokingston <dkingston@fi.uba.ar>
Co-authored-by: Diego K <43053772+diegokingston@users.noreply.github.com>
Co-authored-by: Joaquin Carletti <56092489+ColoCarletti@users.noreply.github.com>
---
 math/benches/criterion_field.rs               |   4 +-
 math/benches/fields/mersenne31.rs             |  64 +-
 math/src/field/fields/mersenne31/extension.rs | 304 ---------
 .../src/field/fields/mersenne31/extensions.rs | 627 ++++++++++++++++++
 math/src/field/fields/mersenne31/field.rs     | 272 ++++----
 math/src/field/fields/mersenne31/mod.rs       |   2 +-
 .../fields/p448_goldilocks_prime_field.rs     |   2 +-
 math/src/field/traits.rs                      |   8 +-
 math/src/unsigned_integer/element.rs          |   8 +
 9 files changed, 863 insertions(+), 428 deletions(-)
 delete mode 100644 math/src/field/fields/mersenne31/extension.rs
 create mode 100644 math/src/field/fields/mersenne31/extensions.rs

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/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/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)]

From 86b2789c95c9a02f1227b17b9019e1ad33f5290e Mon Sep 17 00:00:00 2001
From: jotabulacios <45471455+jotabulacios@users.noreply.github.com>
Date: Thu, 3 Oct 2024 12:48:10 -0300
Subject: [PATCH 3/3] Add Missing Components for BLS12-377 Curve (#922)

* save work. g1 subgroup check works but not the check for g2

* save work. both subgroup check works

* save work. curve is complete

* small refactor

* non-residue as a constant

* fix wrong constant

* fp2 residue as a constant

---------

Co-authored-by: Diego K <43053772+diegokingston@users.noreply.github.com>
---
 .../curves/bls12_377/curve.rs                 | 139 ++++++-
 .../curves/bls12_377/field_extension.rs       | 342 +++++++++++++++++-
 .../short_weierstrass/curves/bls12_377/mod.rs |   1 +
 .../curves/bls12_377/twist.rs                 | 172 +++++++++
 4 files changed, 640 insertions(+), 14 deletions(-)
 create mode 100644 math/src/elliptic_curve/short_weierstrass/curves/bls12_377/twist.rs

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);
+    }
+}