Finish merging in old_polynomial

src/polynomial/division.rs

@@ -3,73 +3,140 @@ use crate::field::field::Field;
 use crate::polynomial::polynomial::PolynomialCoeffs;
 use crate::util::log2_strict;
 
-/// Takes a polynomial `a` in coefficient form, and divides it by `Z_H = X^n - 1`.
+impl<F: Field> PolynomialCoeffs<F> {
-///
+    /// Polynomial division.
-/// This assumes `Z_H | a`, otherwise result is meaningless.
+    /// Returns `(q, r)`, the quotient and remainder of the polynomial division of `a` by `b`.
-pub(crate) fn divide_by_z_h<F: Field>(mut a: PolynomialCoeffs<F>, n: usize) -> PolynomialCoeffs<F> {
+    pub(crate) fn div_rem(&self, b: &Self) -> (Self, Self) {
-    // TODO: Is this special case needed?
+        let (a_degree_plug_1, b_degree_plus_1) = (self.degree_plus_one(), b.degree_plus_one());
-    if a.coeffs.iter().all(|p| *p == F::ZERO) {
+        if self.is_zero() {
-        return a;
+            (Self::zero(1), Self::empty())
+        } else if b.is_zero() {
+            panic!("Division by zero polynomial");
+        } else if a_degree_plug_1 < b_degree_plus_1 {
+            (Self::zero(1), self.clone())
+        } else if b_degree_plus_1 == 1 {
+            (self * b.coeffs[0].inverse(), Self::empty())
+        } else {
+            let rev_b = b.rev();
+            let rev_b_inv = rev_b.inv_mod_xn(a_degree_plug_1 - b_degree_plus_1 + 1);
+            let rhs: Self = self.rev().coeffs[..=a_degree_plug_1 - b_degree_plus_1]
+                .to_vec()
+                .into();
+            let rev_q: Self = (&rev_b_inv * &rhs).coeffs[..=a_degree_plug_1 - b_degree_plus_1]
+                .to_vec()
+                .into();
+            let mut q = rev_q.rev();
+            let mut qb = &q * b;
+            let mut r = self - &qb;
+            q.trim();
+            r.trim();
+            (q, r)
+        }
     }
 
-    let g = F::MULTIPLICATIVE_GROUP_GENERATOR;
+    /// Polynomial long division.
-    let mut g_pow = F::ONE;
+    /// Returns `(q, r)`, the quotient and remainder of the polynomial division of `a` by `b`.
-    // Multiply the i-th coefficient of `a` by `g^i`. Then `new_a(w^j) = old_a(g.w^j)`.
+    /// Generally slower that the equivalent function `Polynomial::polynomial_division`.
-    a.coeffs.iter_mut().for_each(|x| {
+    pub fn div_rem_long_division(&self, b: &Self) -> (Self, Self) {
-        *x *= g_pow;
+        let b = b.trimmed();
-        g_pow *= g;
-    });
 
-    let root = F::primitive_root_of_unity(log2_strict(a.len()));
+        let (a_degree_plus_1, b_degree_plus_1) = (self.degree_plus_one(), b.degree_plus_one());
-    // Equals to the evaluation of `a` on `{g.w^i}`.
+        if self.is_zero() {
-    let mut a_eval = fft(a);
+            (Self::zero(1), Self::empty())
-    // Compute the denominators `1/(g^n.w^(n*i) - 1)` using batch inversion.
+        } else if b.is_zero() {
-    let denominator_g = g.exp_usize(n);
+            panic!("Division by zero polynomial");
-    let root_n = root.exp_usize(n);
+        } else if a_degree_plus_1 < b_degree_plus_1 {
-    let mut root_pow = F::ONE;
+            (Self::zero(1), self.clone())
-    let denominators = (0..a_eval.len())
+        } else {
-        .map(|i| {
+            // Now we know that self.degree() >= divisor.degree();
-            if i != 0 {
+            let mut quotient = Self::zero(a_degree_plus_1 - b_degree_plus_1 + 1);
-                root_pow = root_pow * root_n;
+            let mut remainder = self.clone();
+            // Can unwrap here because we know self is not zero.
+            let divisor_leading_inv = b.lead().inverse();
+            while !remainder.is_zero() && remainder.degree_plus_one() >= b_degree_plus_1 {
+                let cur_q_coeff = remainder.lead() * divisor_leading_inv;
+                let cur_q_degree = remainder.degree_plus_one() - b_degree_plus_1;
+                quotient.coeffs[cur_q_degree] = cur_q_coeff;
+
+                for (i, &div_coeff) in b.coeffs.iter().enumerate() {
+                    remainder.coeffs[cur_q_degree + i] =
+                        remainder.coeffs[cur_q_degree + i] - (cur_q_coeff * div_coeff);
+                }
+                remainder.trim();
             }
-            denominator_g * root_pow - F::ONE
+            (quotient, remainder)
-        })
+        }
-        .collect::<Vec<_>>();
+    }
-    let denominators_inv = F::batch_multiplicative_inverse(&denominators);
+
-    // Divide every element of `a_eval` by the corresponding denominator.
+    /// Takes a polynomial `a` in coefficient form, and divides it by `Z_H = X^n - 1`.
-    // Then, `a_eval` is the evaluation of `a/Z_H` on `{g.w^i}`.
+    ///
-    a_eval
+    /// This assumes `Z_H | a`, otherwise result is meaningless.
-        .values
+    pub(crate) fn divide_by_z_h(&self, n: usize) -> PolynomialCoeffs<F> {
-        .iter_mut()
+        let mut a = self.clone();
-        .zip(denominators_inv.iter())
+
-        .for_each(|(x, &d)| {
+        // TODO: Is this special case needed?
-            *x *= d;
+        if a.coeffs.iter().all(|p| *p == F::ZERO) {
+            return a;
+        }
+
+        let g = F::MULTIPLICATIVE_GROUP_GENERATOR;
+        let mut g_pow = F::ONE;
+        // Multiply the i-th coefficient of `a` by `g^i`. Then `new_a(w^j) = old_a(g.w^j)`.
+        a.coeffs.iter_mut().for_each(|x| {
+            *x *= g_pow;
+            g_pow *= g;
         });
-    // `p` is the interpolating polynomial of `a_eval` on `{w^i}`.
+
-    let mut p = ifft(a_eval);
+        let root = F::primitive_root_of_unity(log2_strict(a.len()));
-    // We need to scale it by `g^(-i)` to get the interpolating polynomial of `a_eval` on `{g.w^i}`,
+        // Equals to the evaluation of `a` on `{g.w^i}`.
-    // a.k.a `a/Z_H`.
+        let mut a_eval = fft(a);
-    let g_inv = g.inverse();
+        // Compute the denominators `1/(g^n.w^(n*i) - 1)` using batch inversion.
-    let mut g_inv_pow = F::ONE;
+        let denominator_g = g.exp_usize(n);
-    p.coeffs.iter_mut().for_each(|x| {
+        let root_n = root.exp_usize(n);
-        *x *= g_inv_pow;
+        let mut root_pow = F::ONE;
-        g_inv_pow *= g_inv;
+        let denominators = (0..a_eval.len())
-    });
+            .map(|i| {
-    p
+                if i != 0 {
+                    root_pow = root_pow * root_n;
+                }
+                denominator_g * root_pow - F::ONE
+            })
+            .collect::<Vec<_>>();
+        let denominators_inv = F::batch_multiplicative_inverse(&denominators);
+        // Divide every element of `a_eval` by the corresponding denominator.
+        // Then, `a_eval` is the evaluation of `a/Z_H` on `{g.w^i}`.
+        a_eval
+            .values
+            .iter_mut()
+            .zip(denominators_inv.iter())
+            .for_each(|(x, &d)| {
+                *x *= d;
+            });
+        // `p` is the interpolating polynomial of `a_eval` on `{w^i}`.
+        let mut p = ifft(a_eval);
+        // We need to scale it by `g^(-i)` to get the interpolating polynomial of `a_eval` on `{g.w^i}`,
+        // a.k.a `a/Z_H`.
+        let g_inv = g.inverse();
+        let mut g_inv_pow = F::ONE;
+        p.coeffs.iter_mut().for_each(|x| {
+            *x *= g_inv_pow;
+            g_inv_pow *= g_inv;
+        });
+        p
+    }
 }
 
 #[cfg(test)]
 mod tests {
     use crate::field::crandall_field::CrandallField;
     use crate::field::field::Field;
-    use crate::polynomial::division::divide_by_z_h;
     use crate::polynomial::polynomial::PolynomialCoeffs;
 
     #[test]
     fn zero_div_z_h() {
         type F = CrandallField;
         let zero = PolynomialCoeffs::<F>::zero(16);
-        let quotient = divide_by_z_h(zero.clone(), 4);
+        let quotient = zero.divide_by_z_h(4);
         assert_eq!(quotient, zero);
     }
 
@@ -91,7 +158,7 @@ mod tests {
         let a = PolynomialCoeffs::new(vec![-six, -five, -four, -three, six, five, four, three]);
         let q = PolynomialCoeffs::new(vec![six, five, four, three, zero, zero, zero, zero]);
 
-        let computed_q = divide_by_z_h(a, 4);
+        let computed_q = a.divide_by_z_h(4);
         assert_eq!(computed_q, q);
     }
 }

src/polynomial/mod.rs

@@ -1,4 +1,3 @@
 pub mod commitment;
 pub(crate) mod division;
-mod old_polynomial;
 pub mod polynomial;

src/polynomial/old_polynomial.rs

@@ -1,580 +0,0 @@
-#![allow(clippy::many_single_char_names)]
-use crate::field::fft::{
-    fft_precompute, fft_with_precomputation_power_of_2, ifft_with_precomputation_power_of_2,
-    FftPrecomputation,
-};
-use crate::field::field::Field;
-use crate::polynomial::polynomial::PolynomialCoeffs;
-use crate::util::log2_ceil;
-use std::cmp::Ordering;
-use std::ops::{Index, IndexMut, RangeBounds};
-use std::slice::{Iter, IterMut, SliceIndex};
-
-/// Polynomial struct holding a polynomial in coefficient form.
-// TODO: Finish merging this with `PolynomialCoeffs`.
-#[deprecated]
-#[derive(Debug, Clone)]
-pub struct Polynomial<F: Field>(Vec<F>);
-
-impl<F: Field> PartialEq for Polynomial<F> {
-    fn eq(&self, other: &Self) -> bool {
-        let max_terms = self.0.len().max(other.0.len());
-        for i in 0..max_terms {
-            let self_i = self.0.get(i).cloned().unwrap_or(F::ZERO);
-            let other_i = other.0.get(i).cloned().unwrap_or(F::ZERO);
-            if self_i != other_i {
-                return false;
-            }
-        }
-        true
-    }
-}
-
-impl<F: Field> Eq for Polynomial<F> {}
-
-impl<F: Field> From<Vec<F>> for Polynomial<F> {
-    /// Takes a vector of coefficients and returns the corresponding polynomial.
-    fn from(coeffs: Vec<F>) -> Self {
-        Self(coeffs)
-    }
-}
-
-impl<F: Field> From<PolynomialCoeffs<F>> for Polynomial<F> {
-    fn from(coeffs: PolynomialCoeffs<F>) -> Self {
-        Self(coeffs.coeffs)
-    }
-}
-impl<F: Field> From<Polynomial<F>> for PolynomialCoeffs<F> {
-    fn from(poly: Polynomial<F>) -> Self {
-        Self::new(poly.0)
-    }
-}
-
-impl<F, I> Index<I> for Polynomial<F>
-where
-    F: Field,
-    I: SliceIndex<[F]>,
-{
-    type Output = I::Output;
-
-    /// Indexing on the coefficients.
-    fn index(&self, index: I) -> &Self::Output {
-        &self.0[index]
-    }
-}
-
-impl<F, I> IndexMut<I> for Polynomial<F>
-where
-    F: Field,
-    I: SliceIndex<[F]>,
-{
-    fn index_mut(&mut self, index: I) -> &mut <Self as Index<I>>::Output {
-        &mut self.0[index]
-    }
-}
-
-impl<F: Field> Polynomial<F> {
-    /// Takes a slice of coefficients and returns the corresponding polynomial.
-    pub fn from_coeffs(coeffs: &[F]) -> Self {
-        Self(coeffs.to_vec())
-    }
-
-    /// Returns the coefficient vector.
-    pub fn coeffs(&self) -> &[F] {
-        &self.0
-    }
-
-    /// Empty polynomial;
-    pub fn empty() -> Self {
-        Self(Vec::new())
-    }
-
-    /// Zero polynomial with length `len`.
-    /// `len = 1` is the standard representation, but sometimes it's useful to set `len > 1`
-    /// to have polynomials with uniform length.
-    pub fn zero(len: usize) -> Self {
-        Self(vec![F::ZERO; len])
-    }
-
-    pub fn iter(&self) -> Iter<F> {
-        self.0.iter()
-    }
-
-    pub fn iter_mut(&mut self) -> IterMut<F> {
-        self.0.iter_mut()
-    }
-
-    pub fn is_zero(&self) -> bool {
-        self.0.iter().all(|x| x.is_zero())
-    }
-
-    /// Number of coefficients held by the polynomial. Is NOT equal to the degree in general.
-    pub fn len(&self) -> usize {
-        self.0.len()
-    }
-
-    /// Degree of the polynomial.
-    /// Panics on zero polynomial.
-    pub fn degree(&self) -> usize {
-        (0usize..self.len())
-            .rev()
-            .find(|&i| self[i].is_nonzero())
-            .expect("Zero polynomial")
-    }
-
-    /// Degree of the polynomial + 1.
-    fn degree_plus_one(&self) -> usize {
-        (0usize..self.len())
-            .rev()
-            .find(|&i| self[i].is_nonzero())
-            .map_or(0, |i| i + 1)
-    }
-
-    pub fn is_empty(&self) -> bool {
-        self.0.is_empty()
-    }
-
-    /// Removes the coefficients in the range `range`.
-    fn drain<R: RangeBounds<usize>>(&mut self, range: R) {
-        self.0.drain(range);
-    }
-
-    /// Evaluates the polynomial at a point `x`.
-    pub fn eval(&self, x: F) -> F {
-        self.iter().rev().fold(F::ZERO, |acc, &c| acc * x + c)
-    }
-
-    /// Evaluates the polynomial at a point `x`, given the list of powers of `x`.
-    /// Assumes that `self.len() == x_pow.len()`.
-    pub fn eval_from_power(&self, x_pow: &[F]) -> F {
-        self.iter()
-            .zip(x_pow)
-            .fold(F::ZERO, |acc, (&c, &p)| acc + c * p)
-    }
-
-    /// Evaluates the polynomial on subgroup of `F^*` with a given FFT precomputation.
-    pub(crate) fn eval_domain(&self, fft_precomputation: &FftPrecomputation<F>) -> Vec<F> {
-        let domain_size = fft_precomputation.size();
-        if self.len() < domain_size {
-            // Need to pad the polynomial to have the same length as the domain.
-            fft_with_precomputation_power_of_2(
-                self.padded(domain_size).coeffs().to_vec().into(),
-                fft_precomputation,
-            )
-            .values
-        } else {
-            fft_with_precomputation_power_of_2(self.coeffs().to_vec().into(), fft_precomputation)
-                .values
-        }
-    }
-
-    /// Computes the interpolating polynomial of a list of `values` on a subgroup of `F^*`.
-    pub(crate) fn from_evaluations(
-        values: &[F],
-        fft_precomputation: &FftPrecomputation<F>,
-    ) -> Self {
-        Self(ifft_with_precomputation_power_of_2(values.to_vec().into(), fft_precomputation).coeffs)
-    }
-
-    /// Leading coefficient.
-    pub fn lead(&self) -> F {
-        self.iter()
-            .rev()
-            .find(|x| x.is_nonzero())
-            .map_or(F::ZERO, |x| *x)
-    }
-
-    /// Reverse the order of the coefficients, not taking into account the leading zero coefficients.
-    fn rev(&self) -> Self {
-        let d = self.degree();
-        Self(self.0[..=d].iter().rev().copied().collect())
-    }
-
-    /// Negates the polynomial's coefficients.
-    pub fn neg(&self) -> Self {
-        Self(self.iter().map(|&x| -x).collect())
-    }
-
-    /// Multiply the polynomial's coefficients by a scalar.
-    pub(crate) fn scalar_mul(&self, c: F) -> Self {
-        Self(self.iter().map(|&x| c * x).collect())
-    }
-
-    /// Removes leading zero coefficients.
-    pub fn trim(&mut self) {
-        self.0.drain(self.degree_plus_one()..);
-    }
-
-    /// Polynomial addition.
-    pub fn add(&self, other: &Self) -> Self {
-        let (mut a, mut b) = (self.clone(), other.clone());
-        match a.len().cmp(&b.len()) {
-            Ordering::Less => a.pad(b.len()),
-            Ordering::Greater => b.pad(a.len()),
-            _ => (),
-        }
-        Self(a.iter().zip(b.iter()).map(|(&x, &y)| x + y).collect())
-    }
-
-    /// Zero-pad the coefficients to have a given length.
-    pub fn pad(&mut self, len: usize) {
-        self.trim();
-        assert!(self.len() <= len);
-        self.0.extend((self.len()..len).map(|_| F::ZERO));
-    }
-
-    /// Returns the zero-padded polynomial.
-    pub fn padded(&self, len: usize) -> Self {
-        let mut a = self.clone();
-        a.pad(len);
-        a
-    }
-
-    /// Polynomial multiplication.
-    pub fn mul(&self, b: &Self) -> Self {
-        if self.is_zero() || b.is_zero() {
-            return Self::zero(1);
-        }
-        let a_deg = self.degree();
-        let b_deg = b.degree();
-        let new_deg = (a_deg + b_deg + 1).next_power_of_two();
-        let a_pad = self.padded(new_deg);
-        let b_pad = b.padded(new_deg);
-
-        let precomputation = fft_precompute(new_deg);
-        let a_evals = fft_with_precomputation_power_of_2(a_pad.0.to_vec().into(), &precomputation);
-        let b_evals = fft_with_precomputation_power_of_2(b_pad.0.to_vec().into(), &precomputation);
-
-        let mul_evals: Vec<F> = a_evals
-            .values
-            .iter()
-            .zip(b_evals.values.iter())
-            .map(|(&pa, &pb)| pa * pb)
-            .collect();
-        ifft_with_precomputation_power_of_2(mul_evals.to_vec().into(), &precomputation)
-            .coeffs
-            .into()
-    }
-
-    /// Polynomial long division.
-    /// Returns `(q,r)` the quotient and remainder of the polynomial division of `a` by `b`.
-    /// Generally slower that the equivalent function `Polynomial::polynomial_division`.
-    pub fn polynomial_long_division(&self, b: &Self) -> (Self, Self) {
-        // Trim b.
-        let mut b = b.clone();
-        b.trim();
-
-        let (a_degree, b_degree) = (self.degree(), b.degree());
-        if self.is_zero() {
-            (Self::zero(1), Self::empty())
-        } else if b.is_zero() {
-            panic!("Division by zero polynomial");
-        } else if a_degree < b_degree {
-            (Self::zero(1), self.clone())
-        } else {
-            // Now we know that self.degree() >= divisor.degree();
-            let mut quotient = Self::zero(a_degree - b_degree + 1);
-            let mut remainder = self.clone();
-            // Can unwrap here because we know self is not zero.
-            let divisor_leading_inv = b.lead().inverse();
-            while !remainder.is_zero() && remainder.degree() >= b_degree {
-                let cur_q_coeff = remainder.lead() * divisor_leading_inv;
-                let cur_q_degree = remainder.degree() - b_degree;
-                quotient[cur_q_degree] = cur_q_coeff;
-
-                for (i, &div_coeff) in b.iter().enumerate() {
-                    remainder[cur_q_degree + i] =
-                        remainder[cur_q_degree + i] - (cur_q_coeff * div_coeff);
-                }
-                remainder.trim();
-            }
-            (quotient, remainder)
-        }
-    }
-
-    /// Computes the inverse of `self` modulo `x^n`.
-    fn inv_mod_xn(&self, n: usize) -> Self {
-        assert!(self[0].is_nonzero(), "Inverse doesn't exist.");
-        let mut h = self.clone();
-        if h.len() < n {
-            h.pad(n);
-        }
-        let mut a = Self::empty();
-        a.0.push(h[0].inverse());
-        for i in 0..log2_ceil(n) {
-            let l = 1 << i;
-            let h0 = h[..l].to_vec().into();
-            let mut h1: Polynomial<F> = h[l..].to_vec().into();
-            let mut c = a.mul(&h0);
-            if l == c.len() {
-                c = Self::zero(1);
-            } else {
-                c.drain(0..l);
-            }
-            h1.trim();
-            let mut tmp = a.mul(&h1);
-            tmp = tmp.add(&c);
-            tmp.iter_mut().for_each(|x| *x = -(*x));
-            tmp.trim();
-            let mut b = a.mul(&tmp);
-            b.trim();
-            if b.len() > l {
-                b.drain(l..);
-            }
-            a.0.extend_from_slice(&b[..]);
-        }
-        a.drain(n..);
-        a
-    }
-
-    /// Polynomial division.
-    /// Returns `(q,r)` the quotient and remainder of the polynomial division of `a` by `b`.
-    /// Algorithm from http://people.csail.mit.edu/madhu/ST12/scribe/lect06.pdf
-    pub fn polynomial_division(&self, b: &Self) -> (Self, Self) {
-        let (a_degree, b_degree) = (self.degree(), b.degree());
-        if self.is_zero() {
-            (Self::zero(1), Self::empty())
-        } else if b.is_zero() {
-            panic!("Division by zero polynomial");
-        } else if a_degree < b_degree {
-            (Self::zero(1), self.clone())
-        } else if b_degree == 0 {
-            (self.scalar_mul(b[0].inverse()), Self::empty())
-        } else {
-            let rev_b = b.rev();
-            let rev_b_inv = rev_b.inv_mod_xn(a_degree - b_degree + 1);
-            let rev_q: Polynomial<F> = rev_b_inv
-                .mul(&self.rev()[..=a_degree - b_degree].to_vec().into())[..=a_degree - b_degree]
-                .to_vec()
-                .into();
-            let mut q = rev_q.rev();
-            let mut qb = q.mul(b);
-            qb.pad(self.len());
-            let mut r = self.add(&qb.neg());
-            q.trim();
-            r.trim();
-            (q, r)
-        }
-    }
-
-    // Divides a polynomial `a` by `Z_H = X^n - 1`. Assumes `Z_H | a`, otherwise result is meaningless.
-    pub fn divide_by_z_h(&self, n: usize) -> Self {
-        if self.is_zero() {
-            return self.clone();
-        }
-        let mut a_trim = self.clone();
-        a_trim.trim();
-        let g = F::MULTIPLICATIVE_GROUP_GENERATOR;
-        let mut g_pow = F::ONE;
-        // Multiply the i-th coefficient of `a` by `g^i`. Then `new_a(w^j) = old_a(g.w^j)`.
-        a_trim.iter_mut().for_each(|x| {
-            *x = (*x) * g_pow;
-            g_pow = g * g_pow;
-        });
-        let d = a_trim.degree();
-        let root = F::primitive_root_of_unity(log2_ceil(d + 1));
-        let precomputation = fft_precompute(d + 1);
-        // Equals to the evaluation of `a` on `{g.w^i}`.
-        let mut a_eval = a_trim.eval_domain(&precomputation);
-        // Compute the denominators `1/(g^n.w^(n*i) - 1)` using batch inversion.
-        let denominator_g = g.exp_usize(n);
-        let root_n = root.exp_usize(n);
-        let mut root_pow = F::ONE;
-        let denominators = (0..a_eval.len())
-            .map(|i| {
-                if i != 0 {
-                    root_pow = root_pow * root_n;
-                }
-                denominator_g * root_pow - F::ONE
-            })
-            .collect::<Vec<_>>();
-        let denominators_inv = F::batch_multiplicative_inverse(&denominators);
-        // Divide every element of `a_eval` by the corresponding denominator.
-        // Then, `a_eval` is the evaluation of `a/Z_H` on `{g.w^i}`.
-        a_eval
-            .iter_mut()
-            .zip(denominators_inv.iter())
-            .for_each(|(x, &d)| {
-                *x = (*x) * d;
-            });
-        // `p` is the interpolating polynomial of `a_eval` on `{w^i}`.
-        let mut p = Self::from_evaluations(&a_eval, &precomputation);
-        // We need to scale it by `g^(-i)` to get the interpolating polynomial of `a_eval` on `{g.w^i}`,
-        // a.k.a `a/Z_H`.
-        let g_inv = g.inverse();
-        let mut g_inv_pow = F::ONE;
-        p.iter_mut().for_each(|x| {
-            *x = (*x) * g_inv_pow;
-            g_inv_pow = g_inv_pow * g_inv;
-        });
-        p
-    }
-}
-
-#[cfg(test)]
-mod test {
-    use super::*;
-    use crate::field::crandall_field::CrandallField;
-    use rand::{thread_rng, Rng};
-    use std::time::Instant;
-
-    #[test]
-    fn test_polynomial_multiplication() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000));
-        let a = Polynomial(F::rand_vec(a_deg));
-        let b = Polynomial(F::rand_vec(b_deg));
-        let m1 = a.mul(&b);
-        let m2 = a.mul(&b);
-        for _ in 0..1000 {
-            let x = F::rand();
-            assert_eq!(m1.eval(x), a.eval(x) * b.eval(x));
-            assert_eq!(m2.eval(x), a.eval(x) * b.eval(x));
-        }
-    }
-
-    #[test]
-    fn test_inv_mod_xn() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let a_deg = rng.gen_range(1, 1_000);
-        let n = rng.gen_range(1, 1_000);
-        let a = Polynomial(F::rand_vec(a_deg));
-        let b = a.inv_mod_xn(n);
-        let mut m = a.mul(&b);
-        m.drain(n..);
-        m.trim();
-        assert_eq!(
-            m,
-            Polynomial(vec![F::ONE]),
-            "a: {:#?}, b:{:#?}, n:{:#?}, m:{:#?}",
-            a,
-            b,
-            n,
-            m
-        );
-    }
-
-    #[test]
-    fn test_polynomial_long_division() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000));
-        let a = Polynomial(F::rand_vec(a_deg));
-        let b = Polynomial(F::rand_vec(b_deg));
-        let (q, r) = a.polynomial_long_division(&b);
-        for _ in 0..1000 {
-            let x = F::rand();
-            assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x));
-        }
-    }
-
-    #[test]
-    fn test_polynomial_division() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000));
-        let a = Polynomial(F::rand_vec(a_deg));
-        let b = Polynomial(F::rand_vec(b_deg));
-        let (q, r) = a.polynomial_division(&b);
-        for _ in 0..1000 {
-            let x = F::rand();
-            assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x));
-        }
-    }
-
-    #[test]
-    fn test_polynomial_division_by_constant() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let a_deg = rng.gen_range(1, 10_000);
-        let a = Polynomial(F::rand_vec(a_deg));
-        let b = Polynomial::from(vec![F::rand()]);
-        let (q, r) = a.polynomial_division(&b);
-        for _ in 0..1000 {
-            let x = F::rand();
-            assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x));
-        }
-    }
-
-    #[test]
-    fn test_division_by_z_h() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let a_deg = rng.gen_range(1, 10_000);
-        let n = rng.gen_range(1, a_deg);
-        let mut a = Polynomial(F::rand_vec(a_deg));
-        a.trim();
-        let z_h = {
-            let mut z_h_vec = vec![F::ZERO; n + 1];
-            z_h_vec[n] = F::ONE;
-            z_h_vec[0] = F::NEG_ONE;
-            Polynomial(z_h_vec)
-        };
-        let m = a.mul(&z_h);
-        let now = Instant::now();
-        let mut a_test = m.divide_by_z_h(n);
-        a_test.trim();
-        println!("Division time: {:?}", now.elapsed());
-        assert_eq!(a, a_test);
-    }
-
-    #[test]
-    fn divide_zero_poly_by_z_h() {
-        let zero_poly = Polynomial::<CrandallField>::empty();
-        zero_poly.divide_by_z_h(16);
-    }
-
-    // Test to see which polynomial division method is faster for divisions of the type
-    // `(X^n - 1)/(X - a)
-    #[test]
-    fn test_division_linear() {
-        type F = CrandallField;
-        let mut rng = thread_rng();
-        let l = 14;
-        let n = 1 << l;
-        let g = F::primitive_root_of_unity(l);
-        let xn_minus_one = {
-            let mut xn_min_one_vec = vec![F::ZERO; n + 1];
-            xn_min_one_vec[n] = F::ONE;
-            xn_min_one_vec[0] = F::NEG_ONE;
-            Polynomial(xn_min_one_vec)
-        };
-
-        let a = g.exp_usize(rng.gen_range(0, n));
-        let denom = Polynomial(vec![-a, F::ONE]);
-        let now = Instant::now();
-        xn_minus_one.polynomial_division(&denom);
-        println!("Division time: {:?}", now.elapsed());
-        let now = Instant::now();
-        xn_minus_one.polynomial_long_division(&denom);
-        println!("Division time: {:?}", now.elapsed());
-    }
-
-    #[test]
-    fn eq() {
-        type F = CrandallField;
-        assert_eq!(Polynomial::<F>(vec![]), Polynomial(vec![]));
-        assert_eq!(Polynomial::<F>(vec![F::ZERO]), Polynomial(vec![F::ZERO]));
-        assert_eq!(Polynomial::<F>(vec![]), Polynomial(vec![F::ZERO]));
-        assert_eq!(Polynomial::<F>(vec![F::ZERO]), Polynomial(vec![]));
-        assert_eq!(
-            Polynomial::<F>(vec![F::ZERO]),
-            Polynomial(vec![F::ZERO, F::ZERO])
-        );
-        assert_eq!(
-            Polynomial::<F>(vec![F::ONE]),
-            Polynomial(vec![F::ONE, F::ZERO])
-        );
-        assert_ne!(Polynomial::<F>(vec![]), Polynomial(vec![F::ONE]));
-        assert_ne!(
-            Polynomial::<F>(vec![F::ZERO]),
-            Polynomial(vec![F::ZERO, F::ONE])
-        );
-        assert_ne!(
-            Polynomial::<F>(vec![F::ZERO]),
-            Polynomial(vec![F::ONE, F::ZERO])
-        );
-    }
-}

