From 079eb67e917494a41d3bc85323ece1031d7b2528 Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Thu, 31 Mar 2022 23:17:53 +0200
Subject: [PATCH 1/7] Modular inversion using binary extended GCD

---
 math/Cargo.toml             |  9 +-------
 math/benches/inv.rs         | 42 +++++++++++++++++++++++++++++++++++++
 math/src/field/f64/mod.rs   | 34 ++++++++++++++++++++++++++++++
 math/src/field/f64/tests.rs | 11 ++++++++++
 4 files changed, 88 insertions(+), 8 deletions(-)
 create mode 100644 math/benches/inv.rs

diff --git a/math/Cargo.toml b/math/Cargo.toml
index f5cebc3a7..09b48e3d1 100644
--- a/math/Cargo.toml
+++ b/math/Cargo.toml
@@ -16,16 +16,9 @@ rust-version = "1.60"
 bench = false
 
 [[bench]]
-name = "fft"
+name = "inv"
 harness = false
 
-[[bench]]
-name = "field"
-harness = false
-
-[[bench]]
-name = "polynom"
-harness = false
 
 [features]
 concurrent = ["utils/concurrent", "std"]
diff --git a/math/benches/inv.rs b/math/benches/inv.rs
new file mode 100644
index 000000000..65b7d9b49
--- /dev/null
+++ b/math/benches/inv.rs
@@ -0,0 +1,42 @@
+// Copyright (c) Facebook, Inc. and its affiliates.
+//
+// This source code is licensed under the MIT license found in the
+// LICENSE file in the root directory of this source tree.
+
+use criterion::{criterion_group, criterion_main, Criterion};
+use rand_utils::rand_value;
+use winter_math::{
+    fields::f64::BaseElement as BaseElement64,
+    FieldElement,
+};
+
+// SEQUENTIAL OPS
+// ================================================================================================
+pub fn field_ops(c: &mut Criterion, field_name: &str) {
+    let mut group = c.benchmark_group(format!("op/{}", field_name));
+
+    group.bench_function("inv64", |bench| {
+        let x: BaseElement64 = rand_value();
+        bench.iter(|| BaseElement64::inv(x))
+    });
+
+    group.bench_function("inv64_gcd", |bench| {
+        let mut x: BaseElement64 = rand_value();
+        bench.iter(|| x.inv_gcd())
+    });
+
+
+}
+
+// GENERIC BENCHMARK RUNNER
+// ================================================================================================
+
+fn bench_field_ops(c: &mut Criterion) {
+    field_ops(c, "inversion");
+}
+
+// CRITERION BOILERPLATE
+// ================================================================================================
+
+criterion_group!(field_group, bench_field_ops);
+criterion_main!(field_group);
\ No newline at end of file
diff --git a/math/src/field/f64/mod.rs b/math/src/field/f64/mod.rs
index 287569fc3..43d836963 100644
--- a/math/src/field/f64/mod.rs
+++ b/math/src/field/f64/mod.rs
@@ -36,6 +36,9 @@ mod tests;
 // Field modulus = 2^64 - 2^32 + 1
 const M: u64 = 0xFFFFFFFF00000001;
 
+// (p+1)/2
+pub const MOD_2: u64 = (0xFFFFFFFF00000001 + 1u64) >> 1;
+
 // Epsilon = 2^32 - 1;
 const E: u64 = 0xFFFFFFFF;
 
