CPU Barycentric Evaluation and Test

crates/prover/src/core/constraints.rs

@@ -5,6 +5,7 @@ use super::fields::m31::BaseField;
 use super::fields::qm31::SecureField;
 use super::fields::ExtensionOf;
 use super::pcs::quotients::PointSample;
+use super::poly::circle::CanonicCoset;
 use crate::core::fields::ComplexConjugate;
 
 /// Evaluates a vanishing polynomial of the coset at a point.
@@ -33,6 +34,27 @@ pub fn coset_vanishing<F: ExtensionOf<BaseField>>(coset: Coset, mut p: CirclePoi
     x
 }
 
+/// Evaluates the formal derivative of a vanishing polynomial of the coset at a point.
+pub fn coset_vanishing_derivative<F: ExtensionOf<BaseField>>(coset: Coset, p: CirclePoint<F>) -> F {
+    // The formal derivative is computed by multiplying the product of the vanishing polynomials of
+    // of all cosets of smallet size than the given coset at the given point by 4^(log_size-1),
+    // where log_size is the log size of the coset. This is because of the chain rule (also see
+    // remark 15 in circle stark article).
+
+    let field_four = F::one() + F::one() + F::one() + F::one();
+    let mut exp = F::one();
+    for _ in 1..coset.log_size {
+        exp *= field_four;
+    }
+
+    let mut vanishing = F::one();
+    for i in 1..coset.log_size {
+        vanishing *= coset_vanishing(CanonicCoset::new(i).coset, p)
+    }
+
+    exp * vanishing
+}
+
 /// Evaluates the polynomial that is used to exclude the excluded point at point
 /// p. Note that this polynomial has a zero of multiplicity 2 at the excluded
 /// point.

crates/prover/src/core/poly/circle/evaluation.rs

@@ -2,14 +2,19 @@ use std::marker::PhantomData;
 use std::ops::{Deref, Index};
 
 use educe::Educe;
+use itertools::Itertools;
+use num_traits::{One, Zero};
 
 use super::{CircleDomain, CirclePoly, PolyOps};
 use crate::core::backend::cpu::CpuCircleEvaluation;
 use crate::core::backend::simd::SimdBackend;
 use crate::core::backend::{Col, Column, ColumnOps, CpuBackend};
-use crate::core::circle::{CirclePointIndex, Coset};
+use crate::core::circle::{CirclePoint, CirclePointIndex, Coset};
+use crate::core::constraints::{coset_vanishing, coset_vanishing_derivative, point_vanishing};
 use crate::core::fields::m31::BaseField;
-use crate::core::fields::ExtensionOf;
+use crate::core::fields::qm31::SecureField;
+use crate::core::fields::{batch_inverse_in_place, ExtensionOf};
+use crate::core::poly::circle::CanonicCoset;
 use crate::core::poly::twiddles::TwiddleTree;
 use crate::core::poly::{BitReversedOrder, NaturalOrder};
 use crate::core::utils::bit_reverse_index;
@@ -155,10 +160,124 @@ impl<F: ExtensionOf<BaseField>> Index<usize> for CosetSubEvaluation<'_, F> {
     }
 }
 
