From 2d8892e2d7b0e1e586ad49f16202521e46278cc4 Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Thu, 10 Apr 2025 17:11:16 +0200
Subject: [PATCH 01/10] feat(fft): add fast multiplication

---
 crates/math/src/fft/polynomial.rs             | 17 ++++++
 .../src/field/fields/mersenne31/extensions.rs | 52 ++++++++++++++++++-
 2 files changed, 68 insertions(+), 1 deletion(-)

diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index 143e0addc..2bb4fcfef 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -97,6 +97,18 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         let scaled = Polynomial::interpolate_fft::<F>(fft_evals)?;
         Ok(scaled.scale(&offset.inv().unwrap()))
     }
+
+    pub fn fast_multiplication<F: IsFFTField + IsSubFieldOf<E>>(
+        &self,
+        other: &Self,
+    ) -> Result<Self, FFTError> {
+        let domain_size = self.degree() + other.degree() + 1;
+        let p = Polynomial::evaluate_fft::<F>(self, 1, Some(domain_size))?;
+        let q = Polynomial::evaluate_fft::<F>(other, 1, Some(domain_size))?;
+        let r = p.into_iter().zip(q).map(|(a, b)| a * b).collect::<Vec<_>>();
+
+        Polynomial::interpolate_fft::<F>(&r)
+    }
 }
 
 pub fn compose_fft<F, E>(
@@ -313,6 +325,11 @@ mod tests {
 
                 prop_assert_eq!(poly, new_poly);
             }
+
+            #[test]
+            fn test_fft_multiplication_works(poly in poly(7), other in poly(7)) {
+                prop_assert_eq!(poly.fast_multiplication::<F>(&other).unwrap(), poly * other);
+            }
         }
 
         #[test]
diff --git a/crates/math/src/field/fields/mersenne31/extensions.rs b/crates/math/src/field/fields/mersenne31/extensions.rs
index 081548901..37773c165 100644
--- a/crates/math/src/field/fields/mersenne31/extensions.rs
+++ b/crates/math/src/field/fields/mersenne31/extensions.rs
@@ -2,7 +2,7 @@ use super::field::Mersenne31Field;
 use crate::field::{
     element::FieldElement,
     errors::FieldError,
-    traits::{IsField, IsSubFieldOf},
+    traits::{IsFFTField, IsField, IsSubFieldOf},
 };
 #[cfg(feature = "alloc")]
 use alloc::vec::Vec;
@@ -93,6 +93,56 @@ impl IsField for Degree2ExtensionField {
     }
 }
 
