Optimize some polynomial operations

src/field/lagrange.rs

@@ -65,11 +65,35 @@ pub fn barycentric_weights<F: Field>(points: &[(F, F)]) -> Vec<F> {
     )
 }
 
+/// Interpolate the linear polynomial passing through `points` on `x`.
+pub fn interpolate2<F: Field>(points: [(F, F); 2], x: F) -> F {
+    // a0 -> a1
+    // b0 -> b1
+    // x  -> a1 + (x-a0)*(b1-a1)/(b0-a0)
+    let (a0, a1) = points[0];
+    let (b0, b1) = points[1];
+    assert_ne!(a0, b0);
+    a1 + (x - a0) * (b1 - a1) / (b0 - a0)
+}
+
+/// Returns the linear polynomial passing through `points`.
+pub fn interpolant2<F: Field>(points: [(F, F); 2]) -> PolynomialCoeffs<F> {
+    // a0 -> a1
+    // b0 -> b1
+    // x  -> a1 + (x-a0)*(b1-a1)/(b0-a0)
+    let (a0, a1) = points[0];
+    let (b0, b1) = points[1];
+    assert_ne!(a0, b0);
+    let mult = (b1 - a1) / (b0 - a0);
+    vec![a1 - a0 * mult, mult].into()
+}
+
 #[cfg(test)]
 mod tests {
+    use super::*;
     use crate::field::crandall_field::CrandallField;
+    use crate::field::extension_field::quartic::QuarticCrandallField;
     use crate::field::field::Field;
-    use crate::field::lagrange::interpolant;
     use crate::polynomial::polynomial::PolynomialCoeffs;
 
     #[test]
@@ -120,4 +144,19 @@ mod tests {
     fn eval_naive<F: Field>(coeffs: &PolynomialCoeffs<F>, domain: &[F]) -> Vec<(F, F)> {
         domain.iter().map(|&x| (x, coeffs.eval(x))).collect()
     }
+
+    #[test]
+    fn test_interpolant2() {
+        type F = QuarticCrandallField;
+        let points = [(F::rand(), F::rand()), (F::rand(), F::rand())];
+        let x = F::rand();
+
+        let intepol0 = interpolant(&points);
+        let intepol1 = interpolant2(points);
+        assert_eq!(intepol0.trimmed(), intepol1.trimmed());
+
+        let ev0 = interpolate(&points, x, &barycentric_weights(&points));
+        let ev1 = interpolate2(points, x);
+        assert_eq!(ev0, ev1);
+    }
 }

src/polynomial/commitment.rs

@@ -4,7 +4,7 @@ use rayon::prelude::*;
 use crate::field::extension_field::Extendable;
 use crate::field::extension_field::{FieldExtension, Frobenius};
 use crate::field::field::Field;
-use crate::field::lagrange::interpolant;
+use crate::field::lagrange::{interpolant, interpolant2};
 use crate::fri::{prover::fri_proof, verifier::verify_fri_proof, FriConfig};
 use crate::merkle_tree::MerkleTree;
 use crate::plonk_challenger::Challenger;
@@ -122,15 +122,10 @@ impl<F: Field> ListPolynomialCommitment<F> {
             .map(|p| p.to_extension());
         let single_os = [&os.constants, &os.plonk_s_sigmas, &os.quotient_polys];
         let single_evals = single_os.iter().flat_map(|v| v.iter());
-        let single_composition_poly = reduce_polys_with_iter(single_polys, alpha_powers.clone());
+        let single_composition_poly = reduce_polys_with_iter(single_polys, &mut alpha_powers);
-        let single_composition_eval = reduce_with_iter(single_evals, &mut alpha_powers);
 
-        let single_quotient = Self::compute_quotient(
+        let single_quotient = Self::compute_quotient1(zeta, single_composition_poly);
-            &[zeta],
+        final_poly += single_quotient;
-            &[single_composition_eval],
-            &single_composition_poly,
-        );
-        final_poly = &final_poly + &single_quotient;
 
         // Zs polynomials are opened at `zeta` and `g*zeta`.
         let zs_polys = commitments[3].polynomials.iter().map(|p| p.to_extension());
@@ -140,12 +135,9 @@ impl<F: Field> ListPolynomialCommitment<F> {
             reduce_with_iter(&os.plonk_zs_right, &mut alpha_powers),
         ];
 
-        let zs_quotient = Self::compute_quotient(
+        let zs_quotient =
-            &[zeta, g * zeta],
+            Self::compute_quotient2([zeta, g * zeta], zs_composition_evals, zs_composition_poly);
-            &zs_composition_evals,
+        final_poly += zs_quotient;
-            &zs_composition_poly,
-        );
-        final_poly = &final_poly + &zs_quotient;
 
         // When working in an extension field, need to check that wires are in the base field.
         // Check this by opening the wires polynomials at `zeta` and `zeta.frobenius()` and using the fact that
@@ -158,12 +150,12 @@ impl<F: Field> ListPolynomialCommitment<F> {
             reduce_with_iter(&wire_evals_frob, alpha_powers),
         ];
 
-        let wires_quotient = Self::compute_quotient(
+        let wires_quotient = Self::compute_quotient2(
-            &[zeta, zeta.frobenius()],
+            [zeta, zeta.frobenius()],
-            &wire_composition_evals,
+            wire_composition_evals,
-            &wire_composition_poly,
+            wire_composition_poly,
         );
-        final_poly = &final_poly + &wires_quotient;
+        final_poly += wires_quotient;
 
         let lde_final_poly = final_poly.lde(config.rate_bits);
         let lde_final_values = lde_final_poly
@@ -192,28 +184,41 @@ impl<F: Field> ListPolynomialCommitment<F> {
         )
     }
 
-    /// Given `points=(x_i)`, `evals=(y_i)` and `poly=P` with `P(x_i)=y_i`, computes the polynomial
+    /// Given `x` and `poly=P(X)`, computes the polynomial `Q=(P-P(x))/(X-x)`.
-    /// `Q=(P-I)/Z` where `I` interpolates `(x_i, y_i)` and `Z` is the vanishing polynomial on `(x_i)`.
+    fn compute_quotient1<const D: usize>(
-    fn compute_quotient<const D: usize>(
+        point: F::Extension,
-        points: &[F::Extension],
+        poly: PolynomialCoeffs<F::Extension>,
-        evals: &[F::Extension],
-        poly: &PolynomialCoeffs<F::Extension>,
     ) -> PolynomialCoeffs<F::Extension>
     where
         F: Extendable<D>,
     {
-        let pairs = points
+        let (quotient, _ev) = poly.divide_by_linear(point);
-            .iter()
+        quotient.padded(quotient.degree_plus_one().next_power_of_two())
-            .zip(evals)
+    }
-            .map(|(&x, &e)| (x, e))
+
-            .collect::<Vec<_>>();
+    /// Given `points=(x_i)`, `evals=(y_i)` and `poly=P` with `P(x_i)=y_i`, computes the polynomial
+    /// `Q=(P-I)/Z` where `I` interpolates `(x_i, y_i)` and `Z` is the vanishing polynomial on `(x_i)`.
+    fn compute_quotient2<const D: usize>(
+        points: [F::Extension; 2],
+        evals: [F::Extension; 2],
+        poly: PolynomialCoeffs<F::Extension>,
+    ) -> PolynomialCoeffs<F::Extension>
+    where
+        F: Extendable<D>,
+    {
+        let pairs = [(points[0], evals[0]), (points[1], evals[1])];
         debug_assert!(pairs.iter().all(|&(x, e)| poly.eval(x) == e));
 
-        let interpolant = interpolant(&pairs);
+        let interpolant = interpolant2(pairs);
-        let denominator = points.iter().fold(PolynomialCoeffs::one(), |acc, &x| {
+        let denominator = vec![
-            &acc * &PolynomialCoeffs::new(vec![-x, F::Extension::ONE])
+            points[0] * points[1],
-        });
+            -points[0] - points[1],
-        let numerator = poly - &interpolant;
+            F::Extension::ONE,
+        ]
+        .into();
+
+        let mut numerator = poly;
+        numerator -= interpolant;
         let (quotient, rem) = numerator.div_rem(&denominator);
         debug_assert!(rem.is_zero());
 

src/polynomial/division.rs

@@ -125,8 +125,25 @@ impl<F: Field> PolynomialCoeffs<F> {
         p
     }
 
+    /// Let `self=p(X)`, this returns `(p(X)-p(z))/(X-z)` and `p(z)`.
+    /// See https://en.wikipedia.org/wiki/Horner%27s_method
+    pub(crate) fn divide_by_linear(&self, z: F) -> (PolynomialCoeffs<F>, F) {
+        let mut bs = self
+            .coeffs
+            .iter()
+            .rev()
+            .scan(F::ZERO, |acc, &c| {
+                *acc = *acc * z + c;
+                Some(*acc)
+            })
+            .collect::<Vec<_>>();
+        let ev = bs.pop().unwrap_or(F::ZERO);
+        bs.reverse();
+        (Self { coeffs: bs }, ev)
+    }
+
     /// Computes the inverse of `self` modulo `x^n`.
-    pub(crate) fn inv_mod_xn(&self, n: usize) -> Self {
+    pub fn inv_mod_xn(&self, n: usize) -> Self {
         assert!(self.coeffs[0].is_nonzero(), "Inverse doesn't exist.");
 
         let h = if self.len() < n {
@@ -166,7 +183,10 @@ impl<F: Field> PolynomialCoeffs<F> {
 
 #[cfg(test)]
 mod tests {
+    use std::time::Instant;
+
     use crate::field::crandall_field::CrandallField;
+    use crate::field::extension_field::quartic::QuarticCrandallField;
     use crate::field::field::Field;
     use crate::polynomial::polynomial::PolynomialCoeffs;
 
@@ -199,4 +219,49 @@ mod tests {
         let computed_q = a.divide_by_z_h(4);
         assert_eq!(computed_q, q);
     }
+
+    #[test]
+    #[ignore]
+    fn test_division_by_linear() {
+        type F = QuarticCrandallField;
+        let n = 1_000_000;
+        let poly = PolynomialCoeffs::new(F::rand_vec(n));
+        let z = F::rand();
+        let ev = poly.eval(z);
+
+        let timer = Instant::now();
+        let (quotient, ev2) = poly.div_rem(&PolynomialCoeffs::new(vec![-z, F::ONE]));
+        println!("{:.3}s for usual", timer.elapsed().as_secs_f32());
+        assert_eq!(ev2.trimmed().coeffs, vec![ev]);
+
+        let timer = Instant::now();
+        let (quotient, ev3) = poly.div_rem_long_division(&PolynomialCoeffs::new(vec![-z, F::ONE]));
+        println!("{:.3}s for long division", timer.elapsed().as_secs_f32());
+        assert_eq!(ev3.trimmed().coeffs, vec![ev]);
+
+        let timer = Instant::now();
+        let horn = poly.divide_by_linear(z);
+        println!("{:.3}s for Horner", timer.elapsed().as_secs_f32());
+        assert_eq!((quotient, ev), horn);
+    }
+
+    #[test]
+    #[ignore]
+    fn test_division_by_quadratic() {
+        type F = QuarticCrandallField;
+        let n = 1_000_000;
+        let poly = PolynomialCoeffs::new(F::rand_vec(n));
+        let quad = PolynomialCoeffs::new(F::rand_vec(2));
+
+        let timer = Instant::now();
+        let (quotient0, rem0) = poly.div_rem(&quad);
+        println!("{:.3}s for usual", timer.elapsed().as_secs_f32());
+
+        let timer = Instant::now();
+        let (quotient1, rem1) = poly.div_rem_long_division(&quad);
+        println!("{:.3}s for long division", timer.elapsed().as_secs_f32());
+
+        assert_eq!(quotient0.trimmed(), quotient1.trimmed());
+        assert_eq!(rem0.trimmed(), rem1.trimmed());
+    }
 }

src/polynomial/polynomial.rs

@@ -1,6 +1,6 @@
 use std::cmp::max;
 use std::iter::Sum;
-use std::ops::{Add, Mul, Sub};
+use std::ops::{Add, AddAssign, Mul, Sub, SubAssign};
 
 use crate::field::extension_field::Extendable;
 use crate::field::fft::{fft, ifft};
@@ -243,6 +243,26 @@ impl<F: Field> Sub for &PolynomialCoeffs<F> {
     }
 }
 
+impl<F: Field> AddAssign for PolynomialCoeffs<F> {
+    fn add_assign(&mut self, rhs: Self) {
+        let len = max(self.len(), rhs.len());
+        self.coeffs.resize(len, F::ZERO);
+        for (l, r) in self.coeffs.iter_mut().zip(rhs.coeffs) {
+            *l += r;
+        }
+    }
+}
+
+impl<F: Field> SubAssign for PolynomialCoeffs<F> {
+    fn sub_assign(&mut self, rhs: Self) {
+        let len = max(self.len(), rhs.len());
+        self.coeffs.resize(len, F::ZERO);
+        for (l, r) in self.coeffs.iter_mut().zip(rhs.coeffs) {
+            *l -= r;
+        }
+    }
+}
+
 impl<F: Field> Mul<F> for &PolynomialCoeffs<F> {
     type Output = PolynomialCoeffs<F>;