diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2acaf9643..237a72262 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,7 @@
 - [\#772](https://github.com/arkworks-rs/algebra/pull/772) (`ark-ff`) Implementation of `mul` method for `BigInteger`.
 - [\#794](https://github.com/arkworks-rs/algebra/pull/794) (`ark-ff`) Fix `wasm` compilation.
 - [\#837](https://github.com/arkworks-rs/algebra/pull/837) (`ark-serialize`) Fix array deserialization panic.
+- [\#845](https://github.com/arkworks-rs/algebra/pull/845) (`Algebra`) Implementation of `mul` method for `DenseMultilinearExtension<F> * F`.
 
 ### Breaking changes
 
diff --git a/Cargo.toml b/Cargo.toml
index c0d3197ac..72b72eb1c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -56,10 +56,10 @@ arrayvec = { version = "0.7", default-features = false }
 criterion = "0.5.0"
 educe = "0.6.0"
 digest = { version = "0.10", default-features = false }
-hashbrown = { version = "0.14", default-features = false, features = ["inline-more", "allocator-api2"] }
+hashbrown = { version = "0.15", default-features = false, features = ["inline-more", "allocator-api2"] }
 hex = "0.4"
 itertools = { version = "0.13", default-features = false }
-libtest-mimic = "0.7.0"
+libtest-mimic = "0.8.1"
 paste = "1.0"
 rayon = "1"
 serde = "1.0"
diff --git a/curves/bn254/src/curves/mod.rs b/curves/bn254/src/curves/mod.rs
index 349488786..68e15233b 100644
--- a/curves/bn254/src/curves/mod.rs
+++ b/curves/bn254/src/curves/mod.rs
@@ -19,9 +19,9 @@ impl BnConfig for Config {
     /// `x` is positive.
     const X_IS_NEGATIVE: bool = false;
     const ATE_LOOP_COUNT: &'static [i8] = &[
-        0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, -1, 0, 0, 0,
-        0, 1, 1, 1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 0, 1, 1, 0,
-        -1, 0, 0, 1, 0, 1, 1,
+        0, 0, 0, 1, 0, 1, 0, -1, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, -1, 0, 0, 0,
+        0, -1, 0, 0, 1, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 1, 0, -1, 0, 0, 0, -1, 0,
+        -1, 0, 0, 0, 1, 0, 1, 1,
     ];
 
     const TWIST_MUL_BY_Q_X: Fq2 = Fq2::new(
diff --git a/ff/src/fields/models/cubic_extension.rs b/ff/src/fields/models/cubic_extension.rs
index 8908a5afa..1049908bf 100644
--- a/ff/src/fields/models/cubic_extension.rs
+++ b/ff/src/fields/models/cubic_extension.rs
@@ -1,6 +1,7 @@
 use crate::{
     fields::{Field, PrimeField},
-    AdditiveGroup, LegendreSymbol, One, SqrtPrecomputation, ToConstraintField, UniformRand, Zero,
+    AdditiveGroup, FftField, LegendreSymbol, One, SqrtPrecomputation, ToConstraintField,
+    UniformRand, Zero,
 };
 use ark_serialize::{
     CanonicalDeserialize, CanonicalDeserializeWithFlags, CanonicalSerialize,
@@ -770,3 +771,28 @@ mod cube_ext_tests {
         }
     }
 }
+
+impl<P: CubicExtConfig> FftField for CubicExtField<P>
+where
+    P::BaseField: FftField,
+{
+    const GENERATOR: Self = Self::new(
+        P::BaseField::GENERATOR,
+        P::BaseField::ZERO,
+        P::BaseField::ZERO,
+    );
+    const TWO_ADICITY: u32 = P::BaseField::TWO_ADICITY;
+    const TWO_ADIC_ROOT_OF_UNITY: Self = Self::new(
+        P::BaseField::TWO_ADIC_ROOT_OF_UNITY,
+        P::BaseField::ZERO,
+        P::BaseField::ZERO,
+    );
+    const SMALL_SUBGROUP_BASE: Option<u32> = P::BaseField::SMALL_SUBGROUP_BASE;
+    const SMALL_SUBGROUP_BASE_ADICITY: Option<u32> = P::BaseField::SMALL_SUBGROUP_BASE_ADICITY;
+    const LARGE_SUBGROUP_ROOT_OF_UNITY: Option<Self> =
+        if let Some(x) = P::BaseField::LARGE_SUBGROUP_ROOT_OF_UNITY {
+            Some(Self::new(x, P::BaseField::ZERO, P::BaseField::ZERO))
+        } else {
+            None
+        };
+}
diff --git a/ff/src/fields/models/fp3.rs b/ff/src/fields/models/fp3.rs
index 3d44f8fe3..96007f3fa 100644
--- a/ff/src/fields/models/fp3.rs
+++ b/ff/src/fields/models/fp3.rs
@@ -39,7 +39,6 @@ impl<P: Fp3Config> CubicExtConfig for Fp3ConfigWrapper<P> {
     type FrobCoeff = P::Fp;
 
     const DEGREE_OVER_BASE_PRIME_FIELD: usize = 3;
-
     const NONRESIDUE: Self::BaseField = P::NONRESIDUE;
 
     const SQRT_PRECOMP: Option<SqrtPrecomputation<CubicExtField<Self>>> =
diff --git a/ff/src/fields/models/quadratic_extension.rs b/ff/src/fields/models/quadratic_extension.rs
index d89feddae..4ebbb5a62 100644
--- a/ff/src/fields/models/quadratic_extension.rs
+++ b/ff/src/fields/models/quadratic_extension.rs
@@ -1,7 +1,7 @@
 use crate::{
     biginteger::BigInteger,
     fields::{Field, LegendreSymbol, PrimeField},
-    AdditiveGroup, One, SqrtPrecomputation, ToConstraintField, UniformRand, Zero,
+    AdditiveGroup, FftField, One, SqrtPrecomputation, ToConstraintField, UniformRand, Zero,
 };
 use ark_serialize::{
     CanonicalDeserialize, CanonicalDeserializeWithFlags, CanonicalSerialize,
@@ -825,3 +825,21 @@ mod quad_ext_tests {
         }
     }
 }
+
+impl<P: QuadExtConfig> FftField for QuadExtField<P>
+where
+    P::BaseField: FftField,
+{
+    const GENERATOR: Self = Self::new(P::BaseField::GENERATOR, P::BaseField::ZERO);
+    const TWO_ADICITY: u32 = P::BaseField::TWO_ADICITY;
+    const TWO_ADIC_ROOT_OF_UNITY: Self =
+        Self::new(P::BaseField::TWO_ADIC_ROOT_OF_UNITY, P::BaseField::ZERO);
+    const SMALL_SUBGROUP_BASE: Option<u32> = P::BaseField::SMALL_SUBGROUP_BASE;
+    const SMALL_SUBGROUP_BASE_ADICITY: Option<u32> = P::BaseField::SMALL_SUBGROUP_BASE_ADICITY;
+    const LARGE_SUBGROUP_ROOT_OF_UNITY: Option<Self> =
+        if let Some(x) = P::BaseField::LARGE_SUBGROUP_ROOT_OF_UNITY {
+            Some(Self::new(x, P::BaseField::ZERO))
+        } else {
+            None
+        };
+}
diff --git a/poly/benches/fft.rs b/poly/benches/fft.rs
index 0697f768e..47ccb0dbc 100644
--- a/poly/benches/fft.rs
+++ b/poly/benches/fft.rs
@@ -12,7 +12,7 @@ use criterion::{criterion_group, criterion_main, Bencher, BenchmarkId, Criterion
 
 // degree bounds to benchmark on
 // e.g. degree bound of 2^{15}, means we do an FFT for a degree (2^{15} - 1) polynomial
-const BENCHMARK_MIN_DEGREE: usize = 1 << 15;
+const BENCHMARK_MIN_DEGREE: usize = 1 << 4;
 const BENCHMARK_MAX_DEGREE_BLS12_381: usize = 1 << 22;
 const BENCHMARK_MAX_DEGREE_MNT6_753: usize = 1 << 17;
 const BENCHMARK_LOG_INTERVAL_DEGREE: usize = 1;
diff --git a/poly/src/domain/radix2/fft.rs b/poly/src/domain/radix2/fft.rs
index 15d5451a6..875a12414 100644
--- a/poly/src/domain/radix2/fft.rs
+++ b/poly/src/domain/radix2/fft.rs
@@ -224,23 +224,33 @@ impl<F: FftField> Radix2EvaluationDomain<F> {
         max_threads: usize,
         gap: usize,
     ) {
-        cfg_chunks_mut!(xi, chunk_size).for_each(|cxi| {
-            let (lo, hi) = cxi.split_at_mut(gap);
-            // If the chunk is sufficiently big that parallelism helps,
-            // we parallelize the butterfly operation within the chunk.
-
-            if gap > MIN_GAP_SIZE_FOR_PARALLELISATION && num_chunks < max_threads {
-                cfg_iter_mut!(lo)
-                    .zip(hi)
-                    .zip(cfg_iter!(roots).step_by(step))
-                    .for_each(g);
-            } else {
+        if xi.len() <= MIN_INPUT_SIZE_FOR_PARALLELIZATION {
+            xi.chunks_mut(chunk_size).for_each(|cxi| {
+                let (lo, hi) = cxi.split_at_mut(gap);
                 lo.iter_mut()
                     .zip(hi)
                     .zip(roots.iter().step_by(step))
                     .for_each(g);
-            }
-        });
+            });
+        } else {
+            cfg_chunks_mut!(xi, chunk_size).for_each(|cxi| {
+                let (lo, hi) = cxi.split_at_mut(gap);
+                // If the chunk is sufficiently big that parallelism helps,
+                // we parallelize the butterfly operation within the chunk.
+
+                if gap > MIN_GAP_SIZE_FOR_PARALLELIZATION && num_chunks < max_threads {
+                    cfg_iter_mut!(lo)
+                        .zip(hi)
+                        .zip(cfg_iter!(roots).step_by(step))
+                        .for_each(g);
+                } else {
+                    lo.iter_mut()
+                        .zip(hi)
+                        .zip(roots.iter().step_by(step))
+                        .for_each(g);
+                }
+            });
+        }
     }
 
     fn io_helper<T: DomainCoeff<F>>(&self, xi: &mut [T], root: F) {
@@ -349,7 +359,11 @@ const MIN_NUM_CHUNKS_FOR_COMPACTION: usize = 1 << 7;
 
 /// The minimum size of a chunk at which parallelization of `butterfly`s is
 /// beneficial. This value was chosen empirically.
-const MIN_GAP_SIZE_FOR_PARALLELISATION: usize = 1 << 10;
+const MIN_GAP_SIZE_FOR_PARALLELIZATION: usize = 1 << 10;
+
+/// The minimum size of a chunk at which parallelization of `butterfly`s is
+/// beneficial. This value was chosen empirically.
+const MIN_INPUT_SIZE_FOR_PARALLELIZATION: usize = 1 << 10;
 
 // minimum size at which to parallelize.
 #[cfg(feature = "parallel")]
diff --git a/poly/src/evaluations/multivariate/multilinear/dense.rs b/poly/src/evaluations/multivariate/multilinear/dense.rs
index 9e7fe482e..315c54e28 100644
--- a/poly/src/evaluations/multivariate/multilinear/dense.rs
+++ b/poly/src/evaluations/multivariate/multilinear/dense.rs
@@ -11,7 +11,7 @@ use ark_std::{
     fmt::Formatter,
     iter::IntoIterator,
     log2,
-    ops::{Add, AddAssign, Index, Neg, Sub, SubAssign},
+    ops::{Add, AddAssign, Index, Mul, MulAssign, Neg, Sub, SubAssign},
     rand::Rng,
     slice::{Iter, IterMut},
     vec::*,
@@ -331,6 +331,44 @@ impl<'a, F: Field> SubAssign<&'a DenseMultilinearExtension<F>> for DenseMultilin
     }
 }
 
+impl<F: Field> Mul<F> for DenseMultilinearExtension<F> {
+    type Output = DenseMultilinearExtension<F>;
+
+    fn mul(self, scalar: F) -> Self::Output {
+        &self * &scalar
+    }
+}
+
+impl<'a, 'b, F: Field> Mul<&'a F> for &'b DenseMultilinearExtension<F> {
+    type Output = DenseMultilinearExtension<F>;
+
+    fn mul(self, scalar: &'a F) -> Self::Output {
+        if scalar.is_zero() {
+            return DenseMultilinearExtension::zero();
+        } else if scalar.is_one() {
+            return self.clone();
+        }
+        let result: Vec<F> = self.evaluations.iter().map(|&x| x * scalar).collect();
+
+        DenseMultilinearExtension {
+            num_vars: self.num_vars,
+            evaluations: result,
+        }
+    }
+}
+
+impl<F: Field> MulAssign<F> for DenseMultilinearExtension<F> {
+    fn mul_assign(&mut self, scalar: F) {
+        *self = &*self * &scalar
+    }
+}
+
+impl<'a, F: Field> MulAssign<&'a F> for DenseMultilinearExtension<F> {
+    fn mul_assign(&mut self, scalar: &'a F) {
+        *self = &*self * scalar
+    }
+}
+
 impl<F: Field> fmt::Debug for DenseMultilinearExtension<F> {
     fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), fmt::Error> {
         write!(f, "DenseML(nv = {}, evaluations = [", self.num_vars)?;
@@ -394,7 +432,7 @@ impl<F: Field> Polynomial<F> for DenseMultilinearExtension<F> {
 #[cfg(test)]
 mod tests {
     use crate::{DenseMultilinearExtension, MultilinearExtension, Polynomial};
-    use ark_ff::{Field, Zero};
+    use ark_ff::{Field, One, Zero};
     use ark_std::{ops::Neg, test_rng, vec::*, UniformRand};
     use ark_test_curves::bls12_381::Fr;
 
@@ -471,6 +509,7 @@ mod tests {
         const NV: usize = 10;
         let mut rng = test_rng();
         for _ in 0..20 {
+            let scalar = Fr::rand(&mut rng);
             let point: Vec<_> = (0..NV).map(|_| Fr::rand(&mut rng)).collect();
             let poly1 = DenseMultilinearExtension::rand(NV, &mut rng);
             let poly2 = DenseMultilinearExtension::rand(NV, &mut rng);
@@ -482,6 +521,8 @@ mod tests {
             assert_eq!((&poly1 - &poly2).evaluate(&point), v1 - v2);
             // test negate
             assert_eq!(poly1.clone().neg().evaluate(&point), -v1);
+            // test mul poly by scalar
+            assert_eq!((&poly1 * &scalar).evaluate(&point), v1 * scalar);
             // test add assign
             {
                 let mut poly1 = poly1.clone();
@@ -515,6 +556,16 @@ mod tests {
                     assert_eq!(zero.evaluate(&point), scalar * v1);
                 }
             }
+            // test mul_assign for poly * scalar
+            {
+                let mut poly1_cloned = poly1.clone();
+                poly1_cloned *= Fr::one();
+                assert_eq!(poly1_cloned.evaluate(&point), v1);
+                poly1_cloned *= scalar;
+                assert_eq!(poly1_cloned.evaluate(&point), v1 * scalar);
+                poly1_cloned *= Fr::zero();
+                assert_eq!(poly1_cloned, DenseMultilinearExtension::zero());
+            }
         }
     }
 
diff --git a/poly/src/polynomial/univariate/dense.rs b/poly/src/polynomial/univariate/dense.rs
index 1df019d58..5ea41b586 100644
--- a/poly/src/polynomial/univariate/dense.rs
+++ b/poly/src/polynomial/univariate/dense.rs
@@ -170,12 +170,12 @@ impl<F: FftField> DensePolynomial<F> {
     pub fn divide_by_vanishing_poly<D: EvaluationDomain<F>>(
         &self,
         domain: D,
-    ) -> Option<(DensePolynomial<F>, DensePolynomial<F>)> {
+    ) -> (DensePolynomial<F>, DensePolynomial<F>) {
         let domain_size = domain.size();
 
         if self.coeffs.len() < domain_size {
             // If degree(self) < len(Domain), then the quotient is zero, and the entire polynomial is the remainder
-            Some((DensePolynomial::<F>::zero(), self.clone()))
+            (DensePolynomial::<F>::zero(), self.clone())
         } else {
             // Compute the quotient
             //
@@ -211,7 +211,7 @@ impl<F: FftField> DensePolynomial<F> {
 
             let quotient = DensePolynomial::<F>::from_coefficients_vec(quotient_vec);
             let remainder = DensePolynomial::<F>::from_coefficients_vec(remainder_vec);
-            Some((quotient, remainder))
+            (quotient, remainder)
         }
     }
 }
@@ -285,14 +285,6 @@ impl<F: Field> DerefMut for DensePolynomial<F> {
     }
 }
 
-impl<F: Field> Add for DensePolynomial<F> {
-    type Output = DensePolynomial<F>;
-
-    fn add(self, other: DensePolynomial<F>) -> Self {
-        &self + &other
-    }
-}
-
 impl<'a, 'b, F: Field> Add<&'a DensePolynomial<F>> for &'b DensePolynomial<F> {
     type Output = DensePolynomial<F>;
 
@@ -601,6 +593,15 @@ impl<'b, F: Field> Mul<F> for &'b DensePolynomial<F> {
     }
 }
 
+impl<F: Field> Mul<F> for DensePolynomial<F> {
+    type Output = DensePolynomial<F>;
+
+    #[inline]
+    fn mul(self, elem: F) -> DensePolynomial<F> {
+        &self * elem
+    }
+}
+
 /// Performs O(nlogn) multiplication of polynomials if F is smooth.
 impl<'a, 'b, F: FftField> Mul<&'a DensePolynomial<F>> for &'b DensePolynomial<F> {
     type Output = DensePolynomial<F>;
@@ -620,6 +621,37 @@ impl<'a, 'b, F: FftField> Mul<&'a DensePolynomial<F>> for &'b DensePolynomial<F>
     }
 }
 
+macro_rules! impl_op {
+    ($trait:ident, $method:ident, $field_bound:ident) => {
+        impl<F: $field_bound> $trait<DensePolynomial<F>> for DensePolynomial<F> {
+            type Output = DensePolynomial<F>;
+
+            #[inline]
+            fn $method(self, other: DensePolynomial<F>) -> DensePolynomial<F> {
+                (&self).$method(&other)
+            }
+        }
+
+        impl<'a, F: $field_bound> $trait<&'a DensePolynomial<F>> for DensePolynomial<F> {
+            type Output = DensePolynomial<F>;
+
+            #[inline]
+            fn $method(self, other: &'a DensePolynomial<F>) -> DensePolynomial<F> {
+                (&self).$method(other)
+            }
+        }
+
+        impl<'a, F: $field_bound> $trait<DensePolynomial<F>> for &'a DensePolynomial<F> {
+            type Output = DensePolynomial<F>;
+
+            #[inline]
+            fn $method(self, other: DensePolynomial<F>) -> DensePolynomial<F> {
+                self.$method(&other)
+            }
+        }
+    };
+}
+
 impl<F: Field> Zero for DensePolynomial<F> {
     /// Returns the zero polynomial.
     fn zero() -> Self {
@@ -632,6 +664,11 @@ impl<F: Field> Zero for DensePolynomial<F> {
     }
 }
 
+impl_op!(Add, add, Field);
+impl_op!(Sub, sub, Field);
+impl_op!(Mul, mul, FftField);
+impl_op!(Div, div, Field);
+
 #[cfg(test)]
 mod tests {
     use crate::{polynomial::univariate::*, GeneralEvaluationDomain};
@@ -899,7 +936,7 @@ mod tests {
             let domain = GeneralEvaluationDomain::new(1 << size).unwrap();
             for degree in 0..12 {
                 let p = DensePolynomial::<Fr>::rand(degree * 100, rng);
-                let (quotient, remainder) = p.divide_by_vanishing_poly(domain).unwrap();
+                let (quotient, remainder) = p.divide_by_vanishing_poly(domain);
                 let p_recovered = quotient.mul_by_vanishing_poly(domain) + remainder;
                 assert_eq!(p, p_recovered);
             }
diff --git a/poly/src/polynomial/univariate/sparse.rs b/poly/src/polynomial/univariate/sparse.rs
index a69bf3fb2..313bb1f5a 100644
--- a/poly/src/polynomial/univariate/sparse.rs
+++ b/poly/src/polynomial/univariate/sparse.rs
@@ -538,8 +538,7 @@ mod tests {
 
                 // Test interpolation works, by checking that interpolated polynomial agrees with the original on the domain
                 let (_q, r) = (dense_poly.clone() + -sparse_evals.interpolate())
-                    .divide_by_vanishing_poly(domain)
-                    .unwrap();
+                    .divide_by_vanishing_poly(domain);
                 assert_eq!(
                     r,
                     DensePolynomial::<Fr>::zero(),
@@ -550,8 +549,7 @@ mod tests {
 
                 // Consistency check that the dense polynomials interpolation is correct.
                 let (_q, r) = (dense_poly.clone() + -dense_evals.interpolate())
-                    .divide_by_vanishing_poly(domain)
-                    .unwrap();
+                    .divide_by_vanishing_poly(domain);
                 assert_eq!(
                     r,
                     DensePolynomial::<Fr>::zero(),