+// TODO(Gali): Remove.
+#[allow(dead_code)]
+/// Computes the weights for Barycentric Lagrange interpolation on a circle domain, given a sample
+/// point p and domain size. The domain is the circle domain corresponding to the canonic coset of
+/// the size provided.
+fn weights(log_size: u32, sample_point: CirclePoint<SecureField>) -> Col<CpuBackend, SecureField> {
+    // For a point p not in the domain, the weight at domain point i is computed as:
+    // ```text
+    // w_i = S_i(p) / S_i(i) = V(p) / (-2 * V'_n(i_x) * i_y * V_i(p))
+    // ```
+    // using the following identities from the circle stark article:
+    // ```text
+    // S_i(p) = V(p) / V_i(p)
+    // S_i(i) = -2 * V'(i_x) * i_y
+    // ```
+    // where:
+    // - S_i(point) is the vanishing polynomial on the domain except i, evaluated at a point.
+    // - V(p) is the vanishing polynomial on the domain, evaluated at p.
+    // - V_i(p) is the vanishing polynomial on point i, evaluated at p.
+    // - V'(i_x) is the derivative of V(i) (evaluated at that point), see
+    //   [`coset_vanishing_derivative`].
+
+    // TODO(Gali): Change weights order to bit-reverse order.
+
+    let domain = CanonicCoset::new(log_size).circle_domain();
+
+    // If p is in the domain at position i, then w_j = δ_ij
+    for i in 0..domain.size() {
+        if domain.at(i).into_ef() == sample_point {
+            let mut weights = vec![SecureField::zero(); domain.size()];
+            weights[i] = SecureField::one();
+            return weights;
+        }
+    }
+
+    // S_i(i) is invariant under G_(n−1) and alternate under J, so we can calculate only 2 values
+    let p_0 = domain.at(0).into_ef::<SecureField>();
+    let weights_first_half = SecureField::one()
+        / (-(p_0.y + p_0.y)
+            * coset_vanishing_derivative(
+                Coset::new(CirclePointIndex::generator(), domain.log_size()),
+                p_0,
+            ));
+    let p_0_inverse = domain.at(domain.half_coset.size()).into_ef::<SecureField>();
+    let weights_second_half = SecureField::one()
+        / (-(p_0_inverse.y + p_0_inverse.y)
+            * coset_vanishing_derivative(
+                Coset::new(CirclePointIndex::generator(), domain.log_size()),
+                p_0_inverse,
+            ));
+
+    // TODO(Gali): Optimize to a batched point_vanishing()
+    let domain_points_vanishing_evaluated_at_point = (0..domain.size())
+        .map(|i| {
+            point_vanishing(
+                domain.at(i).into_ef::<SecureField>(),
+                sample_point.into_ef::<SecureField>(),
+            )
+        })
+        .collect_vec();
+    let mut inversed_domain_points_vanishing_evaluated_at_point =
+        vec![unsafe { std::mem::zeroed() }; domain.size()];
+
+    batch_inverse_in_place(
+        &domain_points_vanishing_evaluated_at_point,
+        &mut inversed_domain_points_vanishing_evaluated_at_point,
+    );
+
+    let coset_vanishing_evaluated_at_point: SecureField = coset_vanishing(
+        CanonicCoset::new(domain.log_size()).coset,
+        sample_point.into_ef::<SecureField>(),
+    );
+
+    (0..domain.size())
+        .map(|i| {
+            if i < domain.half_coset.size() {
+                weights_first_half
+                    * inversed_domain_points_vanishing_evaluated_at_point[i]
+                    * coset_vanishing_evaluated_at_point
+            } else {
+                weights_second_half
+                    * inversed_domain_points_vanishing_evaluated_at_point[i]
+                    * coset_vanishing_evaluated_at_point
+            }
+        })
+        .collect_vec()
+}
+
+// TODO(Gali): Remove.
+#[allow(dead_code)]
+/// Evaluates a polynomial at a point using the barycentric interpolation formula,
+/// given its evaluations on a domain and precomputed barycentric weights for
+/// the domain.
+fn barycentric_eval_at_point(
+    evals: &CircleEvaluation<CpuBackend, BaseField, BitReversedOrder>,
+    point: CirclePoint<SecureField>,
+    weights: &Col<CpuBackend, SecureField>,
+) -> SecureField {
+    // ```text
+    // Evaluation = Σ W_i * Poly(i) for all i in the evaluation domain.
+    // ```
+    // For more information on barycentric weights calculation see [`weights`]
+    for i in 0..evals.domain.size() {
+        if point == evals.domain.at(i).into_ef() {
+            return evals.values[bit_reverse_index(i, evals.domain.log_size())].into();
+        }
+    }
+
+    (0..evals.domain.size()).fold(SecureField::zero(), |acc, i| {
+        acc + (evals.values[bit_reverse_index(i, evals.domain.log_size())] * weights[i])
+    })
+}
+
 #[cfg(test)]
 mod tests {
-    use crate::core::backend::cpu::CpuCircleEvaluation;
-    use crate::core::circle::Coset;
+    use super::*;
+    use crate::core::backend::cpu::{CpuCircleEvaluation, CpuCirclePoly};
+    use crate::core::circle::{CirclePoint, Coset};
     use crate::core::fields::m31::BaseField;
     use crate::core::poly::circle::CanonicCoset;
     use crate::core::poly::NaturalOrder;
@@ -204,4 +323,40 @@ mod tests {
             assert_eq!(sub_eval[i], circle_evaluation.get_at(coset.index_at(i)));
         }
     }
+
+    #[test]
+    fn test_barycentric_evaluation() {
+        let poly = CpuCirclePoly::new(
+            [691, 805673, 5, 435684, 4832, 23876431, 197, 897346068]
+                .map(BaseField::from)
+                .to_vec(),
+        );
+        let s = CanonicCoset::new(3);
+        let domain = s.circle_domain();
+        let eval = poly.evaluate(domain);
+        let sampled_points = [
+            CirclePoint::get_point(348),
+            CirclePoint::get_point(9736524),
+            CirclePoint::get_point(13),
+            CirclePoint::get_point(346752),
+            domain.at(0).into_ef(),
+            domain.at(3).into_ef(),
+        ];
+        let sampled_values = sampled_points
+            .iter()
+            .map(|point| poly.eval_at_point(*point))
+            .collect::<Vec<_>>();
+
+        let sampled_barycentric_values = sampled_points
+            .iter()
+            .map(|point| {
+                barycentric_eval_at_point(&eval, *point, &weights(eval.domain.log_size(), *point))
+            })
+            .collect::<Vec<_>>();
+
+        assert_eq!(
+            sampled_barycentric_values, sampled_values,
+            "Barycentric evaluation should be equal to the polynomial evaluation"
+        );
+    }
 }