diff --git a/crates/prover/src/core/constraints.rs b/crates/prover/src/core/constraints.rs index f66c8d93d..0963db79c 100644 --- a/crates/prover/src/core/constraints.rs +++ b/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,26 @@ pub fn coset_vanishing>(coset: Coset, mut p: CirclePoi x } +/// Evaluates the formal derivative of a vanishing polynomial of a `coset` at a point `p`. +pub fn coset_vanishing_derivative>(coset: Coset, p: CirclePoint) -> F { + // The formal derivative is computed by multiplying the product of the vanishing polynomials of + // of all cosets of smaller size than the given coset at the given point by 4^(log_size-1). + // See Remark 15 in circle stark paper. + + let field_four = F::from(BaseField::from_u32_unchecked(4)); + 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. diff --git a/crates/prover/src/core/poly/circle/evaluation.rs b/crates/prover/src/core/poly/circle/evaluation.rs index 755ddc459..a527f975c 100644 --- a/crates/prover/src/core/poly/circle/evaluation.rs +++ b/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::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::qm31::SecureField; use crate::core::fields::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,80 @@ impl> Index for CosetSubEvaluation<'_, F> { } } +// TODO(Gali): Remove. +#[allow(dead_code)] +/// Computes the weights for Barycentric Lagrange interpolation for point `p` on `coset`. +/// `p` must not be in the domain. +fn barycentric_weights( + coset: CanonicCoset, + p: CirclePoint, +) -> Col { + // For a canonic coset `coset` of size 2^n and a point `p` not in `coset`, the weight at a coset + // point i is computed as: + // + // W_i = S_i(p) / S_i(i) = V_n(p) / (-2 * V'_n(i_x) * i_y * V_i(p)) + // + // using the following identities from the circle stark paper: + // + // S_i(p) = V_n(p) / V_i(p) + // S_i(i) = -2 * V'(i_x) * i_y + // + // where: + // - S_i(point) is the vanishing polynomial on the coset except i, evaluated at a point. + // - V_n(p) is the vanishing polynomial on the coset, 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`]. + + let domain = coset.circle_domain(); + + let (si_i, vi_p): (Vec<_>, Vec<_>) = (0..domain.size()) + .map(|i| { + let coset_point = domain.at(i).into_ef::(); + let minus_two_coset_point_y = coset_point.y * SecureField::from(-2); + ( + minus_two_coset_point_y + * coset_vanishing_derivative( + Coset::new(CirclePointIndex::generator(), domain.log_size()), + coset_point, + ), + point_vanishing(coset_point, p.into_ef::()), + ) + }) + .unzip(); + + let vn_p: SecureField = coset_vanishing( + CanonicCoset::new(domain.log_size()).coset, + p.into_ef::(), + ); + + // TODO(Gali): Change weights order to bit-reverse order. + (0..domain.size()) + .map(|i| vn_p / (si_i[i] * vi_p[i])) + .collect_vec() +} + +// TODO(Gali): Remove. +#[allow(dead_code)] +/// Evaluates a polynomial at a point using the barycentric interpolation formula, +/// given its evaluations on a circle domain and precomputed barycentric weights for the domain +/// at the sampled point. +fn barycentric_eval_at_point( + evals: &CircleEvaluation, + weights: &Col, +) -> SecureField { + // Evaluation = Σ W_i * Poly(i) for all i in the evaluation domain. + // For more information on barycentric weights calculation see [`barycentric_weights`] + (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 +279,36 @@ 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(10); + 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), + ]; + 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, &barycentric_weights(s, *point))) + .collect_vec(); + + assert_eq!( + sampled_barycentric_values, sampled_values, + "Barycentric evaluation should be equal to the polynomial evaluation" + ); + } }