+impl IsFFTField for Degree2ExtensionField {
+    const TWO_ADICITY: u64 = 31;
+    const TWO_ADIC_PRIMITVE_ROOT_OF_UNITY: Self::BaseType =
+        [FpE::const_from_raw(2), FpE::const_from_raw(1268011823)];
+}
+
+impl IsSubFieldOf<Degree4ExtensionField> for Degree2ExtensionField {
+    fn add(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        [Fp2E::from(a) + &b[0], b[1].clone()]
+    }
+
+    fn sub(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        [Fp2E::from(a) - &b[0], -&b[1]]
+    }
+
+    fn mul(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> <Degree4ExtensionField as IsField>::BaseType {
+        [Fp2E::from(a) * &b[0], Fp2E::from(a) * &b[1]]
+    }
+
+    fn div(
+        a: &Self::BaseType,
+        b: &<Degree4ExtensionField as IsField>::BaseType,
+    ) -> Result<<Degree4ExtensionField as IsField>::BaseType, FieldError> {
+        let b_inv = Degree4ExtensionField::inv(b).map_err(|_| FieldError::DivisionByZero)?;
+        Ok(<Self as IsSubFieldOf<Degree4ExtensionField>>::mul(
+            a, &b_inv,
+        ))
+    }
+
+    fn embed(a: Self::BaseType) -> <Degree4ExtensionField as IsField>::BaseType {
+        [FieldElement::from_raw(a), FieldElement::zero()]
+    }
+
+    #[cfg(feature = "alloc")]
+    fn to_subfield_vec(
+        b: <Degree4ExtensionField as IsField>::BaseType,
+    ) -> alloc::vec::Vec<Self::BaseType> {
+        b.into_iter().map(|x| x.to_raw()).collect()
+    }
+}
+
 impl IsSubFieldOf<Degree2ExtensionField> for Mersenne31Field {
     fn add(
         a: &Self::BaseType,

From efb389e23bd1cdf72566103490b842a7293b9a0e Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Thu, 10 Apr 2025 18:18:18 +0200
Subject: [PATCH 02/10] chore: add fast_mul benchmark

---
 crates/math/benches/polynomials/polynomial.rs | 25 ++++++++++-
 crates/math/benches/polynomials/utils.rs      | 30 ++++++++++++-
 .../src/field/fields/mersenne31/extensions.rs | 44 -------------------
 3 files changed, 52 insertions(+), 47 deletions(-)

diff --git a/crates/math/benches/polynomials/polynomial.rs b/crates/math/benches/polynomials/polynomial.rs
index ee63f4b74..038cde97b 100644
--- a/crates/math/benches/polynomials/polynomial.rs
+++ b/crates/math/benches/polynomials/polynomial.rs
@@ -1,8 +1,10 @@
-use super::utils::{rand_field_elements, rand_poly, FE};
+use super::utils::{rand_complex_mersenne_poly, rand_field_elements, rand_poly, FE};
 use const_random::const_random;
 use core::hint::black_box;
 use criterion::Criterion;
-use lambdaworks_math::polynomial::Polynomial;
+use lambdaworks_math::{
+    field::fields::mersenne31::extensions::Degree2ExtensionField, polynomial::Polynomial,
+};
 
 pub fn polynomial_benchmarks(c: &mut Criterion) {
     let mut group = c.benchmark_group("Polynomial");
@@ -42,6 +44,25 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
         let y_poly = rand_poly(order);
         bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
     });
+    group.bench_function("fast mul", |bench| {
+        let x_poly = rand_complex_mersenne_poly(order as u32);
+        let y_poly = rand_complex_mersenne_poly(order as u32);
+        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
+    });
+
+    let big_order = 9;
+    let x_poly = rand_complex_mersenne_poly(big_order);
+    let y_poly = rand_complex_mersenne_poly(big_order);
+    group.bench_function("fast_mul big poly", |bench| {
+        bench.iter(|| {
+            black_box(&x_poly)
+                .fast_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
+                .unwrap()
+        });
+    });
+    group.bench_function("slow mul big poly", |bench| {
+        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
+    });
 
     group.bench_function("div", |bench| {
         let x_poly = rand_poly(order);
diff --git a/crates/math/benches/polynomials/utils.rs b/crates/math/benches/polynomials/utils.rs
index 1616906de..f4ae16f9b 100644
--- a/crates/math/benches/polynomials/utils.rs
+++ b/crates/math/benches/polynomials/utils.rs
@@ -1,6 +1,12 @@
 use const_random::const_random;
 use lambdaworks_math::{
-    field::fields::u64_prime_field::{U64FieldElement, U64PrimeField},
+    field::{
+        element::FieldElement,
+        fields::{
+            mersenne31::{extensions::Degree2ExtensionField, field::Mersenne31Field},
+            u64_prime_field::{U64FieldElement, U64PrimeField},
+        },
+    },
     polynomial::{
         dense_multilinear_poly::DenseMultilinearPolynomial,
         sparse_multilinear_poly::SparseMultilinearPolynomial, Polynomial,
@@ -36,6 +42,28 @@ pub fn rand_field_elements(order: u64) -> Vec<FE> {
 pub fn rand_poly(order: u64) -> Polynomial<FE> {
     Polynomial::new(&rand_field_elements(order))
 }
+#[allow(dead_code)]
+#[inline(never)]
+#[export_name = "u64_utils::rand_complex_mersenne_field_elements"]
+pub fn rand_complex_mersenne_field_elements(
+    order: u32,
+) -> Vec<FieldElement<Degree2ExtensionField>> {
+    let mut result = Vec::with_capacity(1 << order);
+    for _ in 0..result.capacity() {
+        result.push(FieldElement::<Degree2ExtensionField>::new([
+            FieldElement::<Mersenne31Field>::new(random()),
+            FieldElement::<Mersenne31Field>::new(random()),
+        ]));
+    }
+    result
+}
+
+#[allow(dead_code)]
+#[inline(never)]
+#[export_name = "u64_utils::rand_complex_mersenne_poly"]
+pub fn rand_complex_mersenne_poly(order: u32) -> Polynomial<FieldElement<Degree2ExtensionField>> {
+    Polynomial::new(&rand_complex_mersenne_field_elements(order))
+}
 
 #[allow(dead_code)]
 #[inline(never)]
diff --git a/crates/math/src/field/fields/mersenne31/extensions.rs b/crates/math/src/field/fields/mersenne31/extensions.rs
index 37773c165..0aad446e8 100644
--- a/crates/math/src/field/fields/mersenne31/extensions.rs
+++ b/crates/math/src/field/fields/mersenne31/extensions.rs
@@ -99,50 +99,6 @@ impl IsFFTField for Degree2ExtensionField {
         [FpE::const_from_raw(2), FpE::const_from_raw(1268011823)];
 }
 
-impl IsSubFieldOf<Degree4ExtensionField> for Degree2ExtensionField {
-    fn add(
-        a: &Self::BaseType,
-        b: &<Degree4ExtensionField as IsField>::BaseType,
-    ) -> <Degree4ExtensionField as IsField>::BaseType {
-        [Fp2E::from(a) + &b[0], b[1].clone()]
-    }
-
-    fn sub(
-        a: &Self::BaseType,
-        b: &<Degree4ExtensionField as IsField>::BaseType,
-    ) -> <Degree4ExtensionField as IsField>::BaseType {
-        [Fp2E::from(a) - &b[0], -&b[1]]
-    }
-
-    fn mul(
-        a: &Self::BaseType,
-        b: &<Degree4ExtensionField as IsField>::BaseType,
-    ) -> <Degree4ExtensionField as IsField>::BaseType {
-        [Fp2E::from(a) * &b[0], Fp2E::from(a) * &b[1]]
-    }
-
-    fn div(
-        a: &Self::BaseType,
-        b: &<Degree4ExtensionField as IsField>::BaseType,
-    ) -> Result<<Degree4ExtensionField as IsField>::BaseType, FieldError> {
-        let b_inv = Degree4ExtensionField::inv(b).map_err(|_| FieldError::DivisionByZero)?;
-        Ok(<Self as IsSubFieldOf<Degree4ExtensionField>>::mul(
-            a, &b_inv,
-        ))
-    }
-
-    fn embed(a: Self::BaseType) -> <Degree4ExtensionField as IsField>::BaseType {
-        [FieldElement::from_raw(a), FieldElement::zero()]
-    }
-
-    #[cfg(feature = "alloc")]
-    fn to_subfield_vec(
-        b: <Degree4ExtensionField as IsField>::BaseType,
-    ) -> alloc::vec::Vec<Self::BaseType> {
-        b.into_iter().map(|x| x.to_raw()).collect()
-    }
-}
-
 impl IsSubFieldOf<Degree2ExtensionField> for Mersenne31Field {
     fn add(
         a: &Self::BaseType,

From d3e5e181ec01b042caa20a0bc1b7a392af6ed2f2 Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Fri, 11 Apr 2025 16:41:23 +0200
Subject: [PATCH 03/10] impl suggestion

---
 crates/math/benches/polynomials/polynomial.rs         | 5 -----
 crates/math/src/fft/errors.rs                         | 4 ++++
 crates/math/src/fft/polynomial.rs                     | 9 ++++++++-
 crates/math/src/field/fields/mersenne31/extensions.rs | 2 ++
 4 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/crates/math/benches/polynomials/polynomial.rs b/crates/math/benches/polynomials/polynomial.rs
index 038cde97b..b3a20b374 100644
--- a/crates/math/benches/polynomials/polynomial.rs
+++ b/crates/math/benches/polynomials/polynomial.rs
@@ -44,11 +44,6 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
         let y_poly = rand_poly(order);
         bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
     });
-    group.bench_function("fast mul", |bench| {
-        let x_poly = rand_complex_mersenne_poly(order as u32);
-        let y_poly = rand_complex_mersenne_poly(order as u32);
-        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
-    });
 
     let big_order = 9;
     let x_poly = rand_complex_mersenne_poly(big_order);
diff --git a/crates/math/src/fft/errors.rs b/crates/math/src/fft/errors.rs
index 1f94986d8..92a122238 100644
--- a/crates/math/src/fft/errors.rs
+++ b/crates/math/src/fft/errors.rs
@@ -10,6 +10,7 @@ pub enum FFTError {
     RootOfUnityError(u64),
     InputError(usize),
     OrderError(u64),
+    DomainSizeError(usize),
     #[cfg(feature = "cuda")]
     CudaError(CudaError),
 }
@@ -24,6 +25,9 @@ impl Display for FFTError {
             FFTError::OrderError(v) => {
                 write!(f, "Order should be less than or equal to 63, but is {v}")
             }
+            FFTError::DomainSizeError(_) => {
+                write!(f, "Domain size exceeds two adicity of the field")
+            }
             #[cfg(feature = "cuda")]
             FFTError::CudaError(_) => {
                 write!(f, "A CUDA related error has ocurred")
diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index 2bb4fcfef..98ca29a88 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -27,7 +27,9 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
     ) -> Result<Vec<FieldElement<E>>, FFTError> {
         let domain_size = domain_size.unwrap_or(0);
         let len = core::cmp::max(poly.coeff_len(), domain_size).next_power_of_two() * blowup_factor;
-
+        if len.trailing_zeros() as u64 > F::TWO_ADICITY {
+            return Err(FFTError::DomainSizeError(len.trailing_zeros() as usize));
+        }
         if poly.coefficients().is_empty() {
             return Ok(vec![FieldElement::zero(); len]);
         }
@@ -425,6 +427,11 @@ mod tests {
                 let (poly, new_poly) = gen_fft_interpolate_and_evaluate(poly);
                 prop_assert_eq!(poly, new_poly);
             }
+
+            #[test]
+            fn test_fft_multiplication_works(poly in poly(7), other in poly(7)) {
+                prop_assert_eq!(poly.fast_multiplication::<F>(&other).unwrap(), poly * other);
+            }
         }
     }
 
diff --git a/crates/math/src/field/fields/mersenne31/extensions.rs b/crates/math/src/field/fields/mersenne31/extensions.rs
index 0aad446e8..1ad7bf95d 100644
--- a/crates/math/src/field/fields/mersenne31/extensions.rs
+++ b/crates/math/src/field/fields/mersenne31/extensions.rs
@@ -94,6 +94,8 @@ impl IsField for Degree2ExtensionField {
 }
 
 impl IsFFTField for Degree2ExtensionField {
+    // Values taken from stwo
+    // https://github.com/starkware-libs/stwo/blob/dev/crates/prover/src/core/circle.rs#L203-L209
     const TWO_ADICITY: u64 = 31;
     const TWO_ADIC_PRIMITVE_ROOT_OF_UNITY: Self::BaseType =
         [FpE::const_from_raw(2), FpE::const_from_raw(1268011823)];

From 2ad4d1d3946fb2a5349828ef44d4c28339f44ee8 Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Tue, 22 Apr 2025 10:59:55 +0200
Subject: [PATCH 04/10] chore: impl suggestion

---
 crates/math/benches/polynomials/polynomial.rs |  2 +-
 crates/math/src/fft/polynomial.rs             | 10 +++++++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/crates/math/benches/polynomials/polynomial.rs b/crates/math/benches/polynomials/polynomial.rs
index b3a20b374..f08f1dc66 100644
--- a/crates/math/benches/polynomials/polynomial.rs
+++ b/crates/math/benches/polynomials/polynomial.rs
@@ -51,7 +51,7 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
     group.bench_function("fast_mul big poly", |bench| {
         bench.iter(|| {
             black_box(&x_poly)
-                .fast_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
+                .fast_fft_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
                 .unwrap()
         });
     });
diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index 98ca29a88..eea2e8344 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -100,7 +100,11 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         Ok(scaled.scale(&offset.inv().unwrap()))
     }
 
-    pub fn fast_multiplication<F: IsFFTField + IsSubFieldOf<E>>(
+    /// Multiplies two polynomials using FFT.
+    /// It's faster than naive multiplication when the degree of the polynomials is large enough (>=2**6).
+    /// This works best with polynomials whose highest degree is equal to a power of 2 - 1.
+    /// Will return an error if the degree of the resulting polynomial is greater than 2**63.
+    pub fn fast_fft_multiplication<F: IsFFTField + IsSubFieldOf<E>>(
         &self,
         other: &Self,
     ) -> Result<Self, FFTError> {
@@ -330,7 +334,7 @@ mod tests {
 
             #[test]
             fn test_fft_multiplication_works(poly in poly(7), other in poly(7)) {
-                prop_assert_eq!(poly.fast_multiplication::<F>(&other).unwrap(), poly * other);
+                prop_assert_eq!(poly.fast_fft_multiplication::<F>(&other).unwrap(), poly * other);
             }
         }
 
@@ -430,7 +434,7 @@ mod tests {
 
             #[test]
             fn test_fft_multiplication_works(poly in poly(7), other in poly(7)) {
-                prop_assert_eq!(poly.fast_multiplication::<F>(&other).unwrap(), poly * other);
+                prop_assert_eq!(poly.fast_fft_multiplication::<F>(&other).unwrap(), poly * other);
             }
         }
     }

From 0f56fc95d8ced63dbc1422e535f76afbc5559d04 Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Thu, 10 Apr 2025 18:18:18 +0200
Subject: [PATCH 05/10] chore: add fast_mul benchmark

---
 crates/math/benches/polynomials/polynomial.rs | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/crates/math/benches/polynomials/polynomial.rs b/crates/math/benches/polynomials/polynomial.rs
index f08f1dc66..a7b993efb 100644
--- a/crates/math/benches/polynomials/polynomial.rs
+++ b/crates/math/benches/polynomials/polynomial.rs
@@ -44,6 +44,25 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
         let y_poly = rand_poly(order);
         bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
     });
+    group.bench_function("fast mul", |bench| {
+        let x_poly = rand_complex_mersenne_poly(order as u32);
+        let y_poly = rand_complex_mersenne_poly(order as u32);
+        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
+    });
+
+    let big_order = 9;
+    let x_poly = rand_complex_mersenne_poly(big_order);
+    let y_poly = rand_complex_mersenne_poly(big_order);
+    group.bench_function("fast_mul big poly", |bench| {
+        bench.iter(|| {
+            black_box(&x_poly)
+                .fast_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
+                .unwrap()
+        });
+    });
+    group.bench_function("slow mul big poly", |bench| {
+        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
+    });
 
     let big_order = 9;
     let x_poly = rand_complex_mersenne_poly(big_order);

From 1e335fb9799f59bb646d57272a7ed9f3d5a931cf Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Fri, 11 Apr 2025 17:49:53 +0200
Subject: [PATCH 06/10] feat(fft): fast division

---
 crates/math/benches/polynomials/polynomial.rs | 26 +++---
 crates/math/src/fft/polynomial.rs             | 87 +++++++++++++++++++
 crates/math/src/polynomial/mod.rs             | 14 +++
 3 files changed, 111 insertions(+), 16 deletions(-)

diff --git a/crates/math/benches/polynomials/polynomial.rs b/crates/math/benches/polynomials/polynomial.rs
index a7b993efb..5b5d31099 100644
--- a/crates/math/benches/polynomials/polynomial.rs
+++ b/crates/math/benches/polynomials/polynomial.rs
@@ -44,11 +44,6 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
         let y_poly = rand_poly(order);
         bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
     });
-    group.bench_function("fast mul", |bench| {
-        let x_poly = rand_complex_mersenne_poly(order as u32);
-        let y_poly = rand_complex_mersenne_poly(order as u32);
-        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
-    });
 
     let big_order = 9;
     let x_poly = rand_complex_mersenne_poly(big_order);
@@ -64,18 +59,17 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
         bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
     });
 
-    let big_order = 9;
-    let x_poly = rand_complex_mersenne_poly(big_order);
-    let y_poly = rand_complex_mersenne_poly(big_order);
-    group.bench_function("fast_mul big poly", |bench| {
-        bench.iter(|| {
-            black_box(&x_poly)
-                .fast_fft_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
-                .unwrap()
-        });
+    let y_poly = rand_complex_mersenne_poly(big_order - 2);
+
+    group.bench_function("fast div big poly", |bench| {
+        let x_poly = rand_complex_mersenne_poly(order as u32);
+        let y_poly = rand_complex_mersenne_poly(order as u32);
+        bench
+            .iter(|| black_box(&x_poly).fast_division::<Degree2ExtensionField>(black_box(&y_poly)));
     });
-    group.bench_function("slow mul big poly", |bench| {
-        bench.iter(|| black_box(&x_poly) * black_box(&y_poly));
+
+    group.bench_function("slow div big poly", |bench| {
+        bench.iter(|| black_box(x_poly.clone()).long_division_with_remainder(black_box(&y_poly)));
     });
 
     group.bench_function("div", |bench| {
diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index eea2e8344..db2ede446 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -115,6 +115,63 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
 
         Polynomial::interpolate_fft::<F>(&r)
     }
+
+    /// Divides two polynomials with remainder.
+    /// This is faster than the naive division if the degree of the divisor
+    /// is greater than the degree of the dividend and both degrees are large enough.
+    pub fn fast_division<F: IsSubFieldOf<E> + IsFFTField>(
+        &self,
+        divisor: &Self,
+    ) -> Result<(Self, Self), FFTError> {
+        let n = self.degree();
+        let m = divisor.degree();
+        if divisor.coefficients.is_empty() {
+            panic!("Division by zero polynomial");
+        }
+        if n < m {
+            return Ok((Self::zero(), self.clone()));
+        }
+        let d = n - m; // Degree of the quotient
+        let a_rev = self.reverse(n);
+        let b_rev = divisor.reverse(m);
+        let inv_b_rev = Self::invert_polynomial::<F>(&b_rev, d + 1)?;
+        let q = a_rev
+            .fast_multiplication::<F>(&inv_b_rev)?
+            .truncate(d + 1)
+            .reverse(d);
+
+        let r = self - q.fast_multiplication::<F>(divisor)?;
+        Ok((q, r))
+    }
+
+    /// Computes the inverse of polynomial P modulo x^k using Newton iteration.
+    /// P must have an invertible constant term.
+    pub fn invert_polynomial<F: IsSubFieldOf<E> + IsFFTField>(
+        p: &Self,
+        k: usize,
+    ) -> Result<Self, FFTError> {
+        if p.coefficients.is_empty() || p.coefficients.iter().all(|c| c == &FieldElement::zero()) {
+            panic!("Cannot invert polynomial with zero constant term");
+        }
+        let mut q = Self::new(&[p.coefficients[0].inv().unwrap()]);
+        let mut current_precision = 1;
+
+        while current_precision < k {
+            let temp = p
+                .fast_multiplication::<F>(&q)?
+                .truncate(2 * current_precision);
+            let two = Self::new(&[FieldElement::<F>::one() + FieldElement::one()]);
+            let correction = two - temp;
+            q = q
+                .fast_multiplication::<F>(&correction)?
+                .truncate(2 * current_precision);
+
+            current_precision *= 2;
+        }
+
+        // Final truncation to desired degree k
+        Ok(q.truncate(k))
+    }
 }
 
 pub fn compose_fft<F, E>(
@@ -273,6 +330,11 @@ mod tests {
                 vec
             }
         }
+        prop_compose! {
+            fn non_empty_field_vec(max_exp: u8)(vec in collection::vec(field_element(), 1 << max_exp)) -> Vec<FE> {
+                vec
+            }
+        }
         prop_compose! {
             fn non_power_of_two_sized_field_vec(max_exp: u8)(vec in collection::vec(field_element(), 2..1<<max_exp).prop_filter("Avoid polynomials of size power of two", |vec| !vec.len().is_power_of_two())) -> Vec<FE> {
                 vec
@@ -283,6 +345,11 @@ mod tests {
                 Polynomial::new(&coeffs)
             }
         }
+        prop_compose! {
+            fn non_zero_poly(max_exp: u8)(coeffs in non_empty_field_vec(max_exp)) -> Polynomial<FE> {
+                Polynomial::new(&coeffs)
+            }
+        }
         prop_compose! {
             fn poly_with_non_power_of_two_coeffs(max_exp: u8)(coeffs in non_power_of_two_sized_field_vec(max_exp)) -> Polynomial<FE> {
                 Polynomial::new(&coeffs)
@@ -336,6 +403,11 @@ mod tests {
             fn test_fft_multiplication_works(poly in poly(7), other in poly(7)) {
                 prop_assert_eq!(poly.fast_fft_multiplication::<F>(&other).unwrap(), poly * other);
             }
+
+            #[test]
+            fn test_fft_division_works(poly in non_zero_poly(7), other in non_zero_poly(7)) {
+                prop_assert_eq!(poly.fast_division::<F>(&other).unwrap(), poly.long_division_with_remainder(&other));
+            }
         }
 
         #[test]
@@ -371,6 +443,11 @@ mod tests {
                 vec
             }
         }
+        prop_compose! {
+            fn non_empty_field_vec(max_exp: u8)(vec in collection::vec(field_element(), 1 << max_exp)) -> Vec<FE> {
+                vec
+            }
+        }
         prop_compose! {
             fn non_power_of_two_sized_field_vec(max_exp: u8)(vec in collection::vec(field_element(), 2..1<<max_exp).prop_filter("Avoid polynomials of size power of two", |vec| !vec.len().is_power_of_two())) -> Vec<FE> {
                 vec
@@ -381,6 +458,11 @@ mod tests {
                 Polynomial::new(&coeffs)
             }
         }
+        prop_compose! {
+            fn non_zero_poly(max_exp: u8)(coeffs in non_empty_field_vec(max_exp)) -> Polynomial<FE> {
+                Polynomial::new(&coeffs)
+            }
+        }
         prop_compose! {
             fn poly_with_non_power_of_two_coeffs(max_exp: u8)(coeffs in non_power_of_two_sized_field_vec(max_exp)) -> Polynomial<FE> {
                 Polynomial::new(&coeffs)
@@ -436,6 +518,11 @@ mod tests {
             fn test_fft_multiplication_works(poly in poly(7), other in poly(7)) {
                 prop_assert_eq!(poly.fast_fft_multiplication::<F>(&other).unwrap(), poly * other);
             }
+
+            #[test]
+            fn test_fft_division_works(poly in poly(7), other in non_zero_poly(7)) {
+                prop_assert_eq!(poly.fast_division::<F>(&other).unwrap(), poly.long_division_with_remainder(&other));
+            }
         }
     }
 
diff --git a/crates/math/src/polynomial/mod.rs b/crates/math/src/polynomial/mod.rs
index e11cd8214..981a65087 100644
--- a/crates/math/src/polynomial/mod.rs
+++ b/crates/math/src/polynomial/mod.rs
@@ -317,6 +317,20 @@ impl<F: IsField> Polynomial<FieldElement<F>> {
                 .collect(),
         }
     }
+
+    pub fn truncate(&self, k: usize) -> Self {
+        if k == 0 {
+            Self::zero()
+        } else {
+            Self::new(&self.coefficients[0..k.min(self.coefficients.len())])
+        }
+    }
+    pub fn reverse(&self, d: usize) -> Self {
+        let mut coeffs = self.coefficients.clone();
+        coeffs.resize(d + 1, FieldElement::zero());
+        coeffs.reverse();
+        Self::new(&coeffs)
+    }
 }
 
 impl<F: IsPrimeField> Polynomial<FieldElement<F>> {

From 28a868b1798ee930bea618bc76c811bfa3da720c Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Mon, 28 Apr 2025 13:42:07 +0200
Subject: [PATCH 07/10] refacto: proper error handling

---
 crates/math/benches/polynomials/polynomial.rs |  4 +--
 crates/math/src/fft/polynomial.rs             | 33 +++++++++++--------
 2 files changed, 20 insertions(+), 17 deletions(-)

diff --git a/crates/math/benches/polynomials/polynomial.rs b/crates/math/benches/polynomials/polynomial.rs
index 5b5d31099..67cb12c51 100644
--- a/crates/math/benches/polynomials/polynomial.rs
+++ b/crates/math/benches/polynomials/polynomial.rs
@@ -51,7 +51,7 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
     group.bench_function("fast_mul big poly", |bench| {
         bench.iter(|| {
             black_box(&x_poly)
-                .fast_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
+                .fast_fft_multiplication::<Degree2ExtensionField>(black_box(&y_poly))
                 .unwrap()
         });
     });
@@ -62,8 +62,6 @@ pub fn polynomial_benchmarks(c: &mut Criterion) {
     let y_poly = rand_complex_mersenne_poly(big_order - 2);
 
     group.bench_function("fast div big poly", |bench| {
-        let x_poly = rand_complex_mersenne_poly(order as u32);
-        let y_poly = rand_complex_mersenne_poly(order as u32);
         bench
             .iter(|| black_box(&x_poly).fast_division::<Degree2ExtensionField>(black_box(&y_poly)));
     });
diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index db2ede446..41c998faf 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -1,5 +1,6 @@
 use crate::fft::errors::FFTError;
 
+use crate::field::errors::FieldError;
 use crate::field::traits::{IsField, IsSubFieldOf};
 use crate::{
     field::{
@@ -125,8 +126,13 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
     ) -> Result<(Self, Self), FFTError> {
         let n = self.degree();
         let m = divisor.degree();
-        if divisor.coefficients.is_empty() {
-            panic!("Division by zero polynomial");
+        if divisor.coefficients.is_empty()
+            || divisor
+                .coefficients
+                .iter()
+                .all(|c| c == &FieldElement::zero())
+        {
+            return Err(FieldError::DivisionByZero.into());
         }
         if n < m {
             return Ok((Self::zero(), self.clone()));
@@ -136,11 +142,11 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         let b_rev = divisor.reverse(m);
         let inv_b_rev = Self::invert_polynomial::<F>(&b_rev, d + 1)?;
         let q = a_rev
-            .fast_multiplication::<F>(&inv_b_rev)?
+            .fast_fft_multiplication::<F>(&inv_b_rev)?
             .truncate(d + 1)
             .reverse(d);
 
-        let r = self - q.fast_multiplication::<F>(divisor)?;
+        let r = self - q.fast_fft_multiplication::<F>(divisor)?;
         Ok((q, r))
     }
 
@@ -151,22 +157,21 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         k: usize,
     ) -> Result<Self, FFTError> {
         if p.coefficients.is_empty() || p.coefficients.iter().all(|c| c == &FieldElement::zero()) {
-            panic!("Cannot invert polynomial with zero constant term");
+            return Err(FieldError::DivisionByZero.into());
         }
-        let mut q = Self::new(&[p.coefficients[0].inv().unwrap()]);
+        let mut q = Self::new(&[p.coefficients[0].inv()?]);
         let mut current_precision = 1;
 
+        let two = Self::new(&[FieldElement::<F>::one() + FieldElement::one()]);
         while current_precision < k {
+            current_precision *= 2;
             let temp = p
-                .fast_multiplication::<F>(&q)?
-                .truncate(2 * current_precision);
-            let two = Self::new(&[FieldElement::<F>::one() + FieldElement::one()]);
-            let correction = two - temp;
+                .fast_fft_multiplication::<F>(&q)?
+                .truncate(current_precision);
+            let correction = &two - temp;
             q = q
-                .fast_multiplication::<F>(&correction)?
-                .truncate(2 * current_precision);
-
-            current_precision *= 2;
+                .fast_fft_multiplication::<F>(&correction)?
+                .truncate(current_precision);
         }
 
         // Final truncation to desired degree k

From 40611cd3027dece645f1bf9717c4e0391667643f Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Mon, 28 Apr 2025 13:49:42 +0200
Subject: [PATCH 08/10] refacto: rename inversion function

---
 crates/math/src/fft/polynomial.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index 41c998faf..4b1eb5b2f 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -140,7 +140,7 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         let d = n - m; // Degree of the quotient
         let a_rev = self.reverse(n);
         let b_rev = divisor.reverse(m);
-        let inv_b_rev = Self::invert_polynomial::<F>(&b_rev, d + 1)?;
+        let inv_b_rev = Self::invert_polynomial_mod::<F>(&b_rev, d + 1)?;
         let q = a_rev
             .fast_fft_multiplication::<F>(&inv_b_rev)?
             .truncate(d + 1)
@@ -152,7 +152,7 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
 
     /// Computes the inverse of polynomial P modulo x^k using Newton iteration.
     /// P must have an invertible constant term.
-    pub fn invert_polynomial<F: IsSubFieldOf<E> + IsFFTField>(
+    pub fn invert_polynomial_mod<F: IsSubFieldOf<E> + IsFFTField>(
         p: &Self,
         k: usize,
     ) -> Result<Self, FFTError> {

From 8f7614980c55421f912bd462e1395ed39fb00df5 Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Thu, 15 May 2025 13:33:10 +0200
Subject: [PATCH 09/10] test: add tests for fast division helper functions

---
 crates/math/src/fft/polynomial.rs | 24 +++++++++++++++++++-----
 crates/math/src/polynomial/mod.rs | 15 +++++++++++++++
 2 files changed, 34 insertions(+), 5 deletions(-)

diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index 4b1eb5b2f..5e94299b9 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -140,7 +140,7 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         let d = n - m; // Degree of the quotient
         let a_rev = self.reverse(n);
         let b_rev = divisor.reverse(m);
-        let inv_b_rev = Self::invert_polynomial_mod::<F>(&b_rev, d + 1)?;
+        let inv_b_rev = b_rev.invert_polynomial_mod::<F>(d + 1)?;
         let q = a_rev
             .fast_fft_multiplication::<F>(&inv_b_rev)?
             .truncate(d + 1)
@@ -153,19 +153,21 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
     /// Computes the inverse of polynomial P modulo x^k using Newton iteration.
     /// P must have an invertible constant term.
     pub fn invert_polynomial_mod<F: IsSubFieldOf<E> + IsFFTField>(
-        p: &Self,
+        &self,
         k: usize,
     ) -> Result<Self, FFTError> {
-        if p.coefficients.is_empty() || p.coefficients.iter().all(|c| c == &FieldElement::zero()) {
+        if self.coefficients.is_empty()
+            || self.coefficients.iter().all(|c| c == &FieldElement::zero())
+        {
             return Err(FieldError::DivisionByZero.into());
         }
-        let mut q = Self::new(&[p.coefficients[0].inv()?]);
+        let mut q = Self::new(&[self.coefficients[0].inv()?]);
         let mut current_precision = 1;
 
         let two = Self::new(&[FieldElement::<F>::one() + FieldElement::one()]);
         while current_precision < k {
             current_precision *= 2;
-            let temp = p
+            let temp = self
                 .fast_fft_multiplication::<F>(&q)?
                 .truncate(current_precision);
             let correction = &two - temp;
@@ -413,6 +415,12 @@ mod tests {
             fn test_fft_division_works(poly in non_zero_poly(7), other in non_zero_poly(7)) {
                 prop_assert_eq!(poly.fast_division::<F>(&other).unwrap(), poly.long_division_with_remainder(&other));
             }
+
+            #[test]
+            fn test_invert_polynomial_mod_works(poly in non_zero_poly(7), k in powers_of_two(4)) {
+                let inverted_poly = poly.invert_polynomial_mod::<F>(k).unwrap();
+                prop_assert_eq!((poly * inverted_poly).truncate(k), Polynomial::new(&[FE::one()]));
+            }
         }
 
         #[test]
@@ -528,6 +536,12 @@ mod tests {
             fn test_fft_division_works(poly in poly(7), other in non_zero_poly(7)) {
                 prop_assert_eq!(poly.fast_division::<F>(&other).unwrap(), poly.long_division_with_remainder(&other));
             }
+
+            #[test]
+            fn test_invert_polynomial_mod_works(poly in non_zero_poly(7), k in powers_of_two(4)) {
+                let inverted_poly = poly.invert_polynomial_mod::<F>(k).unwrap();
+                prop_assert_eq!((poly * inverted_poly).truncate(k), Polynomial::new(&[FE::one()]));
+            }
         }
     }
 
diff --git a/crates/math/src/polynomial/mod.rs b/crates/math/src/polynomial/mod.rs
index 70607c536..109d54819 100644
--- a/crates/math/src/polynomial/mod.rs
+++ b/crates/math/src/polynomial/mod.rs
@@ -1291,6 +1291,21 @@ mod tests {
         assert_eq!(dpdx, Polynomial::new(&[FE::new(0)]));
     }
 
+    #[test]
+    fn test_reverse() {
+        let p = Polynomial::new(&[FE::new(3), FE::new(2), FE::new(1)]);
+        assert_eq!(
+            p.reverse(3),
+            Polynomial::new(&[FE::new(0), FE::new(1), FE::new(2), FE::new(3)])
+        );
+    }
+
+    #[test]
+    fn test_truncate() {
+        let p = Polynomial::new(&[FE::new(3), FE::new(2), FE::new(1)]);
+        assert_eq!(p.truncate(2), Polynomial::new(&[FE::new(3), FE::new(2)]));
+    }
+
     #[test]
     fn test_print_as_sage_poly() {
         let p = Polynomial::new(&[FE::new(1), FE::new(2), FE::new(3)]);

From bb9e29ccb9c46f5805d6fde5135b6345f6d23fe0 Mon Sep 17 00:00:00 2001
From: 0xLucqs <lucaslevy10@gmail.com>
Date: Thu, 15 May 2025 13:35:11 +0200
Subject: [PATCH 10/10] doc: explain where fast division is taken from

---
 crates/math/src/fft/polynomial.rs | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/crates/math/src/fft/polynomial.rs b/crates/math/src/fft/polynomial.rs
index 5e94299b9..c966d8236 100644
--- a/crates/math/src/fft/polynomial.rs
+++ b/crates/math/src/fft/polynomial.rs
@@ -105,6 +105,10 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
     /// It's faster than naive multiplication when the degree of the polynomials is large enough (>=2**6).
     /// This works best with polynomials whose highest degree is equal to a power of 2 - 1.
     /// Will return an error if the degree of the resulting polynomial is greater than 2**63.
+    ///
+    /// This is an implementation of the fast division algorithm from
+    /// [Gathen's book](https://www.cambridge.org/core/books/modern-computer-algebra/DB3563D4013401734851CF683D2F03F0)
+    /// chapter 9
     pub fn fast_fft_multiplication<F: IsFFTField + IsSubFieldOf<E>>(
         &self,
         other: &Self,