src/polynomial/polynomial.rs

@@ -3,8 +3,7 @@ use std::ops::{Add, Mul, Sub};
 
 use crate::field::fft::{fft, ifft};
 use crate::field::field::Field;
-use crate::polynomial::old_polynomial::Polynomial;
+use crate::util::{log2_ceil, log2_strict};
-use crate::util::log2_strict;
 
 /// A polynomial in point-value form.
 ///
@@ -49,7 +48,7 @@ impl<F: Field> From<Vec<F>> for PolynomialValues<F> {
 }
 
 /// A polynomial in coefficient form.
-#[derive(Clone, Debug, Eq, PartialEq)]
+#[derive(Clone, Debug)]
 pub struct PolynomialCoeffs<F: Field> {
     pub(crate) coeffs: Vec<F>,
 }
@@ -67,6 +66,10 @@ impl<F: Field> PolynomialCoeffs<F> {
         PolynomialCoeffs { coeffs }
     }
 
+    pub(crate) fn empty() -> Self {
+        Self::new(Vec::new())
+    }
+
     pub(crate) fn zero(len: usize) -> Self {
         Self::new(vec![F::ZERO; len])
     }
@@ -112,6 +115,7 @@ impl<F: Field> PolynomialCoeffs<F> {
     }
 
     pub(crate) fn padded(&self, new_len: usize) -> Self {
+        assert!(new_len >= self.len());
         let mut coeffs = self.coeffs.clone();
         coeffs.resize(new_len, F::ZERO);
         Self { coeffs }
@@ -119,10 +123,16 @@ impl<F: Field> PolynomialCoeffs<F> {
 
     /// Removes leading zero coefficients.
     pub fn trim(&mut self) {
-        self.coeffs.drain(self.degree_plus_one()..);
+        self.coeffs.truncate(self.degree_plus_one());
     }
 
-    /// Degree of the polynomial + 1.
+    /// Removes leading zero coefficients.
+    pub fn trimmed(&self) -> Self {
+        let coeffs = self.coeffs[..self.degree_plus_one()].to_vec();
+        Self { coeffs }
+    }
+
+    /// Degree of the polynomial + 1, or 0 for a polynomial with no non-zero coefficients.
     pub(crate) fn degree_plus_one(&self) -> usize {
         (0usize..self.len())
             .rev()
@@ -130,11 +140,50 @@ impl<F: Field> PolynomialCoeffs<F> {
             .map_or(0, |i| i + 1)
     }
 
-    pub(crate) fn div_rem(&self, rhs: &Self) -> (Self, Self) {
+    /// Leading coefficient.
-        let lhs = Polynomial::from(self.clone());
+    pub fn lead(&self) -> F {
-        let rhs = Polynomial::from(rhs.clone());
+        self.coeffs
-        let (q, r) = lhs.polynomial_division(&rhs);
+            .iter()
-        (q.into(), r.into())
+            .rev()
+            .find(|x| x.is_nonzero())
+            .map_or(F::ZERO, |x| *x)
+    }
+
+    /// Computes the inverse of `self` modulo `x^n`.
+    pub(crate) fn inv_mod_xn(&self, n: usize) -> Self {
+        assert!(self.coeffs[0].is_nonzero(), "Inverse doesn't exist.");
+        let mut h = self.padded(n);
+        let mut a = Self::empty();
+        a.coeffs.push(h.coeffs[0].inverse());
+        for i in 0..log2_ceil(n) {
+            let l = 1 << i;
+            let h0 = h.coeffs[..l].to_vec().into();
+            let mut h1: Self = h.coeffs[l..].to_vec().into();
+            let mut c = &a * &h0;
+            if l == c.len() {
+                c = Self::zero(1);
+            } else {
+                c.coeffs.drain(0..l);
+            }
+            h1.trim();
+            let mut tmp = a.mul(&h1);
+            tmp = tmp.add(&c);
+            tmp.coeffs.iter_mut().for_each(|x| *x = -(*x));
+            tmp.trim();
+            let mut b = &a * &tmp;
+            b.trim();
+            if b.len() > l {
+                b.coeffs.drain(l..);
+            }
+            a.coeffs.extend_from_slice(&b.coeffs);
+        }
+        a.coeffs.drain(n..);
+        a
+    }
+
+    /// Reverse the order of the coefficients, not taking into account the leading zero coefficients.
+    pub(crate) fn rev(&self) -> Self {
+        Self::new(self.trimmed().coeffs.into_iter().rev().collect())
     }
 
     pub fn fft(self) -> PolynomialValues<F> {
@@ -152,6 +201,22 @@ impl<F: Field> PolynomialCoeffs<F> {
     }
 }
 
+impl<F: Field> PartialEq for PolynomialCoeffs<F> {
+    fn eq(&self, other: &Self) -> bool {
+        let max_terms = self.coeffs.len().max(other.coeffs.len());
+        for i in 0..max_terms {
+            let self_i = self.coeffs.get(i).cloned().unwrap_or(F::ZERO);
+            let other_i = other.coeffs.get(i).cloned().unwrap_or(F::ZERO);
+            if self_i != other_i {
+                return false;
+            }
+        }
+        true
+    }
+}
+
+impl<F: Field> Eq for PolynomialCoeffs<F> {}
+
 impl<F: Field> From<Vec<F>> for PolynomialCoeffs<F> {
     fn from(coeffs: Vec<F>) -> Self {
         Self::new(coeffs)
@@ -215,10 +280,40 @@ impl<F: Field> Mul for &PolynomialCoeffs<F> {
 
 #[cfg(test)]
 mod tests {
+    use std::time::Instant;
+
+    use rand::{thread_rng, Rng};
+
     use crate::field::crandall_field::CrandallField;
 
     use super::*;
 
+    #[test]
+    fn test_trimmed() {
+        type F = CrandallField;
+
+        assert_eq!(
+            PolynomialCoeffs::<F> { coeffs: vec![] }.trimmed(),
+            PolynomialCoeffs::<F> { coeffs: vec![] }
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F> {
+                coeffs: vec![F::ZERO]
+            }
+            .trimmed(),
+            PolynomialCoeffs::<F> { coeffs: vec![] }
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F> {
+                coeffs: vec![F::ONE, F::TWO, F::ZERO, F::ZERO]
+            }
+            .trimmed(),
+            PolynomialCoeffs::<F> {
+                coeffs: vec![F::ONE, F::TWO]
+            }
+        );
+    }
+
     #[test]
     fn test_coset_fft() {
         type F = CrandallField;
@@ -237,4 +332,179 @@ mod tests {
 
         assert_eq!(coset_evals, naive_coset_evals);
     }
+
+    #[test]
+    fn test_polynomial_multiplication() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000));
+        let a = PolynomialCoeffs::new(F::rand_vec(a_deg));
+        let b = PolynomialCoeffs::new(F::rand_vec(b_deg));
+        let m1 = &a * &b;
+        let m2 = &a * &b;
+        for _ in 0..1000 {
+            let x = F::rand();
+            assert_eq!(m1.eval(x), a.eval(x) * b.eval(x));
+            assert_eq!(m2.eval(x), a.eval(x) * b.eval(x));
+        }
+    }
+
+    #[test]
+    fn test_inv_mod_xn() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let a_deg = rng.gen_range(1, 1_000);
+        let n = rng.gen_range(1, 1_000);
+        let a = PolynomialCoeffs::new(F::rand_vec(a_deg));
+        let b = a.inv_mod_xn(n);
+        let mut m = &a * &b;
+        m.coeffs.drain(n..);
+        m.trim();
+        assert_eq!(
+            m,
+            PolynomialCoeffs::new(vec![F::ONE]),
+            "a: {:#?}, b:{:#?}, n:{:#?}, m:{:#?}",
+            a,
+            b,
+            n,
+            m
+        );
+    }
+
+    #[test]
+    fn test_polynomial_long_division() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000));
+        let a = PolynomialCoeffs::new(F::rand_vec(a_deg));
+        let b = PolynomialCoeffs::new(F::rand_vec(b_deg));
+        let (q, r) = a.div_rem_long_division(&b);
+        for _ in 0..1000 {
+            let x = F::rand();
+            assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x));
+        }
+    }
+
+    #[test]
+    fn test_polynomial_division() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000));
+        let a = PolynomialCoeffs::new(F::rand_vec(a_deg));
+        let b = PolynomialCoeffs::new(F::rand_vec(b_deg));
+        let (q, r) = a.div_rem(&b);
+        for _ in 0..1000 {
+            let x = F::rand();
+            assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x));
+        }
+    }
+
+    #[test]
+    fn test_polynomial_division_by_constant() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let a_deg = rng.gen_range(1, 10_000);
+        let a = PolynomialCoeffs::new(F::rand_vec(a_deg));
+        let b = PolynomialCoeffs::from(vec![F::rand()]);
+        let (q, r) = a.div_rem(&b);
+        for _ in 0..1000 {
+            let x = F::rand();
+            assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x));
+        }
+    }
+
+    #[test]
+    fn test_division_by_z_h() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let a_deg = rng.gen_range(1, 10_000);
+        let n = rng.gen_range(1, a_deg);
+        let mut a = PolynomialCoeffs::new(F::rand_vec(a_deg));
+        a.trim();
+        let z_h = {
+            let mut z_h_vec = vec![F::ZERO; n + 1];
+            z_h_vec[n] = F::ONE;
+            z_h_vec[0] = F::NEG_ONE;
+            PolynomialCoeffs::new(z_h_vec)
+        };
+        let m = &a * &z_h;
+        let now = Instant::now();
+        let mut a_test = m.divide_by_z_h(n);
+        a_test.trim();
+        println!("Division time: {:?}", now.elapsed());
+        assert_eq!(a, a_test);
+    }
+
+    #[test]
+    fn divide_zero_poly_by_z_h() {
+        let zero_poly = PolynomialCoeffs::<CrandallField>::empty();
+        zero_poly.divide_by_z_h(16);
+    }
+
+    // Test to see which polynomial division method is faster for divisions of the type
+    // `(X^n - 1)/(X - a)
+    #[test]
+    fn test_division_linear() {
+        type F = CrandallField;
+        let mut rng = thread_rng();
+        let l = 14;
+        let n = 1 << l;
+        let g = F::primitive_root_of_unity(l);
+        let xn_minus_one = {
+            let mut xn_min_one_vec = vec![F::ZERO; n + 1];
+            xn_min_one_vec[n] = F::ONE;
+            xn_min_one_vec[0] = F::NEG_ONE;
+            PolynomialCoeffs::new(xn_min_one_vec)
+        };
+
+        let a = g.exp_usize(rng.gen_range(0, n));
+        let denom = PolynomialCoeffs::new(vec![-a, F::ONE]);
+        let now = Instant::now();
+        xn_minus_one.div_rem(&denom);
+        println!("Division time: {:?}", now.elapsed());
+        let now = Instant::now();
+        xn_minus_one.div_rem_long_division(&denom);
+        println!("Division time: {:?}", now.elapsed());
+    }
+
+    #[test]
+    fn eq() {
+        type F = CrandallField;
+        assert_eq!(
+            PolynomialCoeffs::<F>::new(vec![]),
+            PolynomialCoeffs::new(vec![])
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F>::new(vec![F::ZERO]),
+            PolynomialCoeffs::new(vec![F::ZERO])
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F>::new(vec![]),
+            PolynomialCoeffs::new(vec![F::ZERO])
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F>::new(vec![F::ZERO]),
+            PolynomialCoeffs::new(vec![])
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F>::new(vec![F::ZERO]),
+            PolynomialCoeffs::new(vec![F::ZERO, F::ZERO])
+        );
+        assert_eq!(
+            PolynomialCoeffs::<F>::new(vec![F::ONE]),
+            PolynomialCoeffs::new(vec![F::ONE, F::ZERO])
+        );
+        assert_ne!(
+            PolynomialCoeffs::<F>::new(vec![]),
+            PolynomialCoeffs::new(vec![F::ONE])
+        );
+        assert_ne!(
+            PolynomialCoeffs::<F>::new(vec![F::ZERO]),
+            PolynomialCoeffs::new(vec![F::ZERO, F::ONE])
+        );
+        assert_ne!(
+            PolynomialCoeffs::<F>::new(vec![F::ZERO]),
+            PolynomialCoeffs::new(vec![F::ONE, F::ZERO])
+        );
+    }
 }

src/prover.rs

@@ -10,7 +10,6 @@ use crate::generator::generate_partial_witness;
 use crate::plonk_challenger::Challenger;
 use crate::plonk_common::{eval_l_1, evaluate_gate_constraints, reduce_with_powers_multi};
 use crate::polynomial::commitment::ListPolynomialCommitment;
-use crate::polynomial::division::divide_by_z_h;
 use crate::polynomial::polynomial::{PolynomialCoeffs, PolynomialValues};
 use crate::proof::Proof;
 use crate::util::transpose;
@@ -106,7 +105,7 @@ pub(crate) fn prove<F: Field>(
             let mut all_quotient_poly_chunks = Vec::with_capacity(num_checks * quotient_degree);
             for vanishing_poly in vanishing_polys.into_iter() {
                 let vanishing_poly_coeff = ifft(vanishing_poly);
-                let quotient_poly_coeff = divide_by_z_h(vanishing_poly_coeff, degree);
+                let quotient_poly_coeff = vanishing_poly_coeff.divide_by_z_h(degree);
                 // Split t into degree-n chunks.
                 let quotient_poly_coeff_chunks = quotient_poly_coeff.chunks(degree);
                 all_quotient_poly_chunks.extend(quotient_poly_coeff_chunks);