@@ -60,6 +63,37 @@ impl BaseElement {
     pub const fn new(value: u64) -> Self {
         Self(value % M)
     }
+    pub const fn half(self) -> Self {
+        // If self is even, then return self/2.
+        // If self is odd, then return (self-1)/2 + (M+1)/2 = (x+M)/2.
+        Self((self.0 >> 1).wrapping_add((self.0 & 1).wrapping_neg() & MOD_2))
+    }
+
+    /// Return the inverse of "self" modulo the prime modulus.
+    /// Based on Algorithm1 in https://eprint.iacr.org/2020/972.pdf
+    pub fn inv_gcd(&mut self) -> BaseElement {
+        let mut a = self.0;
+        let mut u = BaseElement::from(1u64);
+        let mut b = BaseElement::MODULUS;
+        let mut v = BaseElement::from(0u64);
+
+        while a != 0 {
+            if a & 0x1 == 0x0 {
+                a >>= 1;
+                u = u.half();
+            } else {
+                if a < b {
+                    mem::swap(&mut a, &mut b);
+                    mem::swap(&mut u, &mut v);
+                }
+                a -= b;
+                a >>= 1;
+                u -= v;
+                u = u.half();
+            }
+        }
+        return v;
+    }
 }
 
 impl FieldElement for BaseElement {
diff --git a/math/src/field/f64/tests.rs b/math/src/field/f64/tests.rs
index edd15b630..bbfb08b67 100644
--- a/math/src/field/f64/tests.rs
+++ b/math/src/field/f64/tests.rs
@@ -113,6 +113,17 @@ fn inv() {
     assert_eq!(BaseElement::ZERO, BaseElement::inv(BaseElement::ZERO));
 }
 
+
+#[test]
+fn inv_gcd() {
+    assert_eq!(BaseElement::ONE, (BaseElement::from(1u64)).inv_gcd());
+    assert_eq!(BaseElement::ZERO, (BaseElement::from(0u64)).inv_gcd());
+    assert_eq!(
+        BaseElement::inv(BaseElement::from(4232346u64).inv_gcd()),
+        BaseElement::from(4232346u64)
+    );
+}
+
 #[test]
 fn element_as_int() {
     let v = u64::MAX;

From b4ae41db66880a52a7e2f2b5259f5eb240f83286 Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Mon, 13 Jun 2022 23:06:31 +0200
Subject: [PATCH 2/7] feat: Montgomery rep. based 64bit field

---
 math/Cargo.toml                        |   2 +-
 math/benches/inv.rs                    |  42 ----
 math/src/field/extensions/quadratic.rs |   2 +-
 math/src/field/f64/mod.rs              | 275 +++++++++++--------------
 math/src/field/f64/tests.rs            |  27 +--
 5 files changed, 134 insertions(+), 214 deletions(-)
 delete mode 100644 math/benches/inv.rs

diff --git a/math/Cargo.toml b/math/Cargo.toml
index 09b48e3d1..a8b4ee259 100644
--- a/math/Cargo.toml
+++ b/math/Cargo.toml
@@ -16,7 +16,7 @@ rust-version = "1.60"
 bench = false
 
 [[bench]]
-name = "inv"
+name = "fft"
 harness = false
 
 
diff --git a/math/benches/inv.rs b/math/benches/inv.rs
deleted file mode 100644
index 65b7d9b49..000000000
--- a/math/benches/inv.rs
+++ /dev/null
@@ -1,42 +0,0 @@
-// Copyright (c) Facebook, Inc. and its affiliates.
-//
-// This source code is licensed under the MIT license found in the
-// LICENSE file in the root directory of this source tree.
-
-use criterion::{criterion_group, criterion_main, Criterion};
-use rand_utils::rand_value;
-use winter_math::{
-    fields::f64::BaseElement as BaseElement64,
-    FieldElement,
-};
-
-// SEQUENTIAL OPS
-// ================================================================================================
-pub fn field_ops(c: &mut Criterion, field_name: &str) {
-    let mut group = c.benchmark_group(format!("op/{}", field_name));
-
-    group.bench_function("inv64", |bench| {
-        let x: BaseElement64 = rand_value();
-        bench.iter(|| BaseElement64::inv(x))
-    });
-
-    group.bench_function("inv64_gcd", |bench| {
-        let mut x: BaseElement64 = rand_value();
-        bench.iter(|| x.inv_gcd())
-    });
-
-
-}
-
-// GENERIC BENCHMARK RUNNER
-// ================================================================================================
-
-fn bench_field_ops(c: &mut Criterion) {
-    field_ops(c, "inversion");
-}
-
-// CRITERION BOILERPLATE
-// ================================================================================================
-
-criterion_group!(field_group, bench_field_ops);
-criterion_main!(field_group);
\ No newline at end of file
diff --git a/math/src/field/extensions/quadratic.rs b/math/src/field/extensions/quadratic.rs
index e1e8f306c..5eee9a09b 100644
--- a/math/src/field/extensions/quadratic.rs
+++ b/math/src/field/extensions/quadratic.rs
@@ -325,7 +325,7 @@ impl<B: ExtensibleField<2>> Deserializable for QuadExtension<B> {
 #[cfg(test)]
 mod tests {
     use super::{DeserializationError, FieldElement, QuadExtension, Vec};
-    use crate::field::f128::BaseElement;
+    use crate::field::f64::BaseElement;
     use rand_utils::rand_value;
 
     // BASIC ALGEBRA
diff --git a/math/src/field/f64/mod.rs b/math/src/field/f64/mod.rs
index 43d836963..11d8477ea 100644
--- a/math/src/field/f64/mod.rs
+++ b/math/src/field/f64/mod.rs
@@ -36,11 +36,14 @@ mod tests;
 // Field modulus = 2^64 - 2^32 + 1
 const M: u64 = 0xFFFFFFFF00000001;
 
+/// 2^128 mod M; this is used for conversion of elements into Montgomery representation.
+const R2: u64 = 0xFFFFFFFE00000001;
+
 // (p+1)/2
 pub const MOD_2: u64 = (0xFFFFFFFF00000001 + 1u64) >> 1;
 
 // Epsilon = 2^32 - 1;
-const E: u64 = 0xFFFFFFFF;
+//const E: u64 = 0xFFFFFFFF;
 
 // 2^32 root of unity
 const G: u64 = 1753635133440165772;
@@ -56,43 +59,82 @@ const ELEMENT_BYTES: usize = core::mem::size_of::<u64>();
 /// Internal values are stored in the range [0, 2^64). The backing type is `u64`.
 #[derive(Copy, Clone, Debug, Default)]
 pub struct BaseElement(u64);
-
 impl BaseElement {
-    /// Creates a new field element from the provided `value`. If the value is greater than or
-    /// equal to the field modulus, modular reduction is silently performed.
-    pub const fn new(value: u64) -> Self {
-        Self(value % M)
-    }
-    pub const fn half(self) -> Self {
-        // If self is even, then return self/2.
-        // If self is odd, then return (self-1)/2 + (M+1)/2 = (x+M)/2.
-        Self((self.0 >> 1).wrapping_add((self.0 & 1).wrapping_neg() & MOD_2))
-    }
-
-    /// Return the inverse of "self" modulo the prime modulus.
-    /// Based on Algorithm1 in https://eprint.iacr.org/2020/972.pdf
-    pub fn inv_gcd(&mut self) -> BaseElement {
-        let mut a = self.0;
-        let mut u = BaseElement::from(1u64);
-        let mut b = BaseElement::MODULUS;
-        let mut v = BaseElement::from(0u64);
-
-        while a != 0 {
-            if a & 0x1 == 0x0 {
-                a >>= 1;
-                u = u.half();
-            } else {
-                if a < b {
-                    mem::swap(&mut a, &mut b);
-                    mem::swap(&mut u, &mut v);
-                }
-                a -= b;
-                a >>= 1;
-                u -= v;
-                u = u.half();
-            }
+    /// Creates a new field element from the provided `value`; the value is converted into
+    /// Montgomery representation.
+    pub const fn new(value: u64) -> BaseElement {
+        BaseElement(BaseElement::mont_red((value as u128) * (R2 as u128)))
+    }
+    /// Gets the inner value that might not be canonical
+
+    pub const fn inner(self: &Self) -> u64 {
+        return self.0;
+    }
+
+    /// Montgomery reduction
+    #[inline(always)]
+    pub const fn mont_red(x: u128) -> u64 {
+        // See reference above for a description of the following implementation.
+        let xl = x as u64;
+        let xh = (x >> 64) as u64;
+        let (a, e) = xl.overflowing_add(xl << 32);
+
+        let b = a.wrapping_sub(a >> 32).wrapping_sub(e as u64);
+
+        let (r, c) = xh.overflowing_sub(b);
+        r.wrapping_sub(0u32.wrapping_sub(c as u32) as u64)
+    }
+
+    /// Addition in BaseField
+    #[inline(always)]
+    const fn add(self, rhs: Self) -> Self {
+        // We compute a + b = a - (p - b).
+        let (x1, c1) = self.0.overflowing_sub(M - rhs.0);
+        let adj = 0u32.wrapping_sub(c1 as u32);
+        BaseElement(x1.wrapping_sub(adj as u64))
+    }
+
+    /// Subtraction in BaseField
+    #[inline(always)]
+    const fn sub(self, rhs: Self) -> Self {
+        // See reference above for more details.
+        let (x1, c1) = self.0.overflowing_sub(rhs.0);
+        let adj = 0u32.wrapping_sub(c1 as u32);
+        BaseElement(x1.wrapping_sub(adj as u64))
+    }
+
+    /// Multiplication in BaseField
+    #[inline(always)]
+    const fn mul(self, rhs: Self) -> Self {
+        // If x < p and y < p, then x*y <= (p-1)^2, and is thus in
+        // range of mont_red().
+        BaseElement(BaseElement::mont_red((self.0 as u128) * (rhs.0 as u128)))
+    }
+
+    /// Squaring in BaseField
+    #[inline(always)]
+    pub const fn square(self) -> Self {
+        self.mul(self)
+    }
+
+    /// Multiple squarings in BaseField: return x^(2^n)
+    pub fn msquare(self, n: u32) -> Self {
+        let mut x = self;
+        for _ in 0..n {
+            x = x.square();
         }
-        return v;
+        x
+    }
+
+    /// Test of equality between two BaseField elements; return value is
+    /// 0xFFFFFFFFFFFFFFFF if the two values are equal, or 0 otherwise.
+    #[inline(always)]
+    pub const fn equals(self, rhs: Self) -> u64 {
+        // Since internal representation is canonical, we can simply
+        // do a xor between the two operands, and then use the same
+        // expression as iszero().
+        let t = self.0 ^ rhs.0;
+        !((((t | t.wrapping_neg()) as i64) >> 63) as u64)
     }
 }
 
@@ -100,13 +142,12 @@ impl FieldElement for BaseElement {
     type PositiveInteger = u64;
     type BaseField = Self;
 
-    const ZERO: Self = Self::new(0);
-    const ONE: Self = Self::new(1);
+    const ZERO: Self = BaseElement::new(0);
+    const ONE: Self = BaseElement::new(1);
 
     const ELEMENT_BYTES: usize = ELEMENT_BYTES;
     const IS_CANONICAL: bool = false;
 
-    #[inline]
     fn exp(self, power: Self::PositiveInteger) -> Self {
         let mut b = self;
 
@@ -127,40 +168,24 @@ impl FieldElement for BaseElement {
         r
     }
 
-    #[inline]
-    #[allow(clippy::many_single_char_names)]
     fn inv(self) -> Self {
-        // compute base^(M - 2) using 72 multiplications
-        // M - 2 = 0b1111111111111111111111111111111011111111111111111111111111111111
-
-        // compute base^11
-        let t2 = self.square() * self;
-
-        // compute base^111
-        let t3 = t2.square() * self;
-
-        // compute base^111111 (6 ones)
-        let t6 = exp_acc::<3>(t3, t3);
-
-        // compute base^111111111111 (12 ones)
-        let t12 = exp_acc::<6>(t6, t6);
-
-        // compute base^111111111111111111111111 (24 ones)
-        let t24 = exp_acc::<12>(t12, t12);
-
-        // compute base^1111111111111111111111111111111 (31 ones)
-        let t30 = exp_acc::<6>(t24, t6);
-        let t31 = t30.square() * self;
-
-        // compute base^111111111111111111111111111111101111111111111111111111111111111
-        let t63 = exp_acc::<32>(t31, t31);
-
-        // compute base^1111111111111111111111111111111011111111111111111111111111111111
-        t63.square() * self
+        // This uses Fermat's little theorem: 1/x = x^(p-2) mod p.
+        // We have p-2 = 0xFFFFFFFEFFFFFFFF. In the instructions below,
+        // we call 'xj' the value x^(2^j-1).
+        let x = self;
+        let x2 = x * x.square();
+        let x4 = x2 * x2.msquare(2);
+        let x5 = x * x4.square();
+        let x10 = x5 * x5.msquare(5);
+        let x15 = x5 * x10.msquare(5);
+        let x16 = x * x15.square();
+        let x31 = x15 * x16.msquare(15);
+        let x32 = x * x31.square();
+        return x32 * x31.msquare(33);
     }
 
     fn conjugate(&self) -> Self {
-        Self(self.0)
+        BaseElement(self.0)
     }
 
     fn elements_as_bytes(elements: &[Self]) -> &[u8] {
@@ -205,7 +230,6 @@ impl FieldElement for BaseElement {
         unsafe { Vec::from_raw_parts(p as *mut Self, len, cap) }
     }
 
-    #[inline]
     fn as_base_elements(elements: &[Self]) -> &[Self::BaseField] {
         elements
     }
@@ -234,18 +258,12 @@ impl StarkField for BaseElement {
     const TWO_ADIC_ROOT_OF_UNITY: Self = Self::new(G);
 
     fn get_modulus_le_bytes() -> Vec<u8> {
-        M.to_le_bytes().to_vec()
+        Self::MODULUS.to_le_bytes().to_vec()
     }
 
     #[inline]
     fn as_int(&self) -> Self::PositiveInteger {
-        // since the internal value of the element can be in [0, 2^64) range, we do an extra check
-        // here to convert it to the canonical form
-        if self.0 >= M {
-            self.0 - M
-        } else {
-            self.0
-        }
+        BaseElement::mont_red(self.0 as u128)
     }
 }
 
@@ -269,9 +287,7 @@ impl Display for BaseElement {
 impl PartialEq for BaseElement {
     #[inline]
     fn eq(&self, other: &Self) -> bool {
-        // since either of the elements can be in [0, 2^64) range, we first convert them to the
-        // canonical form to ensure that they are in [0, M) range and then compare them.
-        self.as_int() == other.as_int()
+        Self::equals(*self, *other) == 0xFFFFFFFFFFFFFFFF
     }
 }
 
@@ -283,16 +299,12 @@ impl Eq for BaseElement {}
 impl Add for BaseElement {
     type Output = Self;
 
-    #[inline]
-    #[allow(clippy::suspicious_arithmetic_impl)]
     fn add(self, rhs: Self) -> Self {
-        let (result, over) = self.0.overflowing_add(rhs.as_int());
-        Self(result.wrapping_sub(M * (over as u64)))
+        Self::add(self, rhs)
     }
 }
 
 impl AddAssign for BaseElement {
-    #[inline]
     fn add_assign(&mut self, rhs: Self) {
         *self = *self + rhs
     }
@@ -301,16 +313,12 @@ impl AddAssign for BaseElement {
 impl Sub for BaseElement {
     type Output = Self;
 
-    #[inline]
-    #[allow(clippy::suspicious_arithmetic_impl)]
     fn sub(self, rhs: Self) -> Self {
-        let (result, under) = self.0.overflowing_sub(rhs.as_int());
-        Self(result.wrapping_add(M * (under as u64)))
+        Self::sub(self, rhs)
     }
 }
 
 impl SubAssign for BaseElement {
-    #[inline]
     fn sub_assign(&mut self, rhs: Self) {
         *self = *self - rhs;
     }
@@ -319,15 +327,12 @@ impl SubAssign for BaseElement {
 impl Mul for BaseElement {
     type Output = Self;
 
-    #[inline]
     fn mul(self, rhs: Self) -> Self {
-        let z = (self.0 as u128) * (rhs.0 as u128);
-        Self(mod_reduce(z))
+        Self::mul(self, rhs)
     }
 }
 
 impl MulAssign for BaseElement {
-    #[inline]
     fn mul_assign(&mut self, rhs: Self) {
         *self = *self * rhs
     }
@@ -336,15 +341,12 @@ impl MulAssign for BaseElement {
 impl Div for BaseElement {
     type Output = Self;
 
-    #[inline]
-    #[allow(clippy::suspicious_arithmetic_impl)]
     fn div(self, rhs: Self) -> Self {
-        self * rhs.inv()
+        Self::mul(self, Self::inv(rhs))
     }
 }
 
 impl DivAssign for BaseElement {
-    #[inline]
     fn div_assign(&mut self, rhs: Self) {
         *self = *self / rhs
     }
@@ -353,14 +355,8 @@ impl DivAssign for BaseElement {
 impl Neg for BaseElement {
     type Output = Self;
 
-    #[inline]
     fn neg(self) -> Self {
-        let v = self.as_int();
-        if v == 0 {
-            Self::ZERO
-        } else {
-            Self(M - v)
-        }
+        Self::sub(BaseElement::ZERO, self)
     }
 }
 
@@ -455,10 +451,11 @@ impl ExtensibleField<3> for BaseElement {
 // ================================================================================================
 
 impl From<u128> for BaseElement {
-    /// Converts a 128-bit value into a field element. If the value is greater than or equal to
-    /// the field modulus, modular reduction is silently performed.
-    fn from(value: u128) -> Self {
-        Self(mod_reduce(value))
+    /// Converts a 128-bit value into a field element.
+    fn from(x: u128) -> Self {
+        //const R3: u128 = 1 (= 2^192 mod M );// thus we get that mont_reduce((mont_reduce(x) as u128) * R3) becomes
+        //BaseElement(mont_reduce(mont_reduce(x) as u128))  // With branching
+        BaseElement(Self::mont_red(Self::mont_red(x) as u128)) // Constant time
     }
 }
 
@@ -466,28 +463,28 @@ impl From<u64> for BaseElement {
     /// Converts a 64-bit value into a field element. If the value is greater than or equal to
     /// the field modulus, modular reduction is silently performed.
     fn from(value: u64) -> Self {
-        Self::new(value)
+        BaseElement::new(value)
     }
 }
 
 impl From<u32> for BaseElement {
     /// Converts a 32-bit value into a field element.
     fn from(value: u32) -> Self {
-        Self::new(value as u64)
+        BaseElement::new(value as u64)
     }
 }
 
 impl From<u16> for BaseElement {
     /// Converts a 16-bit value into a field element.
     fn from(value: u16) -> Self {
-        Self::new(value as u64)
+        BaseElement::new(value as u64)
     }
 }
 
 impl From<u8> for BaseElement {
     /// Converts an 8-bit value into a field element.
     fn from(value: u8) -> Self {
-        Self::new(value as u64)
+        BaseElement::new(value as u64)
     }
 }
 
@@ -498,7 +495,7 @@ impl From<[u8; 8]> for BaseElement {
     /// performed.
     fn from(bytes: [u8; 8]) -> Self {
         let value = u64::from_le_bytes(bytes);
-        Self::new(value)
+        BaseElement::new(value)
     }
 }
 
@@ -568,41 +565,15 @@ impl Deserializable for BaseElement {
     }
 }
 
-// HELPER FUNCTIONS
-// ================================================================================================
-
-/// Reduces a 128-bit value by M such that the output is in [0, 2^64) range.
-///
-/// Adapted from: <https://github.com/mir-protocol/plonky2/blob/main/src/field/goldilocks_field.rs>
-#[inline(always)]
-fn mod_reduce(x: u128) -> u64 {
-    // assume x consists of four 32-bit values: a, b, c, d such that a contains 32 least
-    // significant bits and d contains 32 most significant bits. we break x into corresponding
-    // values as shown below
-    let ab = x as u64;
-    let cd = (x >> 64) as u64;
-    let c = (cd as u32) as u64;
-    let d = cd >> 32;
-
-    // compute ab - d; because d may be greater than ab we need to handle potential underflow
-    let (tmp0, under) = ab.overflowing_sub(d);
-    let tmp0 = tmp0.wrapping_sub(E * (under as u64));
-
-    // compute c * 2^32 - c; this is guaranteed not to underflow
-    let tmp1 = (c << 32) - c;
-
-    // add temp values and return the result; because each of the temp may be up to 64 bits,
-    // we need to handle potential overflow
-    let (result, over) = tmp0.overflowing_add(tmp1);
-    result.wrapping_add(E * (over as u64))
-}
-
-/// Squares the base N number of times and multiplies the result by the tail value.
 #[inline(always)]
-fn exp_acc<const N: usize>(base: BaseElement, tail: BaseElement) -> BaseElement {
-    let mut result = base;
-    for _ in 0..N {
-        result = result.square();
-    }
-    result * tail
-}
+pub fn mont_reduce(x: u128) -> u64 {
+    const NPRIME: u64 = 4294967297;
+    let q = (((x as u64) as u128) * (NPRIME as u128)) as u64;
+    let m = (q as u128) * (M as u128);
+    let y = (((x as i128).wrapping_sub(m as i128)) >> 64) as i64;
+    if x < m {
+        return (y + (M as i64)) as u64;
+    } else {
+        return y as u64;
+    };
+}
\ No newline at end of file
diff --git a/math/src/field/f64/tests.rs b/math/src/field/f64/tests.rs
index bbfb08b67..033cd2002 100644
--- a/math/src/field/f64/tests.rs
+++ b/math/src/field/f64/tests.rs
@@ -4,7 +4,8 @@
 // LICENSE file in the root directory of this source tree.
 
 use super::{
-    AsBytes, BaseElement, DeserializationError, FieldElement, Serializable, StarkField, E, M,
+    //E,
+    AsBytes, BaseElement, DeserializationError, FieldElement, Serializable, StarkField, M,
 };
 use crate::field::{CubeExtension, ExtensionOf, QuadExtension};
 use core::convert::TryFrom;
@@ -32,10 +33,10 @@ fn add() {
     assert_eq!(BaseElement::ZERO, t + BaseElement::ONE);
     assert_eq!(BaseElement::ONE, t + BaseElement::new(2));
 
-    // test non-canonical representation
-    let a = BaseElement::new(M - 1) + BaseElement::new(E);
-    let expected = ((((M - 1 + E) as u128) * 2) % (M as u128)) as u64;
-    assert_eq!(expected, (a + a).as_int());
+    //// test non-canonical representation
+    //let a = BaseElement::new(M - 1) + BaseElement::new(E);
+    //let expected = ((((M - 1 + E) as u128) * 2) % (M as u128)) as u64;
+    //assert_eq!(expected, (a + a).as_int());
 }
 
 #[test]
@@ -71,7 +72,7 @@ fn mul() {
     assert_eq!(BaseElement::ZERO, r * BaseElement::ZERO);
     assert_eq!(r, r * BaseElement::ONE);
 
-    // test multiplication within bounds
+    // test multifield::extensions::cubic::tests::bytes_as_elementsplication within bounds
     assert_eq!(
         BaseElement::from(15u8),
         BaseElement::from(5u8) * BaseElement::from(3u8)
@@ -114,16 +115,6 @@ fn inv() {
 }
 
 
-#[test]
-fn inv_gcd() {
-    assert_eq!(BaseElement::ONE, (BaseElement::from(1u64)).inv_gcd());
-    assert_eq!(BaseElement::ZERO, (BaseElement::from(0u64)).inv_gcd());
-    assert_eq!(
-        BaseElement::inv(BaseElement::from(4232346u64).inv_gcd()),
-        BaseElement::from(4232346u64)
-    );
-}
-
 #[test]
 fn element_as_int() {
     let v = u64::MAX;
@@ -142,8 +133,8 @@ fn equals() {
     assert_eq!(a.to_bytes(), b.to_bytes());
 
     // but their internal representation is not
-    assert_ne!(a.0, b.0);
-    assert_ne!(a.as_bytes(), b.as_bytes());
+    //assert_ne!(a.0, b.0);
+    //assert_ne!(a.as_bytes(), b.as_bytes());
 }
 
 // ROOTS OF UNITY

From 06d538ccadee529142fe4e47e8dbe047877a015c Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Tue, 14 Jun 2022 21:42:30 +0200
Subject: [PATCH 3/7] feat: Montgomery rep. based 64bit field

---
 math/Cargo.toml           |   8 ++
 math/src/field/f64/mod.rs | 169 ++++++++++++++++++++------------------
 2 files changed, 98 insertions(+), 79 deletions(-)

diff --git a/math/Cargo.toml b/math/Cargo.toml
index a8b4ee259..f3a8b3350 100644
--- a/math/Cargo.toml
+++ b/math/Cargo.toml
@@ -19,6 +19,14 @@ bench = false
 name = "fft"
 harness = false
 
+[[bench]]
+name = "field"
+harness = false
+
+[[bench]]
+name = "polynom"
+harness = false
+
 
 [features]
 concurrent = ["utils/concurrent", "std"]
diff --git a/math/src/field/f64/mod.rs b/math/src/field/f64/mod.rs
index 11d8477ea..027c877b1 100644
--- a/math/src/field/f64/mod.rs
+++ b/math/src/field/f64/mod.rs
@@ -3,7 +3,8 @@
 // This source code is licensed under the MIT license found in the
 // LICENSE file in the root directory of this source tree.
 
-//! An implementation of a 64-bit STARK-friendly prime field with modulus $2^{64} - 2^{32} + 1$.
+//! An implementation of a 64-bit STARK-friendly prime field with modulus $2^{64} - 2^{32} + 1$ using Montgomery representation.
+//! Our implementation follows https://eprint.iacr.org/2022/274.pdf
 //!
 //! This field supports very fast modular arithmetic and has a number of other attractive
 //! properties, including:
@@ -42,9 +43,6 @@ const R2: u64 = 0xFFFFFFFE00000001;
 // (p+1)/2
 pub const MOD_2: u64 = (0xFFFFFFFF00000001 + 1u64) >> 1;
 
-// Epsilon = 2^32 - 1;
-//const E: u64 = 0xFFFFFFFF;
-
 // 2^32 root of unity
 const G: u64 = 1753635133440165772;
 
@@ -63,7 +61,7 @@ impl BaseElement {
     /// Creates a new field element from the provided `value`; the value is converted into
     /// Montgomery representation.
     pub const fn new(value: u64) -> BaseElement {
-        BaseElement(BaseElement::mont_red((value as u128) * (R2 as u128)))
+        BaseElement(mont_red_cst((value as u128) * (R2 as u128)))
     }
     /// Gets the inner value that might not be canonical
 
@@ -71,52 +69,6 @@ impl BaseElement {
         return self.0;
     }
 
-    /// Montgomery reduction
-    #[inline(always)]
-    pub const fn mont_red(x: u128) -> u64 {
-        // See reference above for a description of the following implementation.
-        let xl = x as u64;
-        let xh = (x >> 64) as u64;
-        let (a, e) = xl.overflowing_add(xl << 32);
-
-        let b = a.wrapping_sub(a >> 32).wrapping_sub(e as u64);
-
-        let (r, c) = xh.overflowing_sub(b);
-        r.wrapping_sub(0u32.wrapping_sub(c as u32) as u64)
-    }
-
-    /// Addition in BaseField
-    #[inline(always)]
-    const fn add(self, rhs: Self) -> Self {
-        // We compute a + b = a - (p - b).
-        let (x1, c1) = self.0.overflowing_sub(M - rhs.0);
-        let adj = 0u32.wrapping_sub(c1 as u32);
-        BaseElement(x1.wrapping_sub(adj as u64))
-    }
-
-    /// Subtraction in BaseField
-    #[inline(always)]
-    const fn sub(self, rhs: Self) -> Self {
-        // See reference above for more details.
-        let (x1, c1) = self.0.overflowing_sub(rhs.0);
-        let adj = 0u32.wrapping_sub(c1 as u32);
-        BaseElement(x1.wrapping_sub(adj as u64))
-    }
-
-    /// Multiplication in BaseField
-    #[inline(always)]
-    const fn mul(self, rhs: Self) -> Self {
-        // If x < p and y < p, then x*y <= (p-1)^2, and is thus in
-        // range of mont_red().
-        BaseElement(BaseElement::mont_red((self.0 as u128) * (rhs.0 as u128)))
-    }
-
-    /// Squaring in BaseField
-    #[inline(always)]
-    pub const fn square(self) -> Self {
-        self.mul(self)
-    }
-
     /// Multiple squarings in BaseField: return x^(2^n)
     pub fn msquare(self, n: u32) -> Self {
         let mut x = self;
@@ -142,12 +94,13 @@ impl FieldElement for BaseElement {
     type PositiveInteger = u64;
     type BaseField = Self;
 
-    const ZERO: Self = BaseElement::new(0);
-    const ONE: Self = BaseElement::new(1);
+    const ZERO: Self = Self::new(0);
+    const ONE: Self = Self::new(1);
 
     const ELEMENT_BYTES: usize = ELEMENT_BYTES;
     const IS_CANONICAL: bool = false;
 
+    #[inline(always)]
     fn exp(self, power: Self::PositiveInteger) -> Self {
         let mut b = self;
 
@@ -168,24 +121,40 @@ impl FieldElement for BaseElement {
         r
     }
 
+    #[inline]
+    #[allow(clippy::many_single_char_names)]
     fn inv(self) -> Self {
-        // This uses Fermat's little theorem: 1/x = x^(p-2) mod p.
-        // We have p-2 = 0xFFFFFFFEFFFFFFFF. In the instructions below,
-        // we call 'xj' the value x^(2^j-1).
-        let x = self;
-        let x2 = x * x.square();
-        let x4 = x2 * x2.msquare(2);
-        let x5 = x * x4.square();
-        let x10 = x5 * x5.msquare(5);
-        let x15 = x5 * x10.msquare(5);
-        let x16 = x * x15.square();
-        let x31 = x15 * x16.msquare(15);
-        let x32 = x * x31.square();
-        return x32 * x31.msquare(33);
+        // compute base^(M - 2) using 72 multiplications
+        // M - 2 = 0b1111111111111111111111111111111011111111111111111111111111111111
+
+        // compute base^11
+        let t2 = self.square() * self;
+
+        // compute base^111
+        let t3 = t2.square() * self;
+
+        // compute base^111111 (6 ones)
+        let t6 = exp_acc::<3>(t3, t3);
+
+        // compute base^111111111111 (12 ones)
+        let t12 = exp_acc::<6>(t6, t6);
+
+        // compute base^111111111111111111111111 (24 ones)
+        let t24 = exp_acc::<12>(t12, t12);
+
+        // compute base^1111111111111111111111111111111 (31 ones)
+        let t30 = exp_acc::<6>(t24, t6);
+        let t31 = t30.square() * self;
+
+        // compute base^111111111111111111111111111111101111111111111111111111111111111
+        let t63 = exp_acc::<32>(t31, t31);
+
+        // compute base^1111111111111111111111111111111011111111111111111111111111111111
+        t63.square() * self
     }
 
     fn conjugate(&self) -> Self {
-        BaseElement(self.0)
+        Self(self.0)
     }
 
     fn elements_as_bytes(elements: &[Self]) -> &[u8] {
@@ -258,12 +227,12 @@ impl StarkField for BaseElement {
     const TWO_ADIC_ROOT_OF_UNITY: Self = Self::new(G);
 
     fn get_modulus_le_bytes() -> Vec<u8> {
-        Self::MODULUS.to_le_bytes().to_vec()
+        M.to_le_bytes().to_vec()
     }
 
     #[inline]
     fn as_int(&self) -> Self::PositiveInteger {
-        BaseElement::mont_red(self.0 as u128)
+        mont_red_cst(self.0 as u128)
     }
 }
 
@@ -299,12 +268,19 @@ impl Eq for BaseElement {}
 impl Add for BaseElement {
     type Output = Self;
 
+    /// Addition in BaseField
+    #[inline]
+    #[allow(clippy::suspicious_arithmetic_impl)]
     fn add(self, rhs: Self) -> Self {
-        Self::add(self, rhs)
+        // We compute a + b = a - (p - b).
+        let (x1, c1) = self.0.overflowing_sub(M - rhs.0);
+        let adj = 0u32.wrapping_sub(c1 as u32);
+        BaseElement(x1.wrapping_sub(adj as u64))
     }
 }
 
 impl AddAssign for BaseElement {
+    #[inline]
     fn add_assign(&mut self, rhs: Self) {
         *self = *self + rhs
     }
@@ -313,12 +289,18 @@ impl AddAssign for BaseElement {
 impl Sub for BaseElement {
     type Output = Self;
 
+    #[inline]
+    #[allow(clippy::suspicious_arithmetic_impl)]
     fn sub(self, rhs: Self) -> Self {
-        Self::sub(self, rhs)
+        // See reference above for more details.
+        let (x1, c1) = self.0.overflowing_sub(rhs.0);
+        let adj = 0u32.wrapping_sub(c1 as u32);
+        BaseElement(x1.wrapping_sub(adj as u64))
     }
 }
 
 impl SubAssign for BaseElement {
+    #[inline]
     fn sub_assign(&mut self, rhs: Self) {
         *self = *self - rhs;
     }
@@ -327,12 +309,14 @@ impl SubAssign for BaseElement {
 impl Mul for BaseElement {
     type Output = Self;
 
+    #[inline]
     fn mul(self, rhs: Self) -> Self {
-        Self::mul(self, rhs)
+        BaseElement(mont_red_cst((self.0 as u128) * (rhs.0 as u128)))
     }
 }
 
 impl MulAssign for BaseElement {
+    #[inline]
     fn mul_assign(&mut self, rhs: Self) {
         *self = *self * rhs
     }
@@ -341,12 +325,15 @@ impl MulAssign for BaseElement {
 impl Div for BaseElement {
     type Output = Self;
 
+    #[inline]
+    #[allow(clippy::suspicious_arithmetic_impl)]
     fn div(self, rhs: Self) -> Self {
-        Self::mul(self, Self::inv(rhs))
+        self * Self::inv(rhs)
     }
 }
 
 impl DivAssign for BaseElement {
+    #[inline]
     fn div_assign(&mut self, rhs: Self) {
         *self = *self / rhs
     }
@@ -355,8 +342,9 @@ impl DivAssign for BaseElement {
 impl Neg for BaseElement {
     type Output = Self;
 
+    #[inline]
     fn neg(self) -> Self {
-        Self::sub(BaseElement::ZERO, self)
+        BaseElement::ZERO - self
     }
 }
 
@@ -453,9 +441,9 @@ impl ExtensibleField<3> for BaseElement {
 impl From<u128> for BaseElement {
     /// Converts a 128-bit value into a field element.
     fn from(x: u128) -> Self {
-        //const R3: u128 = 1 (= 2^192 mod M );// thus we get that mont_reduce((mont_reduce(x) as u128) * R3) becomes
-        //BaseElement(mont_reduce(mont_reduce(x) as u128))  // With branching
-        BaseElement(Self::mont_red(Self::mont_red(x) as u128)) // Constant time
+        //const R3: u128 = 1 (= 2^192 mod M );// thus we get that mont_red_var((mont_red_var(x) as u128) * R3) becomes
+        //BaseElement(mont_red_var(mont_red_var(x) as u128))  // Variable time implementation
+        BaseElement(mont_red_cst(mont_red_cst(x) as u128)) // Constant time implementation
     }
 }
 
@@ -565,8 +553,18 @@ impl Deserializable for BaseElement {
     }
 }
 
+/// Squares the base N number of times and multiplies the result by the tail value.
 #[inline(always)]
-pub fn mont_reduce(x: u128) -> u64 {
+fn exp_acc<const N: usize>(base: BaseElement, tail: BaseElement) -> BaseElement {
+    let mut result = base;
+    for _ in 0..N {
+        result = result.square();
+    }
+    result * tail
+}
+/// Montgomery reduction (variable time)
+#[inline(always)]
+pub fn mont_red_var(x: u128) -> u64 {
     const NPRIME: u64 = 4294967297;
     let q = (((x as u64) as u128) * (NPRIME as u128)) as u64;
     let m = (q as u128) * (M as u128);
@@ -576,4 +574,17 @@ pub fn mont_reduce(x: u128) -> u64 {
     } else {
         return y as u64;
     };
-}
\ No newline at end of file
+}
+/// Montgomery reduction (constant time)
+#[inline(always)]
+pub const fn mont_red_cst(x: u128) -> u64 {
+    // See reference above for a description of the following implementation.
+    let xl = x as u64;
+    let xh = (x >> 64) as u64;
+    let (a, e) = xl.overflowing_add(xl << 32);
+
+    let b = a.wrapping_sub(a >> 32).wrapping_sub(e as u64);
+
+    let (r, c) = xh.overflowing_sub(b);
+    r.wrapping_sub(0u32.wrapping_sub(c as u32) as u64)
+}

From 8d47e96c0a89619e0739e6f09338f896b4552991 Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Wed, 15 Jun 2022 01:27:00 +0200
Subject: [PATCH 4/7] feat: Montgomery rep. based 64bit field

---
 math/Cargo.toml           |  1 -
 math/src/field/f64/mod.rs | 89 ++++++++++++++++-----------------------
 2 files changed, 36 insertions(+), 54 deletions(-)

diff --git a/math/Cargo.toml b/math/Cargo.toml
index f3a8b3350..f5cebc3a7 100644
--- a/math/Cargo.toml
+++ b/math/Cargo.toml
@@ -27,7 +27,6 @@ harness = false
 name = "polynom"
 harness = false
 
-
 [features]
 concurrent = ["utils/concurrent", "std"]
 default = ["std"]
diff --git a/math/src/field/f64/mod.rs b/math/src/field/f64/mod.rs
index 027c877b1..9d92f4d8d 100644
--- a/math/src/field/f64/mod.rs
+++ b/math/src/field/f64/mod.rs
@@ -3,7 +3,8 @@
 // This source code is licensed under the MIT license found in the
 // LICENSE file in the root directory of this source tree.
 
-//! An implementation of a 64-bit STARK-friendly prime field with modulus $2^{64} - 2^{32} + 1$ using Montgomery representation.
+//! An implementation of a 64-bit STARK-friendly prime field with modulus $2^{64} - 2^{32} + 1$
+//! using Montgomery representation.
 //! Our implementation follows https://eprint.iacr.org/2022/274.pdf
 //!
 //! This field supports very fast modular arithmetic and has a number of other attractive
@@ -34,16 +35,13 @@ mod tests;
 // CONSTANTS
 // ================================================================================================
 
-// Field modulus = 2^64 - 2^32 + 1
+/// Field modulus = 2^64 - 2^32 + 1
 const M: u64 = 0xFFFFFFFF00000001;
 
 /// 2^128 mod M; this is used for conversion of elements into Montgomery representation.
 const R2: u64 = 0xFFFFFFFE00000001;
 
-// (p+1)/2
-pub const MOD_2: u64 = (0xFFFFFFFF00000001 + 1u64) >> 1;
-
-// 2^32 root of unity
+/// 2^32 root of unity
 const G: u64 = 1753635133440165772;
 
 /// Number of bytes needed to represent field element
@@ -61,32 +59,7 @@ impl BaseElement {
     /// Creates a new field element from the provided `value`; the value is converted into
     /// Montgomery representation.
     pub const fn new(value: u64) -> BaseElement {
-        BaseElement(mont_red_cst((value as u128) * (R2 as u128)))
-    }
-    /// Gets the inner value that might not be canonical
-
-    pub const fn inner(self: &Self) -> u64 {
-        return self.0;
-    }
-
-    /// Multiple squarings in BaseField: return x^(2^n)
-    pub fn msquare(self, n: u32) -> Self {
-        let mut x = self;
-        for _ in 0..n {
-            x = x.square();
-        }
-        x
-    }
-
-    /// Test of equality between two BaseField elements; return value is
-    /// 0xFFFFFFFFFFFFFFFF if the two values are equal, or 0 otherwise.
-    #[inline(always)]
-    pub const fn equals(self, rhs: Self) -> u64 {
-        // Since internal representation is canonical, we can simply
-        // do a xor between the two operands, and then use the same
-        // expression as iszero().
-        let t = self.0 ^ rhs.0;
-        !((((t | t.wrapping_neg()) as i64) >> 63) as u64)
+        Self(mont_red_cst((value as u128) * (R2 as u128)))
     }
 }
 
@@ -256,7 +229,7 @@ impl Display for BaseElement {
 impl PartialEq for BaseElement {
     #[inline]
     fn eq(&self, other: &Self) -> bool {
-        Self::equals(*self, *other) == 0xFFFFFFFFFFFFFFFF
+        equals(self.0, other.0) == 0xFFFFFFFFFFFFFFFF
     }
 }
 
@@ -275,7 +248,7 @@ impl Add for BaseElement {
         // We compute a + b = a - (p - b).
         let (x1, c1) = self.0.overflowing_sub(M - rhs.0);
         let adj = 0u32.wrapping_sub(c1 as u32);
-        BaseElement(x1.wrapping_sub(adj as u64))
+        Self(x1.wrapping_sub(adj as u64))
     }
 }
 
@@ -295,7 +268,7 @@ impl Sub for BaseElement {
         // See reference above for more details.
         let (x1, c1) = self.0.overflowing_sub(rhs.0);
         let adj = 0u32.wrapping_sub(c1 as u32);
-        BaseElement(x1.wrapping_sub(adj as u64))
+        Self(x1.wrapping_sub(adj as u64))
     }
 }
 
@@ -311,7 +284,7 @@ impl Mul for BaseElement {
 
     #[inline]
     fn mul(self, rhs: Self) -> Self {
-        BaseElement(mont_red_cst((self.0 as u128) * (rhs.0 as u128)))
+        Self(mont_red_cst((self.0 as u128) * (rhs.0 as u128)))
     }
 }
 
@@ -344,7 +317,7 @@ impl Neg for BaseElement {
 
     #[inline]
     fn neg(self) -> Self {
-        BaseElement::ZERO - self
+        Self::ZERO - self
     }
 }
 
@@ -425,12 +398,9 @@ impl ExtensibleField<3> for BaseElement {
     fn frobenius(x: [Self; 3]) -> [Self; 3] {
         // coefficients were computed using SageMath
         [
-            x[0] + BaseElement::new(10615703402128488253) * x[1]
-                + BaseElement::new(6700183068485440220) * x[2],
-            BaseElement::new(10050274602728160328) * x[1]
-                + BaseElement::new(14531223735771536287) * x[2],
-            BaseElement::new(11746561000929144102) * x[1]
-                + BaseElement::new(8396469466686423992) * x[2],
+            x[0] + Self::new(10615703402128488253) * x[1] + Self::new(6700183068485440220) * x[2],
+            Self::new(10050274602728160328) * x[1] + Self::new(14531223735771536287) * x[2],
+            Self::new(11746561000929144102) * x[1] + Self::new(8396469466686423992) * x[2],
         ]
     }
 }
@@ -442,8 +412,8 @@ impl From<u128> for BaseElement {
     /// Converts a 128-bit value into a field element.
     fn from(x: u128) -> Self {
         //const R3: u128 = 1 (= 2^192 mod M );// thus we get that mont_red_var((mont_red_var(x) as u128) * R3) becomes
-        //BaseElement(mont_red_var(mont_red_var(x) as u128))  // Variable time implementation
-        BaseElement(mont_red_cst(mont_red_cst(x) as u128)) // Constant time implementation
+        //Self(mont_red_var(mont_red_var(x) as u128))  // Variable time implementation
+        Self(mont_red_cst(mont_red_cst(x) as u128)) // Constant time implementation
     }
 }
 
@@ -451,28 +421,28 @@ impl From<u64> for BaseElement {
     /// Converts a 64-bit value into a field element. If the value is greater than or equal to
     /// the field modulus, modular reduction is silently performed.
     fn from(value: u64) -> Self {
-        BaseElement::new(value)
+        Self::new(value)
     }
 }
 
 impl From<u32> for BaseElement {
     /// Converts a 32-bit value into a field element.
     fn from(value: u32) -> Self {
-        BaseElement::new(value as u64)
+        Self::new(value as u64)
     }
 }
 
 impl From<u16> for BaseElement {
     /// Converts a 16-bit value into a field element.
     fn from(value: u16) -> Self {
-        BaseElement::new(value as u64)
+        Self::new(value as u64)
     }
 }
 
 impl From<u8> for BaseElement {
     /// Converts an 8-bit value into a field element.
     fn from(value: u8) -> Self {
-        BaseElement::new(value as u64)
+        Self::new(value as u64)
     }
 }
 
@@ -483,7 +453,7 @@ impl From<[u8; 8]> for BaseElement {
     /// performed.
     fn from(bytes: [u8; 8]) -> Self {
         let value = u64::from_le_bytes(bytes);
-        BaseElement::new(value)
+        Self::new(value)
     }
 }
 
@@ -518,7 +488,7 @@ impl<'a> TryFrom<&'a [u8]> for BaseElement {
                 value
             )));
         }
-        Ok(BaseElement::new(value))
+        Ok(Self::new(value))
     }
 }
 
@@ -549,7 +519,7 @@ impl Deserializable for BaseElement {
                 value
             )));
         }
-        Ok(BaseElement::new(value))
+        Ok(Self::new(value))
     }
 }
 
@@ -562,9 +532,10 @@ fn exp_acc<const N: usize>(base: BaseElement, tail: BaseElement) -> BaseElement
     }
     result * tail
 }
+
 /// Montgomery reduction (variable time)
 #[inline(always)]
-pub fn mont_red_var(x: u128) -> u64 {
+pub const fn mont_red_var(x: u128) -> u64 {
     const NPRIME: u64 = 4294967297;
     let q = (((x as u64) as u128) * (NPRIME as u128)) as u64;
     let m = (q as u128) * (M as u128);
@@ -575,6 +546,7 @@ pub fn mont_red_var(x: u128) -> u64 {
         return y as u64;
     };
 }
+
 /// Montgomery reduction (constant time)
 #[inline(always)]
 pub const fn mont_red_cst(x: u128) -> u64 {
@@ -588,3 +560,14 @@ pub const fn mont_red_cst(x: u128) -> u64 {
     let (r, c) = xh.overflowing_sub(b);
     r.wrapping_sub(0u32.wrapping_sub(c as u32) as u64)
 }
+
+/// Test of equality between two BaseField elements; return value is
+/// 0xFFFFFFFFFFFFFFFF if the two values are equal, or 0 otherwise.
+#[inline(always)]
+pub fn equals(lhs: u64, rhs: u64) -> u64 {
+    // Since internal representation is canonical, we can simply
+    // do a xor between the two operands, and then use the same
+    // expression as iszero().
+    let t = lhs ^ rhs;
+    !((((t | t.wrapping_neg()) as i64) >> 63) as u64)
+}

From 40e363d35f7327dd63e5f62b19804753978aeb88 Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Wed, 15 Jun 2022 22:22:00 +0200
Subject: [PATCH 5/7] feat: Montgomery rep. based 64bit field

---
 math/src/field/extensions/cubic.rs     | 31 ++++++++++++++++----------
 math/src/field/extensions/quadratic.rs | 31 +++++++++++++-------------
 math/src/field/f64/mod.rs              | 14 +++++++-----
 math/src/field/f64/tests.rs            |  2 +-
 4 files changed, 44 insertions(+), 34 deletions(-)

diff --git a/math/src/field/extensions/cubic.rs b/math/src/field/extensions/cubic.rs
index 3741ee0ae..e3ed7dd41 100644
--- a/math/src/field/extensions/cubic.rs
+++ b/math/src/field/extensions/cubic.rs
@@ -335,7 +335,7 @@ impl<B: ExtensibleField<3>> Deserializable for CubeExtension<B> {
 
 #[cfg(test)]
 mod tests {
-    use super::{CubeExtension, DeserializationError, FieldElement, Vec};
+    use super::{CubeExtension, DeserializationError, FieldElement};
     use crate::field::f64::BaseElement;
     use rand_utils::rand_value;
 
@@ -400,10 +400,13 @@ mod tests {
             ),
         ];
 
-        let expected: Vec<u8> = vec![
-            1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0,
-            0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0,
-        ];
+        let mut expected = vec![];
+        expected.extend_from_slice(&source[0].0.inner().to_le_bytes());
+        expected.extend_from_slice(&source[0].1.inner().to_le_bytes());
+        expected.extend_from_slice(&source[0].2.inner().to_le_bytes());
+        expected.extend_from_slice(&source[1].0.inner().to_le_bytes());
+        expected.extend_from_slice(&source[1].1.inner().to_le_bytes());
+        expected.extend_from_slice(&source[1].2.inner().to_le_bytes());
 
         assert_eq!(
             expected,
@@ -413,12 +416,7 @@ mod tests {
 
     #[test]
     fn bytes_as_elements() {
-        let bytes: Vec<u8> = vec![
-            1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0,
-            0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 7,
-        ];
-
-        let expected = vec![
+        let elements = vec![
             CubeExtension(
                 BaseElement::new(1),
                 BaseElement::new(2),
@@ -431,9 +429,18 @@ mod tests {
             ),
         ];
 
+        let mut bytes = vec![];
+        bytes.extend_from_slice(&elements[0].0.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[0].1.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[0].2.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[1].0.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[1].1.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[1].2.inner().to_le_bytes());
+        bytes.extend_from_slice(&BaseElement::new(5).inner().to_le_bytes());
+
         let result = unsafe { CubeExtension::<BaseElement>::bytes_as_elements(&bytes[..48]) };
         assert!(result.is_ok());
-        assert_eq!(expected, result.unwrap());
+        assert_eq!(elements, result.unwrap());
 
         let result = unsafe { CubeExtension::<BaseElement>::bytes_as_elements(&bytes) };
         assert!(matches!(result, Err(DeserializationError::InvalidValue(_))));
diff --git a/math/src/field/extensions/quadratic.rs b/math/src/field/extensions/quadratic.rs
index 5eee9a09b..d0e7508d2 100644
--- a/math/src/field/extensions/quadratic.rs
+++ b/math/src/field/extensions/quadratic.rs
@@ -324,7 +324,7 @@ impl<B: ExtensibleField<2>> Deserializable for QuadExtension<B> {
 
 #[cfg(test)]
 mod tests {
-    use super::{DeserializationError, FieldElement, QuadExtension, Vec};
+    use super::{DeserializationError, FieldElement, QuadExtension};
     use crate::field::f64::BaseElement;
     use rand_utils::rand_value;
 
@@ -381,11 +381,11 @@ mod tests {
             QuadExtension(BaseElement::new(3), BaseElement::new(4)),
         ];
 
-        let expected: Vec<u8> = vec![
-            1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-            0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-            0, 0, 0, 0, 0, 0,
-        ];
+        let mut expected = vec![];
+        expected.extend_from_slice(&source[0].0.inner().to_le_bytes());
+        expected.extend_from_slice(&source[0].1.inner().to_le_bytes());
+        expected.extend_from_slice(&source[1].0.inner().to_le_bytes());
+        expected.extend_from_slice(&source[1].1.inner().to_le_bytes());
 
         assert_eq!(
             expected,
@@ -395,20 +395,21 @@ mod tests {
 
     #[test]
     fn bytes_as_elements() {
-        let bytes: Vec<u8> = vec![
-            1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-            0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-            0, 0, 0, 0, 0, 0, 5,
-        ];
-
-        let expected = vec![
+        let elements = vec![
             QuadExtension(BaseElement::new(1), BaseElement::new(2)),
             QuadExtension(BaseElement::new(3), BaseElement::new(4)),
         ];
 
-        let result = unsafe { QuadExtension::<BaseElement>::bytes_as_elements(&bytes[..64]) };
+        let mut bytes = vec![];
+        bytes.extend_from_slice(&elements[0].0.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[0].1.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[1].0.inner().to_le_bytes());
+        bytes.extend_from_slice(&elements[1].1.inner().to_le_bytes());
+        bytes.extend_from_slice(&BaseElement::new(5).inner().to_le_bytes());
+        
+        let result = unsafe { QuadExtension::<BaseElement>::bytes_as_elements(&bytes[..32]) };
         assert!(result.is_ok());
-        assert_eq!(expected, result.unwrap());
+        assert_eq!(elements, result.unwrap());
 
         let result = unsafe { QuadExtension::<BaseElement>::bytes_as_elements(&bytes) };
         assert!(matches!(result, Err(DeserializationError::InvalidValue(_))));
diff --git a/math/src/field/f64/mod.rs b/math/src/field/f64/mod.rs
index 9d92f4d8d..7686f953f 100644
--- a/math/src/field/f64/mod.rs
+++ b/math/src/field/f64/mod.rs
@@ -5,7 +5,8 @@
 
 //! An implementation of a 64-bit STARK-friendly prime field with modulus $2^{64} - 2^{32} + 1$
 //! using Montgomery representation.
-//! Our implementation follows https://eprint.iacr.org/2022/274.pdf
+//! Our implementation follows https://eprint.iacr.org/2022/274.pdf and is constant-time except
+//! for the mont_red_var() function.
 //!
 //! This field supports very fast modular arithmetic and has a number of other attractive
 //! properties, including:
@@ -61,6 +62,11 @@ impl BaseElement {
     pub const fn new(value: u64) -> BaseElement {
         Self(mont_red_cst((value as u128) * (R2 as u128)))
     }
+
+    /// Returns the non-canonical u64 inner value. 
+    pub const fn inner(&self) -> u64{
+        self.0
+    }
 }
 
 impl FieldElement for BaseElement {
@@ -265,7 +271,6 @@ impl Sub for BaseElement {
     #[inline]
     #[allow(clippy::suspicious_arithmetic_impl)]
     fn sub(self, rhs: Self) -> Self {
-        // See reference above for more details.
         let (x1, c1) = self.0.overflowing_sub(rhs.0);
         let adj = 0u32.wrapping_sub(c1 as u32);
         Self(x1.wrapping_sub(adj as u64))
@@ -301,7 +306,7 @@ impl Div for BaseElement {
     #[inline]
     #[allow(clippy::suspicious_arithmetic_impl)]
     fn div(self, rhs: Self) -> Self {
-        self * Self::inv(rhs)
+        self * rhs.inv()
     }
 }
 
@@ -565,9 +570,6 @@ pub const fn mont_red_cst(x: u128) -> u64 {
 /// 0xFFFFFFFFFFFFFFFF if the two values are equal, or 0 otherwise.
 #[inline(always)]
 pub fn equals(lhs: u64, rhs: u64) -> u64 {
-    // Since internal representation is canonical, we can simply
-    // do a xor between the two operands, and then use the same
-    // expression as iszero().
     let t = lhs ^ rhs;
     !((((t | t.wrapping_neg()) as i64) >> 63) as u64)
 }
diff --git a/math/src/field/f64/tests.rs b/math/src/field/f64/tests.rs
index 033cd2002..5a3f5259d 100644
--- a/math/src/field/f64/tests.rs
+++ b/math/src/field/f64/tests.rs
@@ -5,7 +5,7 @@
 
 use super::{
     //E,
-    AsBytes, BaseElement, DeserializationError, FieldElement, Serializable, StarkField, M,
+    BaseElement, DeserializationError, FieldElement, Serializable, StarkField, M,
 };
 use crate::field::{CubeExtension, ExtensionOf, QuadExtension};
 use core::convert::TryFrom;

From 67491e0ebf3bc6ebaed6ee948af8d2842e683d48 Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Fri, 17 Jun 2022 01:01:24 +0200
Subject: [PATCH 6/7] chore: fix typos

---
 math/src/field/f64/tests.rs | 12 +-----------
 1 file changed, 1 insertion(+), 11 deletions(-)

diff --git a/math/src/field/f64/tests.rs b/math/src/field/f64/tests.rs
index 5a3f5259d..23a743262 100644
--- a/math/src/field/f64/tests.rs
+++ b/math/src/field/f64/tests.rs
@@ -4,7 +4,6 @@
 // LICENSE file in the root directory of this source tree.
 
 use super::{
-    //E,
     BaseElement, DeserializationError, FieldElement, Serializable, StarkField, M,
 };
 use crate::field::{CubeExtension, ExtensionOf, QuadExtension};
@@ -32,11 +31,6 @@ fn add() {
     let t = BaseElement::new(M - 1);
     assert_eq!(BaseElement::ZERO, t + BaseElement::ONE);
     assert_eq!(BaseElement::ONE, t + BaseElement::new(2));
-
-    //// test non-canonical representation
-    //let a = BaseElement::new(M - 1) + BaseElement::new(E);
-    //let expected = ((((M - 1 + E) as u128) * 2) % (M as u128)) as u64;
-    //assert_eq!(expected, (a + a).as_int());
 }
 
 #[test]
@@ -72,7 +66,7 @@ fn mul() {
     assert_eq!(BaseElement::ZERO, r * BaseElement::ZERO);
     assert_eq!(r, r * BaseElement::ONE);
 
-    // test multifield::extensions::cubic::tests::bytes_as_elementsplication within bounds
+    // test multiplication within bounds
     assert_eq!(
         BaseElement::from(15u8),
         BaseElement::from(5u8) * BaseElement::from(3u8)
@@ -131,10 +125,6 @@ fn equals() {
     assert_eq!(a, b);
     assert_eq!(a.as_int(), b.as_int());
     assert_eq!(a.to_bytes(), b.to_bytes());
-
-    // but their internal representation is not
-    //assert_ne!(a.0, b.0);
-    //assert_ne!(a.as_bytes(), b.as_bytes());
 }
 
 // ROOTS OF UNITY

From f1b7cef217e10101e6f4730cb72dc9d8f2706c89 Mon Sep 17 00:00:00 2001
From: Al-Kindi-0 <oamine0101@gmail.com>
Date: Tue, 28 Jun 2022 22:36:45 +0200
Subject: [PATCH 7/7] feat: Optimized Rescue prime in Montgomery

---
 crypto/src/hash/rescue/rp64_256/mod.rs | 538 ++++++++-----------------
 math/src/field/f64/mod.rs              |   6 +
 2 files changed, 180 insertions(+), 364 deletions(-)

diff --git a/crypto/src/hash/rescue/rp64_256/mod.rs b/crypto/src/hash/rescue/rp64_256/mod.rs
index d6e7f8e96..5e2c2480b 100644
--- a/crypto/src/hash/rescue/rp64_256/mod.rs
+++ b/crypto/src/hash/rescue/rp64_256/mod.rs
@@ -268,12 +268,6 @@ impl Rp64_256 {
     /// The output of the hash function can be read from state elements 4, 5, 6, and 7.
     pub const DIGEST_RANGE: Range<usize> = DIGEST_RANGE;
 
-    /// MDS matrix used for computing the linear layer in a Rescue Prime round.
-    pub const MDS: [[BaseElement; STATE_WIDTH]; STATE_WIDTH] = MDS;
-
-    /// Inverse of the MDS matrix.
-    pub const INV_MDS: [[BaseElement; STATE_WIDTH]; STATE_WIDTH] = INV_MDS;
-
     /// Round constants added to the hasher state in the first half of the Rescue Prime round.
     pub const ARK1: [[BaseElement; STATE_WIDTH]; NUM_ROUNDS] = ARK1;
 
@@ -310,14 +304,34 @@ impl Rp64_256 {
     // --------------------------------------------------------------------------------------------
 
     #[inline(always)]
-    fn apply_mds(state: &mut [BaseElement; STATE_WIDTH]) {
+    fn apply_mds(state_: &mut [BaseElement; STATE_WIDTH]) {
         let mut result = [BaseElement::ZERO; STATE_WIDTH];
-        result.iter_mut().zip(MDS).for_each(|(r, mds_row)| {
-            state.iter().zip(mds_row).for_each(|(&s, m)| {
-                *r += m * s;
-            });
-        });
-        *state = result
+
+        // Using the linearity of the operations we can split the state into a low||high decomposition
+        // and operate on each with no overflow and then combine/reduce the result to a field element.
+        let mut state_l = [0u64; STATE_WIDTH];
+        let mut state_h = [0u64; STATE_WIDTH];
+
+        for r in 0..STATE_WIDTH {
+            let s = state_[r].inner();
+            state_h[r] = s >> 32;
+            state_l[r] = (s as u32) as u64;
+        }
+
+        let state_h = Self::mds_multiply_freq(state_h);
+        let state_l = Self::mds_multiply_freq(state_l);
+
+        for r in 0..STATE_WIDTH {
+            let s = state_l[r] as u128 + ((state_h[r] as u128) << 32);
+            let s_hi = (s >> 64) as u64;
+            let s_lo = s as u64;
+            let z = (s_hi << 32) - s_hi;
+            let (res, over) = s_lo.overflowing_add(z);
+
+            result[r] =
+                BaseElement::new_unsafe(res.wrapping_add(0u32.wrapping_sub(over as u32) as u64));
+        }
+        *state_ = result;
     }
 
     #[inline(always)]
@@ -327,11 +341,26 @@ impl Rp64_256 {
 
     #[inline(always)]
     fn apply_sbox(state: &mut [BaseElement; STATE_WIDTH]) {
-        state.iter_mut().for_each(|v| {
-            let t2 = v.square();
-            let t4 = t2.square();
-            *v *= t2 * t4;
-        });
+        state[0] = Self::pow_7(state[0]);
+        state[1] = Self::pow_7(state[1]);
+        state[2] = Self::pow_7(state[2]);
+        state[3] = Self::pow_7(state[3]);
+        state[4] = Self::pow_7(state[4]);
+        state[5] = Self::pow_7(state[5]);
+        state[6] = Self::pow_7(state[6]);
+        state[7] = Self::pow_7(state[7]);
+        state[8] = Self::pow_7(state[8]);
+        state[9] = Self::pow_7(state[9]);
+        state[10] = Self::pow_7(state[10]);
+        state[11] = Self::pow_7(state[11]);
+    }
+
+    #[inline(always)]
+    fn pow_7(x: BaseElement) -> BaseElement {
+        let x2 = x.square();
+        let x4 = x2.square();
+        let x3 = x * x2;
+        x3 * x4
     }
 
     #[inline(always)]
@@ -369,356 +398,137 @@ impl Rp64_256 {
             *s = a * b;
         }
     }
+
+    // FFT MDS MULTIPLICATION HELPER FUNCTIONS
+    // --------------------------------------------------------------------------------------------
+
+    // We use split 3 x 4 FFT transform in order to transform our vectors into the frequency domain.
+    #[inline(always)]
+    fn mds_multiply_freq(state: [u64; 12]) -> [u64; 12] {
+        let [s0, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11] = state;
+
+        let (u0, u1, u2) = Self::fft4_real([s0, s3, s6, s9]);
+        let (u4, u5, u6) = Self::fft4_real([s1, s4, s7, s10]);
+        let (u8, u9, u10) = Self::fft4_real([s2, s5, s8, s11]);
+
+        // The 4th block is not computed as it is similar to the 2nd one, up to complex conjugation,
+        // and due to the use of the real FFT and iFFT is redundant.
+        let [v0, v4, v8] = Self::block1([u0, u4, u8], MDS_FREQ_BLOCK_ONE);
+        let [v1, v5, v9] = Self::block2([u1, u5, u9], MDS_FREQ_BLOCK_TWO);
+        let [v2, v6, v10] = Self::block3([u2, u6, u10], MDS_FREQ_BLOCK_THREE);
+
+        let [s0, s3, s6, s9] = Self::ifft4_real((v0, v1, v2));
+        let [s1, s4, s7, s10] = Self::ifft4_real((v4, v5, v6));
+        let [s2, s5, s8, s11] = Self::ifft4_real((v8, v9, v10));
+
+        return [s0, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11];
+    }
+
+    // We use the real FFT to avoid redundant computations. See https://www.mdpi.com/2076-3417/12/9/4700
+    #[inline(always)]
+    fn fft2_real(x: [u64; 2]) -> [i64; 2] {
+        return [(x[0] as i64 + x[1] as i64), (x[0] as i64 - x[1] as i64)];
+    }
+
+    #[inline(always)]
+    fn ifft2_real(y: [i64; 2]) -> [u64; 2] {
+        return [(y[0] + y[1]) as u64 >> 1, (y[0] - y[1]) as u64 >> 1];
+    }
+
+    #[inline(always)]
+    fn fft4_real(x: [u64; 4]) -> (i64, (i64, i64), i64) {
+        let [z0, z2] = Self::fft2_real([x[0], x[2]]);
+        let [z1, z3] = Self::fft2_real([x[1], x[3]]);
+        let y0 = z0 + z1;
+        let y1 = (z2, -z3);
+        let y2 = z0 - z1;
+        return (y0, y1, y2);
+    }
+
+    #[inline(always)]
+    fn ifft4_real(y: (i64, (i64, i64), i64)) -> [u64; 4] {
+        let z0 = (y.0 + y.2) >> 1;
+        let z1 = (y.0 - y.2) >> 1;
+        let z2 = y.1 .0;
+        let z3 = -y.1 .1;
+
+        let [x0, x2] = Self::ifft2_real([z0, z2]);
+        let [x1, x3] = Self::ifft2_real([z1, z3]);
+
+        return [x0, x1, x2, x3];
+    }
+
+    #[inline(always)]
+    fn block1(x: [i64; 3], y: [i64; 3]) -> [i64; 3] {
+        let [x0, x1, x2] = x;
+        let [y0, y1, y2] = y;
+        let z0 = x0 * y0 + x1 * y2 + x2 * y1;
+        let z1 = x0 * y1 + x1 * y0 + x2 * y2;
+        let z2 = x0 * y2 + x1 * y1 + x2 * y0;
+
+        return [z0, z1, z2];
+    }
+
+    #[inline(always)]
+    fn block2(x: [(i64, i64); 3], y: [(i64, i64); 3]) -> [(i64, i64); 3] {
+        let [(x0r, x0i), (x1r, x1i), (x2r, x2i)] = x;
+        let [(y0r, y0i), (y1r, y1i), (y2r, y2i)] = y;
+        let x0s = x0r + x0i;
+        let x1s = x1r + x1i;
+        let x2s = x2r + x2i;
+        let y0s = y0r + y0i;
+        let y1s = y1r + y1i;
+        let y2s = y2r + y2i;
+
+        // Compute x0​y0 ​− ix1​y2​ − ix2​y1​ using Karatsuba for complex numbers multiplication
+        let m0 = (x0r * y0r, x0i * y0i);
+        let m1 = (x1r * y2r, x1i * y2i);
+        let m2 = (x2r * y1r, x2i * y1i);
+        let z0r = (m0.0 - m0.1) + (x1s * y2s - m1.0 - m1.1) + (x2s * y1s - m2.0 - m2.1);
+        let z0i = (x0s * y0s - m0.0 - m0.1) + (-m1.0 + m1.1) + (-m2.0 + m2.1);
+        let z0 = (z0r, z0i);
+
+        // Compute x0​y1​ + x1​y0​ − ix2​y2 using Karatsuba for complex numbers multiplication
+        let m0 = (x0r * y1r, x0i * y1i);
+        let m1 = (x1r * y0r, x1i * y0i);
+        let m2 = (x2r * y2r, x2i * y2i);
+        let z1r = (m0.0 - m0.1) + (m1.0 - m1.1) + (x2s * y2s - m2.0 - m2.1);
+        let z1i = (x0s * y1s - m0.0 - m0.1) + (x1s * y0s - m1.0 - m1.1) + (-m2.0 + m2.1);
+        let z1 = (z1r, z1i);
+
+        // Compute x0​y2​ + x1​y1 ​+ x2​y0​ using Karatsuba for complex numbers multiplication
+        let m0 = (x0r * y2r, x0i * y2i);
+        let m1 = (x1r * y1r, x1i * y1i);
+        let m2 = (x2r * y0r, x2i * y0i);
+        let z2r = (m0.0 - m0.1) + (m1.0 - m1.1) + (m2.0 - m2.1);
+        let z2i = (x0s * y2s - m0.0 - m0.1) + (x1s * y1s - m1.0 - m1.1) + (x2s * y0s - m2.0 - m2.1);
+        let z2 = (z2r, z2i);
+
+        return [z0, z1, z2];
+    }
+
+    #[inline(always)]
+    fn block3(x: [i64; 3], y: [i64; 3]) -> [i64; 3] {
+        let [x0, x1, x2] = x;
+        let [y0, y1, y2] = y;
+        let z0 = x0 * y0 - x1 * y2 - x2 * y1;
+        let z1 = x0 * y1 + x1 * y0 - x2 * y2;
+        let z2 = x0 * y2 + x1 * y1 + x2 * y0;
+
+        return [z0, z1, z2];
+    }
 }
 
 // MDS
 // ================================================================================================
-/// Rescue MDS matrix
-/// Computed using algorithm 4 from <https://eprint.iacr.org/2020/1143.pdf>
-const MDS: [[BaseElement; STATE_WIDTH]; STATE_WIDTH] = [
-    [
-        BaseElement::new(2108866337646019936),
-        BaseElement::new(11223275256334781131),
-        BaseElement::new(2318414738826783588),
-        BaseElement::new(11240468238955543594),
-        BaseElement::new(8007389560317667115),
-        BaseElement::new(11080831380224887131),
-        BaseElement::new(3922954383102346493),
-        BaseElement::new(17194066286743901609),
-        BaseElement::new(152620255842323114),
-        BaseElement::new(7203302445933022224),
-        BaseElement::new(17781531460838764471),
-        BaseElement::new(2306881200),
-    ],
-    [
-        BaseElement::new(3368836954250922620),
-        BaseElement::new(5531382716338105518),
-        BaseElement::new(7747104620279034727),
-        BaseElement::new(14164487169476525880),
-        BaseElement::new(4653455932372793639),
-        BaseElement::new(5504123103633670518),
-        BaseElement::new(3376629427948045767),
-        BaseElement::new(1687083899297674997),
-        BaseElement::new(8324288417826065247),
-        BaseElement::new(17651364087632826504),
-        BaseElement::new(15568475755679636039),
-        BaseElement::new(4656488262337620150),
-    ],
-    [
-        BaseElement::new(2560535215714666606),
-        BaseElement::new(10793518538122219186),
-        BaseElement::new(408467828146985886),
-        BaseElement::new(13894393744319723897),
-        BaseElement::new(17856013635663093677),
-        BaseElement::new(14510101432365346218),
-        BaseElement::new(12175743201430386993),
-        BaseElement::new(12012700097100374591),
-        BaseElement::new(976880602086740182),
-        BaseElement::new(3187015135043748111),
-        BaseElement::new(4630899319883688283),
-        BaseElement::new(17674195666610532297),
-    ],
-    [
-        BaseElement::new(10940635879119829731),
-        BaseElement::new(9126204055164541072),
-        BaseElement::new(13441880452578323624),
-        BaseElement::new(13828699194559433302),
-        BaseElement::new(6245685172712904082),
-        BaseElement::new(3117562785727957263),
-        BaseElement::new(17389107632996288753),
-        BaseElement::new(3643151412418457029),
-        BaseElement::new(10484080975961167028),
-        BaseElement::new(4066673631745731889),
-        BaseElement::new(8847974898748751041),
-        BaseElement::new(9548808324754121113),
-    ],
-    [
-        BaseElement::new(15656099696515372126),
-        BaseElement::new(309741777966979967),
-        BaseElement::new(16075523529922094036),
-        BaseElement::new(5384192144218250710),
-        BaseElement::new(15171244241641106028),
-        BaseElement::new(6660319859038124593),
-        BaseElement::new(6595450094003204814),
-        BaseElement::new(15330207556174961057),
-        BaseElement::new(2687301105226976975),
-        BaseElement::new(15907414358067140389),
-        BaseElement::new(2767130804164179683),
-        BaseElement::new(8135839249549115549),
-    ],
-    [
-        BaseElement::new(14687393836444508153),
-        BaseElement::new(8122848807512458890),
-        BaseElement::new(16998154830503301252),
-        BaseElement::new(2904046703764323264),
-        BaseElement::new(11170142989407566484),
-        BaseElement::new(5448553946207765015),
-        BaseElement::new(9766047029091333225),
-        BaseElement::new(3852354853341479440),
-        BaseElement::new(14577128274897891003),
-        BaseElement::new(11994931371916133447),
-        BaseElement::new(8299269445020599466),
-        BaseElement::new(2859592328380146288),
-    ],
-    [
-        BaseElement::new(4920761474064525703),
-        BaseElement::new(13379538658122003618),
-        BaseElement::new(3169184545474588182),
-        BaseElement::new(15753261541491539618),
-        BaseElement::new(622292315133191494),
-        BaseElement::new(14052907820095169428),
-        BaseElement::new(5159844729950547044),
-        BaseElement::new(17439978194716087321),
-        BaseElement::new(9945483003842285313),
-        BaseElement::new(13647273880020281344),
-        BaseElement::new(14750994260825376),
-        BaseElement::new(12575187259316461486),
-    ],
-    [
-        BaseElement::new(3371852905554824605),
-        BaseElement::new(8886257005679683950),
-        BaseElement::new(15677115160380392279),
-        BaseElement::new(13242906482047961505),
-        BaseElement::new(12149996307978507817),
-        BaseElement::new(1427861135554592284),
-        BaseElement::new(4033726302273030373),
-        BaseElement::new(14761176804905342155),
-        BaseElement::new(11465247508084706095),
-        BaseElement::new(12112647677590318112),
-        BaseElement::new(17343938135425110721),
-        BaseElement::new(14654483060427620352),
-    ],
-    [
-        BaseElement::new(5421794552262605237),
-        BaseElement::new(14201164512563303484),
-        BaseElement::new(5290621264363227639),
-        BaseElement::new(1020180205893205576),
-        BaseElement::new(14311345105258400438),
-        BaseElement::new(7828111500457301560),
-        BaseElement::new(9436759291445548340),
-        BaseElement::new(5716067521736967068),
-        BaseElement::new(15357555109169671716),
-        BaseElement::new(4131452666376493252),
-        BaseElement::new(16785275933585465720),
-        BaseElement::new(11180136753375315897),
-    ],
-    [
-        BaseElement::new(10451661389735482801),
-        BaseElement::new(12128852772276583847),
-        BaseElement::new(10630876800354432923),
-        BaseElement::new(6884824371838330777),
-        BaseElement::new(16413552665026570512),
-        BaseElement::new(13637837753341196082),
-        BaseElement::new(2558124068257217718),
-        BaseElement::new(4327919242598628564),
-        BaseElement::new(4236040195908057312),
-        BaseElement::new(2081029262044280559),
-        BaseElement::new(2047510589162918469),
-        BaseElement::new(6835491236529222042),
-    ],
-    [
-        BaseElement::new(5675273097893923172),
-        BaseElement::new(8120839782755215647),
-        BaseElement::new(9856415804450870143),
-        BaseElement::new(1960632704307471239),
-        BaseElement::new(15279057263127523057),
-        BaseElement::new(17999325337309257121),
-        BaseElement::new(72970456904683065),
-        BaseElement::new(8899624805082057509),
-        BaseElement::new(16980481565524365258),
-        BaseElement::new(6412696708929498357),
-        BaseElement::new(13917768671775544479),
-        BaseElement::new(5505378218427096880),
-    ],
-    [
-        BaseElement::new(10318314766641004576),
-        BaseElement::new(17320192463105632563),
-        BaseElement::new(11540812969169097044),
-        BaseElement::new(7270556942018024148),
-        BaseElement::new(4755326086930560682),
-        BaseElement::new(2193604418377108959),
-        BaseElement::new(11681945506511803967),
-        BaseElement::new(8000243866012209465),
-        BaseElement::new(6746478642521594042),
-        BaseElement::new(12096331252283646217),
-        BaseElement::new(13208137848575217268),
-        BaseElement::new(5548519654341606996),
-    ],
-];
-
-/// Rescue Inverse MDS matrix
-/// Computed using algorithm 4 from <https://eprint.iacr.org/2020/1143.pdf> and then
-/// inverting the resulting matrix.
-const INV_MDS: [[BaseElement; STATE_WIDTH]; STATE_WIDTH] = [
-    [
-        BaseElement::new(1025714968950054217),
-        BaseElement::new(2820417286206414279),
-        BaseElement::new(4993698564949207576),
-        BaseElement::new(12970218763715480197),
-        BaseElement::new(15096702659601816313),
-        BaseElement::new(5737881372597660297),
-        BaseElement::new(13327263231927089804),
-        BaseElement::new(4564252978131632277),
-        BaseElement::new(16119054824480892382),
-        BaseElement::new(6613927186172915989),
-        BaseElement::new(6454498710731601655),
-        BaseElement::new(2510089799608156620),
-    ],
-    [
-        BaseElement::new(14311337779007263575),
-        BaseElement::new(10306799626523962951),
-        BaseElement::new(7776331823117795156),
-        BaseElement::new(4922212921326569206),
-        BaseElement::new(8669179866856828412),
-        BaseElement::new(936244772485171410),
-        BaseElement::new(4077406078785759791),
-        BaseElement::new(2938383611938168107),
-        BaseElement::new(16650590241171797614),
-        BaseElement::new(16578411244849432284),
-        BaseElement::new(17600191004694808340),
-        BaseElement::new(5913375445729949081),
-    ],
-    [
-        BaseElement::new(13640353831792923980),
-        BaseElement::new(1583879644687006251),
-        BaseElement::new(17678309436940389401),
-        BaseElement::new(6793918274289159258),
-        BaseElement::new(3594897835134355282),
-        BaseElement::new(2158539885379341689),
-        BaseElement::new(12473871986506720374),
-        BaseElement::new(14874332242561185932),
-        BaseElement::new(16402478875851979683),
-        BaseElement::new(9893468322166516227),
-        BaseElement::new(8142413325661539529),
-        BaseElement::new(3444000755516388321),
-    ],
-    [
-        BaseElement::new(14009777257506018221),
-        BaseElement::new(18218829733847178457),
-        BaseElement::new(11151899210182873569),
-        BaseElement::new(14653120475631972171),
-        BaseElement::new(9591156713922565586),
-        BaseElement::new(16622517275046324812),
-        BaseElement::new(3958136700677573712),
-        BaseElement::new(2193274161734965529),
-        BaseElement::new(15125079516929063010),
-        BaseElement::new(3648852869044193741),
-        BaseElement::new(4405494440143722315),
-        BaseElement::new(15549070131235639125),
-    ],
-    [
-        BaseElement::new(14324333194410783741),
-        BaseElement::new(12565645879378458115),
-        BaseElement::new(4028590290335558535),
-        BaseElement::new(17936155181893467294),
-        BaseElement::new(1833939650657097992),
-        BaseElement::new(14310984655970610026),
-        BaseElement::new(4701042357351086687),
-        BaseElement::new(1226379890265418475),
-        BaseElement::new(2550212856624409740),
-        BaseElement::new(5670703442709406167),
-        BaseElement::new(3281485106506301394),
-        BaseElement::new(9804247840970323440),
-    ],
-    [
-        BaseElement::new(7778523590474814059),
-        BaseElement::new(7154630063229321501),
-        BaseElement::new(17790326505487126055),
-        BaseElement::new(3160574440608126866),
-        BaseElement::new(7292349907185131376),
-        BaseElement::new(1916491575080831825),
-        BaseElement::new(11523142515674812675),
-        BaseElement::new(2162357063341827157),
-        BaseElement::new(6650415936886875699),
-        BaseElement::new(11522955632464608509),
-        BaseElement::new(16740856792338897018),
-        BaseElement::new(16987840393715133187),
-    ],
-    [
-        BaseElement::new(14499296811525152023),
-        BaseElement::new(118549270069446537),
-        BaseElement::new(3041471724857448013),
-        BaseElement::new(3827228106225598612),
-        BaseElement::new(2081369067662751050),
-        BaseElement::new(15406142490454329462),
-        BaseElement::new(8943531526276617760),
-        BaseElement::new(3545513411057560337),
-        BaseElement::new(11433277564645295966),
-        BaseElement::new(9558995950666358829),
-        BaseElement::new(7443251815414752292),
-        BaseElement::new(12335092608217610725),
-    ],
-    [
-        BaseElement::new(184304165023253232),
-        BaseElement::new(11596940249585433199),
-        BaseElement::new(18170668175083122019),
-        BaseElement::new(8318891703682569182),
-        BaseElement::new(4387895409295967519),
-        BaseElement::new(14599228871586336059),
-        BaseElement::new(2861651216488619239),
-        BaseElement::new(567601091253927304),
-        BaseElement::new(10135289435539766316),
-        BaseElement::new(14905738261734377063),
-        BaseElement::new(3345637344934149303),
-        BaseElement::new(3159874422865401171),
-    ],
-    [
-        BaseElement::new(1134458872778032479),
-        BaseElement::new(4102035717681749376),
-        BaseElement::new(14030271225872148070),
-        BaseElement::new(10312336662487337312),
-        BaseElement::new(12938229830489392977),
-        BaseElement::new(17758804398255988457),
-        BaseElement::new(15482323580054918356),
-        BaseElement::new(1010277923244261213),
-        BaseElement::new(12904552397519353856),
-        BaseElement::new(5073478003078459047),
-        BaseElement::new(11514678194579805863),
-        BaseElement::new(4419017610446058921),
-    ],
-    [
-        BaseElement::new(2916054498252226520),
-        BaseElement::new(9880379926449218161),
-        BaseElement::new(15314650755395914465),
-        BaseElement::new(8335514387550394159),
-        BaseElement::new(8955267746483690029),
-        BaseElement::new(16353914237438359160),
-        BaseElement::new(4173425891602463552),
-        BaseElement::new(14892581052359168234),
-        BaseElement::new(17561678290843148035),
-        BaseElement::new(7292975356887551984),
-        BaseElement::new(18039512759118984712),
-        BaseElement::new(5411253583520971237),
-    ],
-    [
-        BaseElement::new(9848042270158364544),
-        BaseElement::new(809689769037458603),
-        BaseElement::new(5884047526712050760),
-        BaseElement::new(12956871945669043745),
-        BaseElement::new(14265127496637532237),
-        BaseElement::new(6211568220597222123),
-        BaseElement::new(678544061771515015),
-        BaseElement::new(16295989318674734123),
-        BaseElement::new(11782767968925152203),
-        BaseElement::new(1359397660819991739),
-        BaseElement::new(16148400912425385689),
-        BaseElement::new(14440017265059055146),
-    ],
-    [
-        BaseElement::new(1634272668217219807),
-        BaseElement::new(16290589064070324125),
-        BaseElement::new(5311838222680798126),
-        BaseElement::new(15044064140936894715),
-        BaseElement::new(15775025788428030421),
-        BaseElement::new(12586374713559327349),
-        BaseElement::new(8118943473454062014),
-        BaseElement::new(13223746794660766349),
-        BaseElement::new(13059674280609257192),
-        BaseElement::new(16605443174349648289),
-        BaseElement::new(13586971219878687822),
-        BaseElement::new(16337009014471658360),
-    ],
-];
+/// Rescue MDS matrix in frequency domain. More precisely, this is the output of the 3 4-point 
+/// (real) FFT i.e. just before the multiplication with the appropriate twiddle factors and 
+/// application of the final 4 3-point FFT to get the represention of the MDS matrix, which is
+/// the circular matrix with first row [17, 15, 41, 16, 2, 28, 13, 13, 39, 18, 34, 20],
+/// in frequency domain.
+const MDS_FREQ_BLOCK_ONE: [i64; 3] = [64, 128, 64];
+const MDS_FREQ_BLOCK_TWO: [(i64, i64); 3] = [(4, -2), (-8, 2), (32, 2)];
+const MDS_FREQ_BLOCK_THREE: [i64; 3] = [-4, -32, 8];
 
 // ROUND CONSTANTS
 // ================================================================================================
diff --git a/math/src/field/f64/mod.rs b/math/src/field/f64/mod.rs
index 7686f953f..c39dd5908 100644
--- a/math/src/field/f64/mod.rs
+++ b/math/src/field/f64/mod.rs
@@ -67,6 +67,12 @@ impl BaseElement {
     pub const fn inner(&self) -> u64{
         self.0
     }
+
+    /// Returns a new field element from the provided 'value'. Assumes that 'value' is already
+    /// in canonical Montgomery form.
+    pub const fn new_unsafe(value: u64) -> BaseElement {
+        BaseElement(value)
+    }
 }
 
 impl FieldElement for BaseElement {