Fix bug

src/fri/verifier.rs

@@ -2,7 +2,7 @@ use anyhow::{ensure, Result};
 
 use crate::field::extension_field::{flatten, Extendable, FieldExtension, Frobenius};
 use crate::field::field::Field;
-use crate::field::lagrange::{barycentric_weights, interpolant, interpolate};
+use crate::field::interpolation::{barycentric_weights, interpolate, interpolate2};
 use crate::fri::FriConfig;
 use crate::hash::hash_n_to_1;
 use crate::merkle_proofs::verify_merkle_proof;
@@ -192,11 +192,14 @@ fn fri_combine_initial<F: Field + Extendable<D>, const D: usize>(
         .map(|&e| F::Extension::from_basefield(e));
     let zs_composition_eval = alpha.clone().reduce(zs_evals);
     let zeta_right = F::Extension::primitive_root_of_unity(degree_log) * zeta;
-    let zs_interpol = interpolant(&[
-        (zeta, alpha.clone().reduce(os.plonk_zs.iter())),
-        (zeta_right, alpha.reduce(os.plonk_zs_right.iter())),
-    ]);
-    let zs_numerator = zs_composition_eval - zs_interpol.eval(subgroup_x);
+    let zs_interpol = interpolate2(
+        [
+            (zeta, alpha.clone().reduce(os.plonk_zs.iter())),
+            (zeta_right, alpha.reduce(os.plonk_zs_right.iter())),
+        ],
+        subgroup_x,
+    );
+    let zs_numerator = zs_composition_eval - zs_interpol;
     let zs_denominator = (subgroup_x - zeta) * (subgroup_x - zeta_right);
     sum = alpha.shift(sum);
     sum += zs_numerator / zs_denominator;
@@ -208,18 +211,18 @@ fn fri_combine_initial<F: Field + Extendable<D>, const D: usize>(
         .map(|&e| F::Extension::from_basefield(e));
     let wire_composition_eval = alpha.clone().reduce(wire_evals);
     let zeta_frob = zeta.frobenius();
-    let wire_eval = alpha.clone().reduce(os.wires.iter());
+    let mut alpha_frob = alpha.repeated_frobenius(D - 1);
+    let wire_eval = alpha.reduce(os.wires.iter());
     // We want to compute `sum a^i*phi(w_i)`, where `phi` denotes the Frobenius automorphism.
     // Since `phi^D=id` and `phi` is a field automorphism, we have the following equalities:
     // `sum a^i*phi(w_i) = sum phi(phi^(D-1)(a^i)*w_i) = phi(sum phi^(D-1)(a)^i*w_i)`
     // So we can compute the original sum using only one call to the `D-1`-repeated Frobenius of alpha,
     // and one call at the end of the sum.
-    let mut alpha_frob = alpha.repeated_frobenius(D - 1);
     let wire_eval_frob = alpha_frob.reduce(os.wires.iter()).frobenius();
-    let wire_interpol = interpolant(&[(zeta, wire_eval), (zeta_frob, wire_eval_frob)]);
-    let wire_numerator = wire_composition_eval - wire_interpol.eval(subgroup_x);
+    let wire_interpol = interpolate2([(zeta, wire_eval), (zeta_frob, wire_eval_frob)], subgroup_x);
+    let wire_numerator = wire_composition_eval - wire_interpol;
     let wire_denominator = (subgroup_x - zeta) * (subgroup_x - zeta_frob);
-    sum = alpha_frob.shift(sum);
+    sum = alpha.shift(sum);
     sum += wire_numerator / wire_denominator;
 
     sum

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::interpolation::interpolate2;
 use crate::fri::{prover::fri_proof, verifier::verify_fri_proof, FriConfig};
 use crate::merkle_tree::MerkleTree;
 use crate::plonk_challenger::Challenger;
@@ -125,17 +125,10 @@ impl<F: Field> ListPolynomialCommitment<F> {
         .iter()
         .flat_map(|&p| &commitments[p.index].polynomials)
         .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 = alpha.clone().reduce_polys(single_polys);
-        let single_composition_eval = alpha.reduce(single_evals);
+        let single_composition_poly = alpha.reduce_polys(single_polys);
 
-        let single_quotient = Self::compute_quotient(
-            &[zeta],
-            &[single_composition_eval],
-            &single_composition_poly,
-        );
-        final_poly = &final_poly + &single_quotient;
+        let single_quotient = Self::compute_quotient([zeta], single_composition_poly);
+        final_poly += single_quotient;
         alpha.reset();
 
         // Zs polynomials are opened at `zeta` and `g*zeta`.
@@ -143,19 +136,11 @@ impl<F: Field> ListPolynomialCommitment<F> {
             .polynomials
             .iter()
             .map(|p| p.to_extension());
-        let zs_composition_poly = alpha.clone().reduce_polys(zs_polys);
-        let zs_composition_evals = [
-            alpha.clone().reduce(os.plonk_zs.iter()),
-            alpha.reduce(os.plonk_zs_right.iter()),
-        ];
+        let zs_composition_poly = alpha.reduce_polys(zs_polys);
 
-        let zs_quotient = Self::compute_quotient(
-            &[zeta, g * zeta],
-            &zs_composition_evals,
-            &zs_composition_poly,
-        );
+        let zs_quotient = Self::compute_quotient([zeta, g * zeta], zs_composition_poly);
         final_poly = alpha.shift_poly(final_poly);
-        final_poly = &final_poly + &zs_quotient;
+        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
@@ -164,20 +149,12 @@ impl<F: Field> ListPolynomialCommitment<F> {
             .polynomials
             .iter()
             .map(|p| p.to_extension());
-        let wire_composition_poly = alpha.clone().reduce_polys(wire_polys);
-        let mut alpha_frob = alpha.repeated_frobenius(D - 1);
-        let wire_composition_evals = [
-            alpha.clone().reduce(os.wires.iter()),
-            alpha_frob.reduce(os.wires.iter()).frobenius(),
-        ];
+        let wire_composition_poly = alpha.reduce_polys(wire_polys);
 
-        let wires_quotient = Self::compute_quotient(
-            &[zeta, zeta.frobenius()],
-            &wire_composition_evals,
-            &wire_composition_poly,
-        );
-        final_poly = alpha_frob.shift_poly(final_poly);
-        final_poly = &final_poly + &wires_quotient;
+        let wires_quotient =
+            Self::compute_quotient([zeta, zeta.frobenius()], wire_composition_poly);
+        final_poly = alpha.shift_poly(final_poly);
+        final_poly += wires_quotient;
 
         let lde_final_poly = final_poly.lde(config.rate_bits);
         let lde_final_values = lde_final_poly
@@ -208,28 +185,27 @@ impl<F: Field> ListPolynomialCommitment<F> {
 
     /// 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_quotient<const D: usize>(
-        points: &[F::Extension],
-        evals: &[F::Extension],
-        poly: &PolynomialCoeffs<F::Extension>,
+    fn compute_quotient<const D: usize, const N: usize>(
+        points: [F::Extension; N],
+        poly: PolynomialCoeffs<F::Extension>,
     ) -> PolynomialCoeffs<F::Extension>
     where
         F: Extendable<D>,
     {
-        let pairs = points
-            .iter()
-            .zip(evals)
-            .map(|(&x, &e)| (x, e))
-            .collect::<Vec<_>>();
-        debug_assert!(pairs.iter().all(|&(x, e)| poly.eval(x) == e));
-
-        let interpolant = interpolant(&pairs);
-        let denominator = points.iter().fold(PolynomialCoeffs::one(), |acc, &x| {
-            &acc * &PolynomialCoeffs::new(vec![-x, F::Extension::ONE])
-        });
-        let numerator = poly - &interpolant;
-        let (quotient, rem) = numerator.div_rem(&denominator);
-        debug_assert!(rem.is_zero());
+        let quotient = if N == 1 {
+            poly.divide_by_linear(points[0]).0
+        } else if N == 2 {
+            // The denominator is `(X - p0)(X - p1) = p0 p1 - (p0 + p1) X + X^2`.
+            let denominator = vec![
+                points[0] * points[1],
+                -points[0] - points[1],
+                F::Extension::ONE,
+            ]
+            .into();
+            poly.div_rem_long_division(&denominator).0 // Could also use `divide_by_linear` twice.
+        } else {
+            unreachable!("This shouldn't happen. Plonk should open polynomials at 1 or 2 points.")
+        };
 
         quotient.padded(quotient.degree_plus_one().next_power_of_two())
     }