From b6f987f8a65feac82ac519b6a56004df3ddfc7be Mon Sep 17 00:00:00 2001 From: Bobbin Threadbare <43513081+bobbinth@users.noreply.github.com> Date: Tue, 25 Jun 2024 23:12:36 -0700 Subject: [PATCH] feat: initial implementation of LogUp-GKR bus MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Al-Kindi-0 <82364884+Al-Kindi-0@users.noreply.github.com> Co-authored-by: Philippe Laferrière --- air/src/trace/main_trace.rs | 12 +- processor/Cargo.toml | 4 +- processor/src/lib.rs | 5 +- processor/src/trace/mod.rs | 3 + .../src/trace/virtual_bus/circuit/error.rs | 24 + .../src/trace/virtual_bus/circuit/mod.rs | 290 +++++++++ .../src/trace/virtual_bus/circuit/prover.rs | 587 ++++++++++++++++++ .../src/trace/virtual_bus/circuit/verifier.rs | 237 +++++++ processor/src/trace/virtual_bus/mod.rs | 39 ++ .../trace/virtual_bus/multilinear/error.rs | 5 + .../virtual_bus/multilinear/lagrange_ker.rs | 103 +++ .../src/trace/virtual_bus/multilinear/mod.rs | 108 ++++ .../src/trace/virtual_bus/sum_check/mod.rs | 250 ++++++++ .../virtual_bus/sum_check/prover/error.rs | 15 + .../trace/virtual_bus/sum_check/prover/mod.rs | 329 ++++++++++ .../virtual_bus/sum_check/verifier/error.rs | 9 + .../virtual_bus/sum_check/verifier/mod.rs | 152 +++++ processor/src/trace/virtual_bus/tests.rs | 81 +++ processor/src/trace/virtual_bus/univariate.rs | 253 ++++++++ 19 files changed, 2500 insertions(+), 6 deletions(-) create mode 100644 processor/src/trace/virtual_bus/circuit/error.rs create mode 100644 processor/src/trace/virtual_bus/circuit/mod.rs create mode 100644 processor/src/trace/virtual_bus/circuit/prover.rs create mode 100644 processor/src/trace/virtual_bus/circuit/verifier.rs create mode 100644 processor/src/trace/virtual_bus/mod.rs create mode 100644 processor/src/trace/virtual_bus/multilinear/error.rs create mode 100644 processor/src/trace/virtual_bus/multilinear/lagrange_ker.rs create mode 100644 processor/src/trace/virtual_bus/multilinear/mod.rs create mode 100644 processor/src/trace/virtual_bus/sum_check/mod.rs create mode 100644 processor/src/trace/virtual_bus/sum_check/prover/error.rs create mode 100644 processor/src/trace/virtual_bus/sum_check/prover/mod.rs create mode 100644 processor/src/trace/virtual_bus/sum_check/verifier/error.rs create mode 100644 processor/src/trace/virtual_bus/sum_check/verifier/mod.rs create mode 100644 processor/src/trace/virtual_bus/tests.rs create mode 100644 processor/src/trace/virtual_bus/univariate.rs diff --git a/air/src/trace/main_trace.rs b/air/src/trace/main_trace.rs index c25d8716c8..d3b23b8001 100644 --- a/air/src/trace/main_trace.rs +++ b/air/src/trace/main_trace.rs @@ -19,9 +19,6 @@ use super::{ use core::ops::{Deref, Range}; use vm_core::{utils::range, Felt, Word, ONE, ZERO}; -#[cfg(any(test, feature = "internals"))] -use alloc::vec::Vec; - // CONSTANTS // ================================================================================================ @@ -43,6 +40,13 @@ impl Deref for MainTrace { } } +#[cfg(any(test, feature = "internals"))] +impl core::ops::DerefMut for MainTrace { + fn deref_mut(&mut self) -> &mut Self::Target { + &mut self.columns + } +} + impl MainTrace { pub fn new(main_trace: ColMatrix) -> Self { Self { @@ -55,7 +59,7 @@ impl MainTrace { } #[cfg(any(test, feature = "internals"))] - pub fn get_column_range(&self, range: Range) -> Vec> { + pub fn get_column_range(&self, range: Range) -> alloc::vec::Vec> { range.fold(vec![], |mut acc, col_idx| { acc.push(self.get_column(col_idx).to_vec()); acc diff --git a/processor/Cargo.toml b/processor/Cargo.toml index 83da03a0cf..1e5f963906 100644 --- a/processor/Cargo.toml +++ b/processor/Cargo.toml @@ -26,8 +26,10 @@ std = ["vm-core/std", "winter-prover/std"] tracing = { version = "0.1", default-features = false, features = [ "attributes", ] } -vm-core = { package = "miden-core", path = "../core", version = "0.9", default-features = false } miden-air = { package = "miden-air", path = "../air", version = "0.9", default-features = false } +static_assertions = "1.1.0" +thiserror = { version = "1.0", default-features = false } +vm-core = { package = "miden-core", path = "../core", version = "0.9", default-features = false } winter-prover = { package = "winter-prover", version = "0.9", default-features = false } [dev-dependencies] diff --git a/processor/src/lib.rs b/processor/src/lib.rs index fc6b136a2b..583be294f4 100644 --- a/processor/src/lib.rs +++ b/processor/src/lib.rs @@ -61,7 +61,10 @@ use chiplets::Chiplets; mod trace; use trace::TraceFragment; -pub use trace::{ChipletsLengths, ExecutionTrace, TraceLenSummary, NUM_RAND_ROWS}; +pub use trace::{ + prove_virtual_bus, verify_virtual_bus, ChipletsLengths, ExecutionTrace, TraceLenSummary, + NUM_RAND_ROWS, +}; mod errors; pub use errors::{ExecutionError, Ext2InttError}; diff --git a/processor/src/trace/mod.rs b/processor/src/trace/mod.rs index 888d6064e1..1d4ace0f73 100644 --- a/processor/src/trace/mod.rs +++ b/processor/src/trace/mod.rs @@ -18,6 +18,9 @@ use winter_prover::{crypto::RandomCoin, EvaluationFrame, Trace, TraceInfo}; mod utils; pub use utils::{AuxColumnBuilder, ChipletsLengths, TraceFragment, TraceLenSummary}; +mod virtual_bus; +pub use virtual_bus::{prove as prove_virtual_bus, verify as verify_virtual_bus}; + #[cfg(test)] mod tests; #[cfg(test)] diff --git a/processor/src/trace/virtual_bus/circuit/error.rs b/processor/src/trace/virtual_bus/circuit/error.rs new file mode 100644 index 0000000000..bfc2a35f91 --- /dev/null +++ b/processor/src/trace/virtual_bus/circuit/error.rs @@ -0,0 +1,24 @@ +use crate::trace::virtual_bus::sum_check::SumCheckProverError; +use crate::trace::virtual_bus::sum_check::SumCheckVerifierError; + +#[derive(Debug, thiserror::Error)] +pub enum ProverError { + #[error("failed to generate multi-linear from the given evaluations")] + FailedToGenerateML, + #[error("failed to generate the sum-check proof")] + FailedToProveSumCheck(#[from] SumCheckProverError), + #[error("failed to generate the random challenge")] + FailedToGenerateChallenge, +} + +#[derive(Debug, thiserror::Error)] +pub enum VerifierError { + #[error("one of the claimed circuit denominators is zero")] + ZeroOutputDenominator, + #[error("the output of the fraction circuit is not equal to the expected value")] + MismatchingCircuitOutput, + #[error("failed to generate the random challenge")] + FailedToGenerateChallenge, + #[error("failed to verify the sum-check proof")] + FailedToVerifySumCheck(#[from] SumCheckVerifierError), +} diff --git a/processor/src/trace/virtual_bus/circuit/mod.rs b/processor/src/trace/virtual_bus/circuit/mod.rs new file mode 100644 index 0000000000..5737693bab --- /dev/null +++ b/processor/src/trace/virtual_bus/circuit/mod.rs @@ -0,0 +1,290 @@ +use core::ops::Add; + +use crate::trace::virtual_bus::multilinear::EqFunction; +use crate::trace::virtual_bus::sum_check::{CompositionPolynomial, RoundProof}; +use alloc::vec::Vec; +use miden_air::trace::chiplets::{MEMORY_D0_COL_IDX, MEMORY_D1_COL_IDX}; +use miden_air::trace::decoder::{DECODER_OP_BITS_OFFSET, DECODER_USER_OP_HELPERS_OFFSET}; +use miden_air::trace::range::{M_COL_IDX, V_COL_IDX}; +use miden_air::trace::{CHIPLETS_OFFSET, TRACE_WIDTH}; +use prover::CircuitLayerPolys; +use static_assertions::const_assert; +use vm_core::{Felt, FieldElement}; + +mod error; +mod prover; +pub use prover::prove; + +mod verifier; +pub use verifier::verify; + +use super::multilinear::MultiLinearPoly; +use super::sum_check::{FinalOpeningClaim, Proof as SumCheckProof}; + +/// Defines the number of wires in the input layer that are generated from a single main trace row. +const NUM_WIRES_PER_TRACE_ROW: usize = 8; +const_assert!(NUM_WIRES_PER_TRACE_ROW.is_power_of_two()); + +// CIRCUIT WIRE +// ================================================================================================ + +/// Represents a fraction `numerator / denominator` as a pair `(numerator, denominator)`. This is +/// the type for the gates' inputs in [`prover::EvaluatedCircuit`]. +/// +/// Hence, addition is defined in the natural way fractions are added together: `a/b + c/d = (ad + +/// bc) / bd`. +#[derive(Debug, Clone, Copy)] +pub struct CircuitWire { + numerator: E, + denominator: E, +} + +impl CircuitWire +where + E: FieldElement, +{ + /// Creates new projective coordinates from a numerator and a denominator. + pub fn new(numerator: E, denominator: E) -> Self { + assert_ne!(denominator, E::ZERO); + + Self { + numerator, + denominator, + } + } +} + +impl Add for CircuitWire +where + E: FieldElement, +{ + type Output = Self; + + fn add(self, other: Self) -> Self { + let numerator = self.numerator * other.denominator + other.numerator * self.denominator; + let denominator = self.denominator * other.denominator; + + Self::new(numerator, denominator) + } +} + +/// Converts a main trace row (or more generally "query") to numerators and denominators of the +/// input layer. +fn evaluate_fractions_at_main_trace_query( + query: &[E], + log_up_randomness: &[E], +) -> [[E; NUM_WIRES_PER_TRACE_ROW]; 2] +where + E: FieldElement, +{ + // numerators + let multiplicity = query[M_COL_IDX]; + let f_m = { + let mem_selec0 = query[CHIPLETS_OFFSET]; + let mem_selec1 = query[CHIPLETS_OFFSET + 1]; + let mem_selec2 = query[CHIPLETS_OFFSET + 2]; + mem_selec0 * mem_selec1 * (E::ONE - mem_selec2) + }; + + let f_rc = { + let op_bit_4 = query[DECODER_OP_BITS_OFFSET + 4]; + let op_bit_5 = query[DECODER_OP_BITS_OFFSET + 5]; + let op_bit_6 = query[DECODER_OP_BITS_OFFSET + 6]; + + (E::ONE - op_bit_4) * (E::ONE - op_bit_5) * op_bit_6 + }; + + // denominators + let alphas = log_up_randomness; + + let table_denom = alphas[0] - query[V_COL_IDX]; + let memory_denom_0 = -(alphas[0] - query[MEMORY_D0_COL_IDX]); + let memory_denom_1 = -(alphas[0] - query[MEMORY_D1_COL_IDX]); + let stack_value_denom_0 = -(alphas[0] - query[DECODER_USER_OP_HELPERS_OFFSET]); + let stack_value_denom_1 = -(alphas[0] - query[DECODER_USER_OP_HELPERS_OFFSET + 1]); + let stack_value_denom_2 = -(alphas[0] - query[DECODER_USER_OP_HELPERS_OFFSET + 2]); + let stack_value_denom_3 = -(alphas[0] - query[DECODER_USER_OP_HELPERS_OFFSET + 3]); + + [ + [multiplicity, f_m, f_m, f_rc, f_rc, f_rc, f_rc, E::ZERO], + [ + table_denom, + memory_denom_0, + memory_denom_1, + stack_value_denom_0, + stack_value_denom_1, + stack_value_denom_2, + stack_value_denom_3, + E::ONE, + ], + ] +} + +/// Computes the wires added to the input layer that come from a given main trace row (or more +/// generally, "query"). +fn compute_input_layer_wires_at_main_trace_query( + query: &[E], + log_up_randomness: &[E], +) -> [CircuitWire; NUM_WIRES_PER_TRACE_ROW] +where + E: FieldElement, +{ + let [numerators, denominators] = + evaluate_fractions_at_main_trace_query(query, log_up_randomness); + let input_gates_values: Vec> = numerators + .into_iter() + .zip(denominators) + .map(|(numerator, denominator)| CircuitWire::new(numerator, denominator)) + .collect(); + input_gates_values.try_into().unwrap() +} + +/// A GKR proof for the correct evaluation of the sum of fractions circuit. +#[derive(Debug)] +pub struct GkrCircuitProof { + circuit_outputs: CircuitLayerPolys, + before_final_layer_proofs: BeforeFinalLayerProof, + final_layer_proof: FinalLayerProof, +} + +impl GkrCircuitProof { + pub fn get_final_opening_claim(&self) -> FinalOpeningClaim { + self.final_layer_proof.after_merge_proof.openings_claim.clone() + } +} + +/// A set of sum-check proofs for all GKR layers but for the input circuit layer. +#[derive(Debug)] +pub struct BeforeFinalLayerProof { + pub proof: Vec>, +} + +/// A proof for the input circuit layer i.e., the final layer in the GKR protocol. +#[derive(Debug)] +pub struct FinalLayerProof { + before_merge_proof: Vec>, + after_merge_proof: SumCheckProof, +} + +/// Represents a claim to be proven by a subsequent call to the sum-check protocol. +#[derive(Debug)] +pub struct GkrClaim { + pub evaluation_point: Vec, + pub claimed_evaluation: (E, E), +} + +/// A composition polynomial used in the GKR protocol for all of its sum-checks except the final +/// one. +#[derive(Clone)] +pub struct GkrComposition +where + E: FieldElement, +{ + pub combining_randomness: E, +} + +impl GkrComposition +where + E: FieldElement, +{ + pub fn new(combining_randomness: E) -> Self { + Self { + combining_randomness, + } + } +} + +impl CompositionPolynomial for GkrComposition +where + E: FieldElement, +{ + fn max_degree(&self) -> u32 { + 3 + } + + fn evaluate(&self, query: &[E]) -> E { + let eval_left_numerator = query[0]; + let eval_right_numerator = query[1]; + let eval_left_denominator = query[2]; + let eval_right_denominator = query[3]; + let eq_eval = query[4]; + eq_eval + * ((eval_left_numerator * eval_right_denominator + + eval_right_numerator * eval_left_denominator) + + eval_left_denominator * eval_right_denominator * self.combining_randomness) + } +} + +/// A composition polynomial used in the GKR protocol for its final sum-check. +#[derive(Clone)] +pub struct GkrCompositionMerge +where + E: FieldElement, +{ + pub sum_check_combining_randomness: E, + pub tensored_merge_randomness: Vec, + pub log_up_randomness: Vec, +} + +impl GkrCompositionMerge +where + E: FieldElement, +{ + pub fn new( + combining_randomness: E, + merge_randomness: Vec, + log_up_randomness: Vec, + ) -> Self { + let tensored_merge_randomness = + EqFunction::ml_at(merge_randomness.clone()).evaluations().to_vec(); + + Self { + sum_check_combining_randomness: combining_randomness, + tensored_merge_randomness, + log_up_randomness, + } + } +} + +impl CompositionPolynomial for GkrCompositionMerge +where + E: FieldElement, +{ + fn max_degree(&self) -> u32 { + // Computed as: + // 1 + max(left_numerator_degree + right_denom_degree, right_numerator_degree + + // left_denom_degree) + 5 + } + + fn evaluate(&self, query: &[E]) -> E { + let [numerators, denominators] = + evaluate_fractions_at_main_trace_query(query, &self.log_up_randomness); + + let numerators = MultiLinearPoly::from_evaluations(numerators.to_vec()).unwrap(); + let denominators = MultiLinearPoly::from_evaluations(denominators.to_vec()).unwrap(); + + let (left_numerators, right_numerators) = numerators.project_least_significant_variable(); + let (left_denominators, right_denominators) = + denominators.project_least_significant_variable(); + + let eval_left_numerators = + left_numerators.evaluate_with_lagrange_kernel(&self.tensored_merge_randomness); + let eval_right_numerators = + right_numerators.evaluate_with_lagrange_kernel(&self.tensored_merge_randomness); + + let eval_left_denominators = + left_denominators.evaluate_with_lagrange_kernel(&self.tensored_merge_randomness); + let eval_right_denominators = + right_denominators.evaluate_with_lagrange_kernel(&self.tensored_merge_randomness); + + let eq_eval = query[TRACE_WIDTH]; + + eq_eval + * ((eval_left_numerators * eval_right_denominators + + eval_right_numerators * eval_left_denominators) + + eval_left_denominators + * eval_right_denominators + * self.sum_check_combining_randomness) + } +} diff --git a/processor/src/trace/virtual_bus/circuit/prover.rs b/processor/src/trace/virtual_bus/circuit/prover.rs new file mode 100644 index 0000000000..083721da44 --- /dev/null +++ b/processor/src/trace/virtual_bus/circuit/prover.rs @@ -0,0 +1,587 @@ +use super::{ + super::sum_check::Proof as SumCheckProof, compute_input_layer_wires_at_main_trace_query, + error::ProverError, BeforeFinalLayerProof, CircuitWire, FinalLayerProof, GkrCircuitProof, + GkrClaim, GkrComposition, GkrCompositionMerge, NUM_WIRES_PER_TRACE_ROW, +}; +use crate::trace::virtual_bus::{ + multilinear::{EqFunction, MultiLinearPoly}, + sum_check::{FinalClaimBuilder, FinalOpeningClaim, RoundClaim, RoundProof}, + SumCheckProver, +}; +use alloc::vec::Vec; +use core::marker::PhantomData; +use miden_air::trace::main_trace::MainTrace; +use vm_core::{Felt, FieldElement}; +use winter_prover::crypto::{ElementHasher, RandomCoin}; + +// EVALUATED CIRCUIT +// ================================================================================================ + +/// Evaluation of a layered circuit for computing a sum of fractions. +/// +/// The circuit computes a sum of fractions based on the formula a / c + b / d = (a * d + b * c) / +/// (c * d) which defines a "gate" ((a, b), (c, d)) --> (a * d + b * c, c * d) upon which the +/// [`EvaluatedCircuit`] is built. Due to the uniformity of the circuit, each of the circuit +/// layers collect all the: +/// +/// 1. `a`'s into a [`MultiLinearPoly`] called `left_numerators`. +/// 2. `b`'s into a [`MultiLinearPoly`] called `right_numerators`. +/// 3. `c`'s into a [`MultiLinearPoly`] called `left_denominators`. +/// 4. `d`'s into a [`MultiLinearPoly`] called `right_denominators`. +/// +/// The relation between two subsequent layers is given by the formula +/// +/// p_0[layer + 1](x_0, x_1, ..., x_{ν - 2}) = p_0[layer](x_0, x_1, ..., x_{ν - 2}, 0) * +/// q_1[layer](x_0, x_1, ..., x_{ν - 2}, 0) +/// + p_1[layer](x_0, x_1, ..., x_{ν - 2}, 0) * q_0[layer](x_0, +/// x_1, ..., x_{ν - 2}, 0) +/// +/// p_1[layer + 1](x_0, x_1, ..., x_{ν - 2}) = p_0[layer](x_0, x_1, ..., x_{ν - 2}, 1) * +/// q_1[layer](x_0, x_1, ..., x_{ν - 2}, 1) +/// + p_1[layer](x_0, x_1, ..., x_{ν - 2}, 1) * q_0[layer](x_0, +/// x_1, ..., x_{ν - 2}, 1) +/// +/// and +/// +/// q_0[layer + 1](x_0, x_1, ..., x_{ν - 2}) = q_0[layer](x_0, x_1, ..., x_{ν - 2}, 0) * +/// q_1[layer](x_0, x_1, ..., x_{ν - 1}, 0) +/// q_1[layer + 1](x_0, x_1, ..., x_{ν - 2}) = q_0[layer](x_0, x_1, ..., x_{ν - 2}, 1) * +/// q_1[layer](x_0, x_1, ..., x_{ν - 1}, 1) +/// +/// This logic is encoded in [`CircuitWire`]. +/// +/// This means that layer ν will be the output layer and will consist of four values +/// (p_0[ν - 1], p_1[ν - 1], p_0[ν - 1], p_1[ν - 1]) ∈ 𝔽^ν. +pub struct EvaluatedCircuit { + layer_polys: Vec>, +} + +impl EvaluatedCircuit { + /// Creates a new [`EvaluatedCircuit`] by evaluating the circuit where the input layer is + /// defined from the main trace columns. + pub fn new( + main_trace_columns: &[MultiLinearPoly], + log_up_randomness: &[E], + ) -> Result { + let mut layer_polys = Vec::new(); + + let mut current_layer = Self::generate_input_layer(main_trace_columns, log_up_randomness); + while current_layer.num_wires() > 1 { + let next_layer = Self::compute_next_layer(¤t_layer); + + layer_polys.push(CircuitLayerPolys::from_circuit_layer(current_layer)); + + current_layer = next_layer; + } + + Ok(Self { layer_polys }) + } + + /// Returns a layer of the evaluated circuit. + /// + /// Note that the return type is [`LayerPolys`] as opposed to [`Layer`], since the evaluated + /// circuit is stored in a representation which can be proved using GKR. + pub fn get_layer(&self, layer_idx: usize) -> &CircuitLayerPolys { + &self.layer_polys[layer_idx] + } + + /// Returns all layers of the evaluated circuit, starting from the input layer. + /// + /// Note that the return type is a slice of [`CircuitLayerPolys`] as opposed to + /// [`CircuitLayer`], since the evaluated layers are stored in a representation which can be + /// proved using GKR. + pub fn layers(&self) -> &[CircuitLayerPolys] { + &self.layer_polys + } + + /// Returns the numerator/denominator polynomials representing the output layer of the circuit. + pub fn output_layer(&self) -> &CircuitLayerPolys { + self.layer_polys.last().expect("circuit has at least one layer") + } + + /// Evaluates the output layer at `query`, where the numerators of the output layer are treated + /// as evaluations of a multilinear polynomial, and similarly for the denominators. + pub fn evaluate_output_layer(&self, query: E) -> (E, E) { + let CircuitLayerPolys { + numerators, + denominators, + } = self.output_layer(); + + (numerators.evaluate(&[query]), denominators.evaluate(&[query])) + } + + // HELPERS + // ------------------------------------------------------------------------------------------- + + /// Generates the input layer of the circuit from the main trace columns and some randomness + /// provided by the verifier. + fn generate_input_layer( + main_trace_columns: &[MultiLinearPoly], + log_up_randomness: &[E], + ) -> CircuitLayer { + let num_evaluations = main_trace_columns[0].num_evaluations(); + let mut input_layer_wires = Vec::with_capacity(num_evaluations * NUM_WIRES_PER_TRACE_ROW); + + for i in 0..num_evaluations { + let wires_from_trace_row = { + let query: Vec = main_trace_columns.iter().map(|ml| ml[i]).collect(); + compute_input_layer_wires_at_main_trace_query(&query, log_up_randomness) + }; + + input_layer_wires.extend(wires_from_trace_row); + } + + CircuitLayer::new(input_layer_wires) + } + + /// Computes the subsequent layer of the circuit from a given layer. + fn compute_next_layer(prev_layer: &CircuitLayer) -> CircuitLayer { + let next_layer_wires = prev_layer + .wires() + .chunks_exact(2) + .map(|input_wires| { + let left_input_wire = input_wires[0]; + let right_input_wire = input_wires[1]; + + // output wire + left_input_wire + right_input_wire + }) + .collect(); + + CircuitLayer::new(next_layer_wires) + } +} + +// CIRCUIT LAYER +// ================================================================================================ + +/// Represents a layer in a [`EvaluatedCircuit`]. +/// +/// A layer is made up of a set of `n` wires, where `n` is a power of two. This is the natural +/// circuit representation of a layer, where each consecutive pair of wires are summed to yield a +/// wire in the subsequent layer of an [`EvaluatedCircuit`]. +/// +/// Note that a [`Layer`] needs to be first converted to a [`LayerPolys`] before the evaluation of +/// the layer can be proved using GKR. +struct CircuitLayer { + wires: Vec>, +} + +impl CircuitLayer { + /// Creates a new [`Layer`] from a set of projective coordinates. + /// + /// Panics if the number of projective coordinates is not a power of two. + pub fn new(wires: Vec>) -> Self { + assert!(wires.len().is_power_of_two()); + + Self { wires } + } + + /// Returns the wires that make up this circuit layer. + pub fn wires(&self) -> &[CircuitWire] { + &self.wires + } + + /// Returns the number of wires in the layer. + pub fn num_wires(&self) -> usize { + self.wires.len() + } +} + +// CIRCUIT LAYER POLYNOMIAL +// ================================================================================================ + +/// Holds a layer of an [`EvaluatedCircuit`] in a representation amenable to proving circuit +/// evaluation using GKR. +#[derive(Clone, Debug)] +pub struct CircuitLayerPolys { + pub numerators: MultiLinearPoly, + pub denominators: MultiLinearPoly, +} + +impl CircuitLayerPolys +where + E: FieldElement, +{ + fn from_circuit_layer(layer: CircuitLayer) -> Self { + Self::from_wires(layer.wires) + } + + pub fn from_wires(wires: Vec>) -> Self { + let mut numerators = Vec::new(); + let mut denominators = Vec::new(); + + for wire in wires { + numerators.push(wire.numerator); + denominators.push(wire.denominator); + } + + Self { + numerators: MultiLinearPoly::from_evaluations(numerators) + .expect("evaluations guaranteed to be a power of two"), + denominators: MultiLinearPoly::from_evaluations(denominators) + .expect("evaluations guaranteed to be a power of two"), + } + } +} + +// PROVER +// ================================================================================================ + +/// Evaluates and proves a fractional sum circuit given a set of composition polynomials. +/// +/// For the input layer of the circuit, each individual component of the quadruple +/// [p_0, p_1, q_0, q_1] is of the form: +/// +/// m(z_0, ... , z_{μ - 1}, x_0, ... , x_{ν - 1}) = \sum_{y ∈ {0,1}^μ} EQ(z, y) * g_{[y]}(f_0(x_0, +/// ... , x_{ν - 1}), ... , f_{κ - 1}(x_0, ... , x_{ν +/// - 1})) +/// +/// where: +/// +/// 1. μ is the log_2 of the number of different numerator/denominator expressions divided by two. +/// 2. [y] := \sum_{j = 0}^{μ - 1} y_j * 2^j +/// 3. κ is the number of multi-linears (i.e., main trace columns) involved in the computation of +/// the circuit (i.e., virtual bus). +/// 4. ν is the log_2 of the trace length. +/// +/// The above `m` is usually referred to as the merge of the individual composed multi-linear +/// polynomials g_{[y]}(f_0(x_0, ... , x_{ν - 1}), ... , f_{κ - 1}(x_0, ... , x_{ν - 1})). +/// +/// The composition polynomials `g` are provided as inputs and then used in order to compute the +/// evaluations of each of the four merge polynomials over {0, 1}^{μ + ν}. The resulting evaluations +/// are then used in order to evaluate the circuit. At this point, the GKR protocol is used to prove +/// the correctness of circuit evaluation. It should be noted that the input layer, which +/// corresponds to the last layer treated by the GKR protocol, is handled differently from the other +/// layers. More specifically, the sum-check protocol used for the input layer is composed of two +/// sum-check protocols, the first one works directly with the evaluations of the `m`'s over {0, +/// 1}^{μ + ν} and runs for μ rounds. After these μ rounds, and using the resulting [`RoundClaim`], +/// we run the second and final sum-check protocol for ν rounds on the composed multi-linear +/// polynomial given by +/// +/// \sum_{y ∈ {0,1}^μ} EQ(ρ', y) * g_{[y]}(f_0(x_0, ... , x_{ν - 1}), ... , f_{κ - 1}(x_0, ... , +/// x_{ν - 1})) +/// +/// where ρ' is the randomness sampled during the first sum-check protocol. +/// +/// As part of the final sum-check protocol, the openings {f_j(ρ)} are provided as part of a +/// [`FinalOpeningClaim`]. This latter claim will be proven by the STARK prover later on using the +/// auxiliary trace. +pub fn prove< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + trace: &MainTrace, + log_up_randomness: Vec, + transcript: &mut C, +) -> Result, ProverError> { + // TODO: Optimize this so that we can work with base field element directly and thus save + // on memory usage. + let main_trace_columns: Vec> = trace + .columns() + .map(|col| { + let mut values: Vec = col.iter().map(|value| E::from(*value)).collect(); + if let Some(value) = values.last_mut() { + *value = E::ZERO + } + MultiLinearPoly::from_evaluations(values).unwrap() + }) + .collect(); + + // evaluate the GKR fractional sum circuit + let mut circuit = EvaluatedCircuit::new(&main_trace_columns, &log_up_randomness)?; + + // run the GKR prover for all layers except the input layer + let (before_final_layer_proofs, gkr_claim) = + prove_before_final_circuit_layers(&mut circuit, transcript)?; + + // run the GKR prover for the input layer + let num_rounds_before_merge = NUM_WIRES_PER_TRACE_ROW.ilog2() as usize - 1; + let final_layer_proof = prove_final_circuit_layer( + log_up_randomness, + main_trace_columns, + num_rounds_before_merge, + gkr_claim, + &mut circuit, + transcript, + )?; + + // include the circuit output as part of the final proof + let circuit_outputs = circuit.output_layer(); + + Ok(GkrCircuitProof { + circuit_outputs: circuit_outputs.clone(), + before_final_layer_proofs, + final_layer_proof, + }) +} + +// HELPER FUNCTIONS +// ================================================================================================ + +/// Proves the final GKR layer which corresponds to the input circuit layer. +fn prove_final_circuit_layer< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + log_up_randomness: Vec, + mut mls: Vec>, + num_rounds_merge: usize, + gkr_claim: GkrClaim, + circuit: &mut EvaluatedCircuit, + transcript: &mut C, +) -> Result, ProverError> { + // parse the [GkrClaim] resulting from the previous GKR layer + let GkrClaim { + evaluation_point, + claimed_evaluation, + } = gkr_claim; + + // compute the EQ function at the evaluation point + let poly_x = EqFunction::ml_at(evaluation_point.clone()); + + // get the multi-linears of the 4 merge polynomials + let layer = circuit.get_layer(0); + let (left_numerators, right_numerators) = layer.numerators.project_least_significant_variable(); + let (left_denominators, right_denominators) = + layer.denominators.project_least_significant_variable(); + let mut merged_mls = + vec![left_numerators, right_numerators, left_denominators, right_denominators, poly_x]; + // run the first sum-check protocol + let ((round_claim, before_merge_proof), r_sum_check) = sum_check_prover_plain_partial( + claimed_evaluation, + num_rounds_merge, + &mut merged_mls, + transcript, + )?; + + // parse the output of the first sum-check protocol + let RoundClaim { + eval_point: rand_merge, + claim, + } = round_claim; + + // create the composed multi-linear for the second sum-check protocol using the randomness + // sampled during the first one + let gkr_composition = GkrCompositionMerge::new(r_sum_check, rand_merge, log_up_randomness); + + // include the partially evaluated at the first sum-check randomness EQ multi-linear + // TODO: Find a better way than to push the evaluation of `EqFunction` here. + mls.push(merged_mls[4].clone()); + + // run the second sum-check protocol + let main_prover = SumCheckProver::new(gkr_composition, SimpleGkrFinalClaimBuilder(PhantomData)); + let after_merge_proof = main_prover.prove(claim, mls, transcript)?; + + Ok(FinalLayerProof { + before_merge_proof, + after_merge_proof, + }) +} + +/// Proves all GKR layers except for input layer. +fn prove_before_final_circuit_layers< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + circuit: &mut EvaluatedCircuit, + transcript: &mut C, +) -> Result<(BeforeFinalLayerProof, GkrClaim), ProverError> { + // absorb the circuit output layer. This corresponds to sending the four values of the output + // layer to the verifier. The verifier then replies with a challenge `r` in order to evaluate + // `p` and `q` at `r` as multi-linears. + let CircuitLayerPolys { + numerators, + denominators, + } = circuit.output_layer(); + let mut evaluations = numerators.evaluations().to_vec(); + evaluations.extend_from_slice(denominators.evaluations()); + transcript.reseed(H::hash_elements(&evaluations)); + + // generate the challenge and reduce [p0, p1, q0, q1] to [pr, qr] + let r = transcript.draw().map_err(|_| ProverError::FailedToGenerateChallenge)?; + let mut claim = circuit.evaluate_output_layer(r); + + let mut proof_layers: Vec> = Vec::new(); + let mut rand = vec![r]; + + // Loop over all inner layers, from output to input. + // + // In a layered circuit, each layer is defined in terms of its predecessor. The first inner + // layer (starting from the output layer) is the first layer that has a predecessor. Here, we + // loop over all inner layers in order to iteratively reduce a layer in terms of its successor + // layer. Note that we don't include the input layer, since its predecessor layer will be + // reduced in terms of the input layer separately in `prove_final_circuit_layer`. + for inner_layer in circuit.layers().iter().skip(1).rev().skip(1) { + // construct the Lagrange kernel evaluated at the previous GKR round randomness + let poly_x = EqFunction::ml_at(rand.clone()); + + // construct the vector of multi-linear polynomials + // TODO: avoid unnecessary allocation + let (left_numerators, right_numerators) = + inner_layer.numerators.project_least_significant_variable(); + let (left_denominators, right_denominators) = + inner_layer.denominators.project_least_significant_variable(); + let mls = + vec![left_numerators, right_numerators, left_denominators, right_denominators, poly_x]; + + // run the sumcheck protocol + let (proof, _) = sum_check_prover_plain_full(claim, mls, transcript)?; + + // sample a random challenge to reduce claims + transcript.reseed(H::hash_elements(&proof.openings_claim.openings)); + let r_layer = transcript.draw().map_err(|_| ProverError::FailedToGenerateChallenge)?; + + // reduce the claim + claim = { + let left_numerators_opening = proof.openings_claim.openings[0]; + let right_numerators_opening = proof.openings_claim.openings[1]; + let left_denominators_opening = proof.openings_claim.openings[2]; + let right_denominators_opening = proof.openings_claim.openings[3]; + + reduce_layer_claim( + left_numerators_opening, + right_numerators_opening, + left_denominators_opening, + right_denominators_opening, + r_layer, + ) + }; + + // collect the randomness used for the current layer + let mut ext = vec![r_layer]; + ext.extend_from_slice(&proof.openings_claim.eval_point); + rand = ext; + + proof_layers.push(proof); + } + + Ok(( + BeforeFinalLayerProof { + proof: proof_layers, + }, + GkrClaim { + evaluation_point: rand, + claimed_evaluation: claim, + }, + )) +} + +/// We receive our 4 multilinear polynomials which were evaluated at a random point: +/// `left_numerators` (or `p0`), `right_numerators` (or `p1`), `left_denominators` (or `q0`), and +/// `right_denominators` (or `q1`). We'll call the 4 evaluations at a random point `p0(r)`, `p1(r)`, +/// `q0(r)`, and `q1(r)`, respectively, where `r` is the random point. Note that `r` is a shorthand +/// for a tuple of random values `(r_0, ... r_{l-1})`, where `2^{l + 1}` is the number of wires in +/// the layer. +/// +/// It is important to recall how `p0` and `p1` were constructed (and analogously for `q0` and +/// `q1`). They are the `numerators` layer polynomial (or `p`) evaluations `p(0, r)` and `p(1, r)`, +/// obtained from [`MultiLinearPoly::project_least_significant_variable`]. Hence, `[p0, p1]` form +/// the evaluations of polynomial `p'(x_0) = p(x_0, r)`. Then, the round claim for `numerators`, +/// defined as `p(r_layer, r)`, is simply `p'(r_layer)`. +fn reduce_layer_claim( + left_numerators_opening: E, + right_numerators_opening: E, + left_denominators_opening: E, + right_denominators_opening: E, + r_layer: E, +) -> (E, E) +where + E: FieldElement, +{ + // This is the `numerators` layer polynomial `f(x_0) = numerators(x_0, rx_0, ..., rx_{l-1})`, + // where `rx_0, ..., rx_{l-1}` are the random variables that were sampled during the sumcheck + // round for this layer. + let numerators_univariate = + MultiLinearPoly::from_evaluations(vec![left_numerators_opening, right_numerators_opening]) + .unwrap(); + + // This is analogous to `numerators_univariate`, but for the `denominators` layer polynomial + let denominators_univariate = MultiLinearPoly::from_evaluations(vec![ + left_denominators_opening, + right_denominators_opening, + ]) + .unwrap(); + + ( + numerators_univariate.evaluate(&[r_layer]), + denominators_univariate.evaluate(&[r_layer]), + ) +} + +/// Runs the first sum-check prover for the input layer. +#[allow(clippy::type_complexity)] +fn sum_check_prover_plain_partial< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + claim: (E, E), + num_rounds: usize, + ml_polys: &mut [MultiLinearPoly], + transcript: &mut C, +) -> Result<((RoundClaim, Vec>), E), ProverError> { + // generate challenge to batch two sumchecks + let data = vec![claim.0, claim.1]; + transcript.reseed(H::hash_elements(&data)); + let r_batch = transcript.draw().map_err(|_| ProverError::FailedToGenerateChallenge)?; + let claim = claim.0 + claim.1 * r_batch; + + // generate the composition polynomial + let composer = GkrComposition::new(r_batch); + + // run the sum-check protocol + let main_prover = SumCheckProver::new(composer, SimpleGkrFinalClaimBuilder(PhantomData)); + let proof = main_prover.prove_rounds(claim, ml_polys, num_rounds, transcript)?; + + Ok((proof, r_batch)) +} + +/// Runs the sum-check prover used in all but the input layer. +fn sum_check_prover_plain_full< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + claim: (E, E), + ml_polys: Vec>, + transcript: &mut C, +) -> Result<(SumCheckProof, E), ProverError> { + // generate challenge to batch two sumchecks + transcript.reseed(H::hash_elements(&[claim.0, claim.1])); + let r_batch = transcript.draw().map_err(|_| ProverError::FailedToGenerateChallenge)?; + let claim_ = claim.0 + claim.1 * r_batch; + + // generate the composition polynomial + let composer = GkrComposition::new(r_batch); + + // run the sum-check protocol + let main_prover = SumCheckProver::new(composer, SimpleGkrFinalClaimBuilder(PhantomData)); + let proof = main_prover.prove(claim_, ml_polys, transcript)?; + + Ok((proof, r_batch)) +} + +/// Constructs [`FinalOpeningClaim`] for the sum-checks used in the GKR protocol. +/// +/// TODO: currently, this just removes the EQ evaluation as it can be computed by the verifier. +/// This should be generalized for other "transparent" multi-linears e.g., periodic columns. +struct SimpleGkrFinalClaimBuilder(PhantomData); + +impl FinalClaimBuilder for SimpleGkrFinalClaimBuilder { + type Field = E; + + fn build_claim( + &self, + openings: Vec, + evaluation_point: &[Self::Field], + ) -> FinalOpeningClaim { + FinalOpeningClaim { + eval_point: evaluation_point.to_vec(), + openings: (openings[..openings.len() - 1]).to_vec(), + } + } +} diff --git a/processor/src/trace/virtual_bus/circuit/verifier.rs b/processor/src/trace/virtual_bus/circuit/verifier.rs new file mode 100644 index 0000000000..9369f7067b --- /dev/null +++ b/processor/src/trace/virtual_bus/circuit/verifier.rs @@ -0,0 +1,237 @@ +use super::{ + error::VerifierError, prover::CircuitLayerPolys, FinalLayerProof, GkrCircuitProof, + GkrComposition, GkrCompositionMerge, +}; +use crate::trace::virtual_bus::{ + multilinear::EqFunction, + sum_check::{ + CompositionPolyQueryBuilder, FinalOpeningClaim, Proof as SumCheckFullProof, RoundClaim, + }, + SumCheckVerifier, +}; +use alloc::{borrow::ToOwned, vec::Vec}; +use vm_core::{Felt, FieldElement}; +use winter_prover::crypto::{ElementHasher, RandomCoin}; + +// VERIFIER +// ================================================================================================ + +/// Verifies the validity of a GKR proof for the correct evaluation of a fractional sum circuit. +pub fn verify< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + claim: E, + proof: GkrCircuitProof, + log_up_randomness: Vec, + transcript: &mut C, +) -> Result, VerifierError> { + let GkrCircuitProof { + circuit_outputs, + before_final_layer_proofs, + final_layer_proof, + } = proof; + + let CircuitLayerPolys { + numerators, + denominators, + } = circuit_outputs; + let p0 = numerators.evaluations()[0]; + let p1 = numerators.evaluations()[1]; + let q0 = denominators.evaluations()[0]; + let q1 = denominators.evaluations()[1]; + + // make sure that both denominators are not equal to E::ZERO + if q0 == E::ZERO || q1 == E::ZERO { + return Err(VerifierError::ZeroOutputDenominator); + } + + // check that the output matches the expected `claim` + if (p0 * q1 + p1 * q0) / (q0 * q1) != claim { + return Err(VerifierError::MismatchingCircuitOutput); + } + + // generate the random challenge to reduce two claims into a single claim + let mut evaluations = numerators.evaluations().to_vec(); + evaluations.extend_from_slice(denominators.evaluations()); + transcript.reseed(H::hash_elements(&evaluations)); + let r = transcript.draw().map_err(|_| VerifierError::FailedToGenerateChallenge)?; + + // reduce the claim + let p_r = p0 + r * (p1 - p0); + let q_r = q0 + r * (q1 - q0); + let mut reduced_claim = (p_r, q_r); + + // verify all GKR layers but for the last one + let num_layers = before_final_layer_proofs.proof.len(); + let mut rand = vec![r]; + for i in 0..num_layers { + let FinalOpeningClaim { + eval_point, + openings, + } = verify_sum_check_proof_before_last( + &before_final_layer_proofs.proof[i], + &rand, + reduced_claim, + transcript, + )?; + + // generate the random challenge to reduce two claims into a single claim + transcript.reseed(H::hash_elements(&openings)); + let r_layer = transcript.draw().map_err(|_| VerifierError::FailedToGenerateChallenge)?; + + let p0 = openings[0]; + let p1 = openings[1]; + let q0 = openings[2]; + let q1 = openings[3]; + reduced_claim = (p0 + r_layer * (p1 - p0), q0 + r_layer * (q1 - q0)); + + // collect the randomness used for the current layer + let rand_sumcheck = eval_point; + let mut ext = vec![r_layer]; + ext.extend_from_slice(&rand_sumcheck); + rand = ext; + } + + // verify the proof of the final GKR layer and pass final opening claim for verification + // to the STARK + verify_sum_check_proof_last( + final_layer_proof, + log_up_randomness, + &rand, + reduced_claim, + transcript, + ) +} + +// HELPER FUNCTIONS +// ================================================================================================ + +/// Verifies sum-check proofs, as part of the GKR proof, for all GKR layers except for the last one +/// i.e., the circuit input layer. +fn verify_sum_check_proof_before_last< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + proof: &SumCheckFullProof, + gkr_eval_point: &[E], + claim: (E, E), + transcript: &mut C, +) -> Result, VerifierError> { + // generate challenge to batch sum-checks + transcript.reseed(H::hash_elements(&[claim.0, claim.1])); + let r_batch: E = transcript.draw().map_err(|_| VerifierError::FailedToGenerateChallenge)?; + + // compute the claim for the batched sum-check + let reduced_claim = claim.0 + claim.1 * r_batch; + + // verify the sum-check protocol + let composition_poly = GkrComposition::new(r_batch); + let verifier = + SumCheckVerifier::new(composition_poly, GkrQueryBuilder::new(gkr_eval_point.to_owned())); + verifier + .verify(reduced_claim, proof.clone(), transcript) + .map_err(VerifierError::FailedToVerifySumCheck) +} + +/// Verifies the final sum-check proof as part of the GKR proof. +fn verify_sum_check_proof_last< + E: FieldElement, + C: RandomCoin, + H: ElementHasher, +>( + proof: FinalLayerProof, + log_up_randomness: Vec, + gkr_eval_point: &[E], + claim: (E, E), + transcript: &mut C, +) -> Result, VerifierError> { + let FinalLayerProof { + before_merge_proof, + after_merge_proof, + } = proof; + + // generate challenge to batch sum-checks + transcript.reseed(H::hash_elements(&[claim.0, claim.1])); + let r_sum_check: E = transcript.draw().map_err(|_| VerifierError::FailedToGenerateChallenge)?; + + // compute the claim for the batched sum-check + let reduced_claim = claim.0 + claim.1 * r_sum_check; + + // verify the first part of the sum-check protocol + let composition_poly = GkrComposition::new(r_sum_check); + let verifier = + SumCheckVerifier::new(composition_poly, GkrQueryBuilder::new(gkr_eval_point.to_owned())); + let RoundClaim { + eval_point: rand_merge, + claim, + } = verifier.verify_rounds(reduced_claim, before_merge_proof, transcript)?; + + // verify the second part of the sum-check protocol + let gkr_composition = + GkrCompositionMerge::new(r_sum_check, rand_merge.clone(), log_up_randomness); + let verifier = SumCheckVerifier::new( + gkr_composition, + GkrMergeQueryBuilder::new(gkr_eval_point.to_vec(), rand_merge), + ); + verifier + .verify(claim, after_merge_proof, transcript) + .map_err(VerifierError::FailedToVerifySumCheck) +} + +// GKR QUERY BUILDERS +// ================================================================================================ + +/// A [`GkrQueryBuilder`] for the sum-check verifier used for all sum-checks but for the final one. +#[derive(Default)] +struct GkrQueryBuilder { + gkr_eval_point: Vec, +} + +impl GkrQueryBuilder { + fn new(gkr_eval_point: Vec) -> Self { + Self { gkr_eval_point } + } +} + +impl CompositionPolyQueryBuilder for GkrQueryBuilder { + fn build_query(&self, openings_claim: &FinalOpeningClaim, evaluation_point: &[E]) -> Vec { + let rand_sumcheck = evaluation_point; + let eq_at_gkr_eval_point = EqFunction::new(self.gkr_eval_point.clone()); + let eq = eq_at_gkr_eval_point.evaluate(rand_sumcheck); + + let mut query = openings_claim.openings.clone(); + query.push(eq); + query + } +} + +/// A [`GkrMergeQueryBuilder`] for the sum-check verifier used for the final sum-check. +#[derive(Default)] +struct GkrMergeQueryBuilder { + gkr_eval_point: Vec, + merge_rand: Vec, +} + +impl GkrMergeQueryBuilder { + fn new(gkr_eval_point: Vec, merge_rand: Vec) -> Self { + Self { + gkr_eval_point, + merge_rand, + } + } +} + +impl CompositionPolyQueryBuilder for GkrMergeQueryBuilder { + fn build_query(&self, openings_claim: &FinalOpeningClaim, evaluation_point: &[E]) -> Vec { + let eq_at_gkr_eval_point = EqFunction::new(self.gkr_eval_point.clone()); + let mut rand_sumcheck = self.merge_rand.clone(); + rand_sumcheck.extend_from_slice(evaluation_point); + let eq = eq_at_gkr_eval_point.evaluate(&rand_sumcheck); + let mut query = openings_claim.openings.clone(); + query.push(eq); + query + } +} diff --git a/processor/src/trace/virtual_bus/mod.rs b/processor/src/trace/virtual_bus/mod.rs new file mode 100644 index 0000000000..02388d46bb --- /dev/null +++ b/processor/src/trace/virtual_bus/mod.rs @@ -0,0 +1,39 @@ +//! Represents a global virtual bus for the Miden VM. +//! +//! A global bus is a single bus which encompasses several sub-buses, each representing +//! a communication channel between two or more components of the VM. +//! +//! A bus represents a client-server relationship between some VM components. The server is +//! usually one specific component of the VM, e.g., hasher chiplet, and the client can be one or +//! several other components of the VM, e.g., the decoder. The communication between the clients +//! and the server is composed of `request` messages made by the clients and corresponding +//! `reply` messages by the server. +//! +//! The purpose of the sub-bus then, from the verifiable computation point of view, is to ensure +//! the consistency between the `request` and `reply` messages exchanged by the clients and +//! the server. +//! +//! The global bus uses a per-sub-bus address in order to ensure correct routing of the `request` +//! messages and their matching `reply` messages. +//! +//! Miden VM uses a virtual global bus in the sense that neither the global bus nor the individual +//! sub-buses are fully materialized as part of the (auxiliary) trace. This is replaced by a layered +//! circuit which computes the global bus relation. The correct evaluation of this circuit is then +//! proved using the GKR protocol of Goldwasser, Kalai and Rothblum [1] using the protocol in +//! GKR-LogUp [2]. +//! +//! [1]: https://dl.acm.org/doi/10.1145/2699436 +//! [2]: https://eprint.iacr.org/2023/1284 + +mod circuit; +pub use circuit::{prove, verify}; + +mod multilinear; + +mod sum_check; +pub use sum_check::{SumCheckProver, SumCheckVerifier}; + +mod univariate; + +#[cfg(test)] +mod tests; diff --git a/processor/src/trace/virtual_bus/multilinear/error.rs b/processor/src/trace/virtual_bus/multilinear/error.rs new file mode 100644 index 0000000000..cd4f3c5789 --- /dev/null +++ b/processor/src/trace/virtual_bus/multilinear/error.rs @@ -0,0 +1,5 @@ +#[derive(Debug, thiserror::Error)] +pub enum Error { + #[error("A multi-linear polynomial should have a power of 2 number of evaluations over the Boolean hyper-cube")] + EvaluationsNotPowerOfTwo, +} diff --git a/processor/src/trace/virtual_bus/multilinear/lagrange_ker.rs b/processor/src/trace/virtual_bus/multilinear/lagrange_ker.rs new file mode 100644 index 0000000000..4c309a6a56 --- /dev/null +++ b/processor/src/trace/virtual_bus/multilinear/lagrange_ker.rs @@ -0,0 +1,103 @@ +use super::{FieldElement, MultiLinearPoly}; +use alloc::vec::Vec; + +/// The EQ (equality) function is the binary function defined by +/// +/// EQ: {0 , 1}^ν ⛌ {0 , 1}^ν ⇾ {0 , 1} +/// ((x_0, ..., x_{ν - 1}), (y_0, ..., y_{ν - 1})) ↦ \prod_{i = 0}^{ν - 1} (x_i * y_i + (1 - x_i) +/// * (1 - y_i)) +/// +/// Taking It's multi-linear extension EQ^{~}, we can define a basis for the set of multi-linear +/// polynomials in ν variables by +/// {EQ^{~}(., (y_0, ..., y_{ν - 1})): (y_0, ..., y_{ν - 1}) ∈ {0 , 1}^ν} +/// where each basis function is a function of its first argument. This is called the Lagrange or +/// evaluation basis with evaluation set {0 , 1}^ν. +/// +/// Given a function (f: {0 , 1}^ν ⇾ 𝔽), its multi-linear extension (i.e., the unique +/// mult-linear polynomial extending f to (f^{~}: 𝔽^ν ⇾ 𝔽) and agrees with it on {0 , 1}^ν) is +/// defined as the summation of the evaluations of f against the Lagrange basis. +/// More specifically, given (r_0, ..., r_{ν - 1}) ∈ 𝔽^ν, then: +/// +/// f^{~}(r_0, ..., r_{ν - 1}) = \sum_{(y_0, ..., y_{ν - 1}) ∈ {0 , 1}^ν} +/// f(y_0, ..., y_{ν - 1}) EQ^{~}((r_0, ..., r_{ν - 1}), (y_0, ..., y_{ν - 1})) +/// +/// We call the Lagrange kernel the evaluation of the EQ^{~} function at +/// ((r_0, ..., r_{ν - 1}), (y_0, ..., y_{ν - 1})) for all (y_0, ..., y_{ν - 1}) ∈ {0 , 1}^ν for +/// a fixed (r_0, ..., r_{ν - 1}) ∈ 𝔽^ν. +/// +/// [`EqFunction`] represents EQ^{~} the mult-linear extension of +/// +/// ((y_0, ..., y_{ν - 1}) ↦ EQ((r_0, ..., r_{ν - 1}), (y_0, ..., y_{ν - 1}))) +/// +/// and contains a method to generate the Lagrange kernel for defining evaluations of multi-linear +/// extensions of arbitrary functions (f: {0 , 1}^ν ⇾ 𝔽) at a given point (r_0, ..., r_{ν - 1}) +/// as well as a method to evaluate EQ^{~}((r_0, ..., r_{ν - 1}), (t_0, ..., t_{ν - 1}))) for +/// (t_0, ..., t_{ν - 1}) ∈ 𝔽^ν. +pub struct EqFunction { + r: Vec, +} + +impl EqFunction { + /// Creates a new [EqFunction]. + pub fn new(r: Vec) -> Self { + let tmp = r.clone(); + EqFunction { r: tmp } + } + + /// Computes EQ((r_0, ..., r_{ν - 1}), (t_0, ..., t_{ν - 1}))). + pub fn evaluate(&self, t: &[E]) -> E { + assert_eq!(self.r.len(), t.len()); + + (0..self.r.len()) + .map(|i| self.r[i] * t[i] + (E::ONE - self.r[i]) * (E::ONE - t[i])) + .fold(E::ONE, |acc, term| acc * term) + } + + /// Computes EQ((r_0, ..., r_{ν - 1}), (y_0, ..., y_{ν - 1})) for all + /// (y_0, ..., y_{ν - 1}) ∈ {0 , 1}^ν i.e., the Lagrange kernel at r = (r_0, ..., r_{ν - 1}). + pub fn evaluations(&self) -> Vec { + compute_lagrange_basis_evals_at(&self.r) + } + + /// Returns the evaluations of + /// ((y_0, ..., y_{ν - 1}) ↦ EQ^{~}((r_0, ..., r_{ν - 1}), (y_0, ..., y_{ν - 1}))) + /// over {0 , 1}^ν. + pub fn ml_at(evaluation_point: Vec) -> MultiLinearPoly { + let eq_evals = EqFunction::new(evaluation_point.clone()).evaluations(); + MultiLinearPoly::from_evaluations(eq_evals) + .expect("should not fail because evaluations is a power of two") + } +} + +// HELPER +// ================================================================================================ + +/// Computes the inner product of two vectors of the same length. +/// +/// Panics if the vectors have different lengths. +pub fn inner_product(evaluations: &[E], tensored_query: &[E]) -> E { + assert_eq!(evaluations.len(), tensored_query.len()); + evaluations + .iter() + .zip(tensored_query.iter()) + .fold(E::ZERO, |acc, (x_i, y_i)| acc + *x_i * *y_i) +} + +/// Computes the evaluations of the Lagrange basis polynomials over the interpolating +/// set {0 , 1}^ν at (r_0, ..., r_{ν - 1}) i.e., the Lagrange kernel at (r_0, ..., r_{ν - 1}). +pub fn compute_lagrange_basis_evals_at(query: &[E]) -> Vec { + let nu = query.len(); + let n = 1 << nu; + + let mut evals: Vec = vec![E::ONE; n]; + let mut size = 1; + for r_i in query.iter().rev() { + size *= 2; + for i in (0..size).rev().step_by(2) { + let scalar = evals[i / 2]; + evals[i] = scalar * *r_i; + evals[i - 1] = scalar - evals[i]; + } + } + evals +} diff --git a/processor/src/trace/virtual_bus/multilinear/mod.rs b/processor/src/trace/virtual_bus/multilinear/mod.rs new file mode 100644 index 0000000000..de892f22d2 --- /dev/null +++ b/processor/src/trace/virtual_bus/multilinear/mod.rs @@ -0,0 +1,108 @@ +use alloc::vec::Vec; +use core::ops::Index; +use vm_core::FieldElement; + +mod lagrange_ker; +pub use lagrange_ker::{inner_product, EqFunction}; + +mod error; +use self::{error::Error, lagrange_ker::compute_lagrange_basis_evals_at}; + +// MULTI-LINEAR POLYNOMIAL +// ================================================================================================ + +/// Represents a multi-linear polynomial. +/// +/// The representation stores the evaluations of the polynomial over the boolean hyper-cube +/// {0 , 1}^ν. +#[derive(Clone, Debug)] +pub struct MultiLinearPoly { + num_variables: usize, + evaluations: Vec, +} + +impl MultiLinearPoly { + /// Constructs a [`MultiLinearPoly`] from its evaluations over the boolean hyper-cube {0 , 1}^ν. + pub fn from_evaluations(evaluations: Vec) -> Result { + if !evaluations.len().is_power_of_two() { + return Err(Error::EvaluationsNotPowerOfTwo); + } + Ok(Self { + num_variables: (evaluations.len().ilog2()) as usize, + evaluations, + }) + } + + /// Returns the number of variables of the multi-linear polynomial. + pub fn num_variables(&self) -> usize { + self.num_variables + } + + /// Returns the evaluations over the boolean hyper-cube. + pub fn evaluations(&self) -> &[E] { + &self.evaluations + } + + /// Returns the number of evaluations. This is equal to the size of the boolean hyper-cube. + pub fn num_evaluations(&self) -> usize { + self.evaluations.len() + } + + /// Evaluate the multi-linear at some query (r_0, ..., r_{ν - 1}) ∈ 𝔽^ν. + /// + /// It first computes the evaluations of the Lagrange basis polynomials over the interpolating + /// set {0 , 1}^ν at (r_0, ..., r_{ν - 1}) i.e., the Lagrange kernel at (r_0, ..., r_{ν - 1}). + /// The evaluation then is the inner product, indexed by {0 , 1}^ν, of the vector of + /// evaluations times the Lagrange kernel. + pub fn evaluate(&self, query: &[E]) -> E { + let tensored_query = compute_lagrange_basis_evals_at(query); + inner_product(&self.evaluations, &tensored_query) + } + + /// Similar to [`Self::evaluate`], except that the query was already turned into the Lagrange + /// kernel (i.e. the [`lagrange_ker::EqFunction`] evaluated at every point in the set + /// `{0 , 1}^ν`). + /// + /// This is more efficient than [`Self::evaluate`] when multiple different [`MultiLinearPoly`] + /// need to be evaluated at the same query point. + pub fn evaluate_with_lagrange_kernel(&self, lagrange_kernel: &[E]) -> E { + inner_product(&self.evaluations, lagrange_kernel) + } + + /// Computes f(r_0, y_1, ..., y_{ν - 1}) using the linear interpolation formula + /// (1 - r_0) * f(0, y_1, ..., y_{ν - 1}) + r_0 * f(1, y_1, ..., y_{ν - 1}) and assigns + /// the resulting multi-linear, defined over a domain of half the size, to `self`. + pub fn bind_least_significant_variable(&mut self, round_challenge: E) { + let mut result = vec![E::ZERO; 1 << (self.num_variables() - 1)]; + for (i, res) in result.iter_mut().enumerate() { + *res = self.evaluations[i << 1] + + round_challenge * (self.evaluations[(i << 1) + 1] - self.evaluations[i << 1]); + } + *self = Self::from_evaluations(result) + .expect("should not fail given that it is a multi-linear"); + } + + /// Given the multilinear polynomial f(y_0, y_1, ..., y_{ν - 1}), returns two polynomials: + /// f(0, y_1, ..., y_{ν - 1}) and f(1, y_1, ..., y_{ν - 1}). + pub fn project_least_significant_variable(&self) -> (Self, Self) { + let mut p0 = Vec::with_capacity(self.num_evaluations() / 2); + let mut p1 = Vec::with_capacity(self.num_evaluations() / 2); + for chunk in self.evaluations.chunks_exact(2) { + p0.push(chunk[0]); + p1.push(chunk[1]); + } + + ( + MultiLinearPoly::from_evaluations(p0).unwrap(), + MultiLinearPoly::from_evaluations(p1).unwrap(), + ) + } +} + +impl Index for MultiLinearPoly { + type Output = E; + + fn index(&self, index: usize) -> &E { + &(self.evaluations[index]) + } +} diff --git a/processor/src/trace/virtual_bus/sum_check/mod.rs b/processor/src/trace/virtual_bus/sum_check/mod.rs new file mode 100644 index 0000000000..0565bc4c45 --- /dev/null +++ b/processor/src/trace/virtual_bus/sum_check/mod.rs @@ -0,0 +1,250 @@ +use super::univariate::UnivariatePolyCoef; +use alloc::vec::Vec; +use vm_core::FieldElement; + +mod prover; +pub use prover::{Error as SumCheckProverError, FinalClaimBuilder, SumCheckProver}; +mod verifier; +pub use verifier::{CompositionPolyQueryBuilder, Error as SumCheckVerifierError, SumCheckVerifier}; + +/// A sum-check round proof. +/// +/// This represents the partial polynomial sent by the Prover during one of the rounds of the +/// sum-check protocol. The polynomial is in coefficient form and excludes the coefficient for +/// the linear term as the Verifier can recover it from the other coefficients and the current +/// (reduced) claim. +#[derive(Debug, Clone)] +pub struct RoundProof { + pub round_poly_coefs: UnivariatePolyCoef, +} + +/// A sum-check proof. +/// +/// Composed of the round proofs i.e., the polynomials sent by the Prover at each round as well as +/// the (claimed) openings of the multi-linear oracles at the evaluation point given by the round +/// challenges. +#[derive(Debug, Clone)] +pub struct Proof { + pub openings_claim: FinalOpeningClaim, + pub round_proofs: Vec>, +} + +/// Contains the round challenges sent by the Verifier up to some round as well as the current +/// reduced claim. +#[derive(Debug)] +pub struct RoundClaim { + pub eval_point: Vec, + pub claim: E, +} + +/// Reduces an old claim to a new claim using the round challenge. +pub fn reduce_claim( + current_poly: &RoundProof, + current_round_claim: RoundClaim, + round_challenge: E, +) -> RoundClaim { + // evaluate the round polynomial at the round challenge to obtain the new claim + let new_claim = current_poly + .round_poly_coefs + .evaluate_using_claim(¤t_round_claim.claim, &round_challenge); + + // update the evaluation point using the round challenge + let mut new_partial_eval_point = current_round_claim.eval_point; + new_partial_eval_point.push(round_challenge); + + RoundClaim { + eval_point: new_partial_eval_point, + claim: new_claim, + } +} + +/// Represents an opening claim at an evaluation point against a batch of oracles. +/// +/// After verifying [`Proof`], the verifier is left with a question on the validity of a final +/// claim on a number of oracles open to a given set of values at some given point. +/// This question is answered either using further interaction with the Prover or using +/// a polynomial commitment opening proof in the compiled protocol. +#[derive(Clone, Debug)] +pub struct FinalOpeningClaim { + pub eval_point: Vec, + pub openings: Vec, +} + +// COMPOSITION POLYNOMIAL +// ================================================================================================ + +/// A multi-variate polynomial for composing individual multi-linear polynomials. +pub trait CompositionPolynomial { + /// Maximum degree in all variables. + fn max_degree(&self) -> u32; + + /// Given a query, of length equal the number of variables, evaluates [Self] at this query. + fn evaluate(&self, query: &[E]) -> E; +} + +// TESTS +// ================================================================================================ + +#[cfg(test)] +mod tests { + use super::{ + prover::{FinalClaimBuilder, SumCheckProver}, + verifier::{CompositionPolyQueryBuilder, SumCheckVerifier}, + CompositionPolynomial, + }; + use crate::trace::virtual_bus::multilinear::MultiLinearPoly; + use alloc::{borrow::ToOwned, vec::Vec}; + use test_utils::rand::rand_vector; + use vm_core::{crypto::random::RpoRandomCoin, Felt, FieldElement, Word, ONE, ZERO}; + + #[test] + fn test_sum_check_sum() { + let num_variables = 14; + let values = rand_vector(1 << num_variables); + let claim = values.iter().fold(ZERO, |acc, &x| x + acc); + + let ml = MultiLinearPoly::from_evaluations(values.to_vec()).expect("should not fail"); + let mls = vec![ml]; + let virtual_poly = ProjectionComposition::new(0); + + // Prover + let prover = SumCheckProver::new(virtual_poly, PlainClaimBuilder); + let mut coin = RpoRandomCoin::new(Word::default()); + let proof = prover.prove(claim, mls, &mut coin).unwrap(); + + // Verifier + let plain_query_builder = ProjectionPolyQueryBuilder::default(); + let verifier = SumCheckVerifier::new(virtual_poly, plain_query_builder); + let mut coin = RpoRandomCoin::new(Word::default()); + let result = verifier.verify(claim, proof, &mut coin); + + assert!(result.is_ok()) + } + + #[test] + fn test_sum_check_product() { + let num_variables = 14; + let values_0 = rand_vector(1 << num_variables); + let values_1 = rand_vector(1 << num_variables); + let claim = values_0.iter().zip(values_1.iter()).fold(ZERO, |acc, (x, y)| *x * *y + acc); + + let ml_0 = MultiLinearPoly::from_evaluations(values_0.to_vec()).expect("should not fail"); + let ml_1 = MultiLinearPoly::from_evaluations(values_1.to_vec()).expect("should not fail"); + let mls = vec![ml_0, ml_1]; + let virtual_poly = ProductComposition; + + // Prover + let prover = SumCheckProver::new(virtual_poly, PlainClaimBuilder); + let mut coin = RpoRandomCoin::new(Word::default()); + let proof = prover.prove(claim, mls, &mut coin).unwrap(); + + // Verifier + let plain_query_builder = ProjectionPolyQueryBuilder::default(); + let verifier = SumCheckVerifier::new(virtual_poly, plain_query_builder); + let mut coin = RpoRandomCoin::new(Word::default()); + let result = verifier.verify(claim, proof, &mut coin); + + assert!(result.is_ok()) + } + + #[test] + fn test_sum_check_product_failure() { + let num_variables = 14; + let values_0 = rand_vector(1 << num_variables); + let values_1 = rand_vector(1 << num_variables); + let mut claim = + values_0.iter().zip(values_1.iter()).fold(ZERO, |acc, (x, y)| *x * *y + acc); + + // modifying the claim should make the Verifier reject the proof + claim += ONE; + + let ml_0 = MultiLinearPoly::from_evaluations(values_0.to_vec()).expect("should not fail"); + let ml_1 = MultiLinearPoly::from_evaluations(values_1.to_vec()).expect("should not fail"); + let mls = vec![ml_0, ml_1]; + let virtual_poly = ProductComposition; + + // Prover + let prover = SumCheckProver::new(virtual_poly, PlainClaimBuilder); + let mut coin = RpoRandomCoin::new(Word::default()); + let proof = prover.prove(claim, mls, &mut coin).unwrap(); + + // Verifier + let plain_query_builder = ProjectionPolyQueryBuilder::default(); + let verifier = SumCheckVerifier::new(virtual_poly, plain_query_builder); + let mut coin = RpoRandomCoin::new(Word::default()); + let result = verifier.verify(claim, proof, &mut coin); + + assert!(result.is_err()) + } + + struct PlainClaimBuilder; + + impl FinalClaimBuilder for PlainClaimBuilder { + type Field = Felt; + + fn build_claim( + &self, + openings: Vec, + evaluation_point: &[Self::Field], + ) -> super::FinalOpeningClaim { + super::FinalOpeningClaim { + eval_point: evaluation_point.to_owned(), + openings, + } + } + } + + #[derive(Default)] + struct ProjectionPolyQueryBuilder; + + impl CompositionPolyQueryBuilder for ProjectionPolyQueryBuilder { + fn build_query( + &self, + openings_claim: &super::FinalOpeningClaim, + _evaluation_point: &[E], + ) -> Vec { + openings_claim.openings.to_vec() + } + } + + #[derive(Clone, Copy, Debug)] + pub struct ProjectionComposition { + coordinate: usize, + } + + impl ProjectionComposition { + pub fn new(coordinate: usize) -> Self { + Self { coordinate } + } + } + + impl CompositionPolynomial for ProjectionComposition + where + E: FieldElement, + { + fn max_degree(&self) -> u32 { + 1 + } + + fn evaluate(&self, query: &[E]) -> E { + query[self.coordinate] + } + } + + #[derive(Clone, Copy, Debug, Default)] + pub struct ProductComposition; + + impl CompositionPolynomial for ProductComposition + where + E: FieldElement, + { + fn max_degree(&self) -> u32 { + 2 + } + + fn evaluate(&self, query: &[E]) -> E { + assert_eq!(query.len(), 2); + query[0] * query[1] + } + } +} diff --git a/processor/src/trace/virtual_bus/sum_check/prover/error.rs b/processor/src/trace/virtual_bus/sum_check/prover/error.rs new file mode 100644 index 0000000000..e425bbfbfd --- /dev/null +++ b/processor/src/trace/virtual_bus/sum_check/prover/error.rs @@ -0,0 +1,15 @@ +#[derive(Debug, thiserror::Error)] +pub enum Error { + #[error("number of rounds for sum-check must be greater than zero")] + NumRoundsZero, + #[error("the number of rounds is greater than the number of variables")] + TooManyRounds, + #[error("should provide at least one multi-linear polynomial as input")] + NoMlsProvided, + #[error("failed to generate round challenge")] + FailedToGenerateChallenge, + #[error("the provided multi-linears have different arities")] + MlesDifferentArities, + #[error("multi-linears should have at least one variable")] + AtLeastOneVariable, +} diff --git a/processor/src/trace/virtual_bus/sum_check/prover/mod.rs b/processor/src/trace/virtual_bus/sum_check/prover/mod.rs new file mode 100644 index 0000000000..ce211ed686 --- /dev/null +++ b/processor/src/trace/virtual_bus/sum_check/prover/mod.rs @@ -0,0 +1,329 @@ +use super::{ + reduce_claim, CompositionPolynomial, FinalOpeningClaim, Proof, RoundClaim, RoundProof, +}; +use crate::trace::virtual_bus::{multilinear::MultiLinearPoly, univariate::UnivariatePolyEvals}; +use alloc::vec::Vec; +use core::marker::PhantomData; +use vm_core::FieldElement; +use winter_prover::crypto::{ElementHasher, RandomCoin}; + +mod error; +pub use self::error::Error; + +/// A struct that contains relevant information for the execution of the multivariate sum-check +/// protocol prover. +/// The sum-check protocol is an interactive protocol (IP) for proving the following relation: +/// +/// v = \sum_{(x_0,\cdots, x_{\nu - 1}) \in \{0 , 1\}^{\nu}} +/// g(f_0((x_0,\cdots, x_{\nu - 1})), \cdots , f_c((x_0,\cdots, x_{\nu - 1}))) +/// +/// where: +/// +/// 1. v ∈ 𝔽 where 𝔽 is a finite field. +/// 2. f_i are multi-linear polynomials i.e., polynomials in 𝔽[X_0, \cdots ,X_{\nu - 1}] with degree +/// at most one in each variable. +/// 3. g is a multivariate polynomial with degree at most d in each variable. +/// +/// The Verifier is given commitments to each `f_i` in addition to the claimed sum `v`. The Prover +/// then engages in an IP to convince the Verifier that the above relation holds for the given +/// `f_i` and `v`. More precisely: +/// +/// 0. Denote by w(x_0,\cdots, x_{\nu - 1}) := g(f_0((x_0,\cdots, x_{\nu - 1})), \cdots , +/// f_c((x_0,\cdots, x_{\nu - 1}))). +/// +/// 1. In the first round, the Prover sends the polynomial defined by: s_0(X_0) := +/// \sum_{(x_{1},\cdots, x_{\nu - 1}) w(X_0, x_{1}, \cdots, x_{\nu - 1}) +/// +/// 2. The Verifier then checks that s_0(0) + s_0(1) = v rejecting if not. +/// +/// 3. The Verifier samples a random challenge `r_0 ∈ 𝔽` and sends it to the Prover. +/// +/// 4. For each i in 1...(\nu - 1): a. The Prover sends the univariate polynomial defined by: +/// +/// s_i(X_i) := \sum_{(x_{i + 1},\cdots, x_{\nu - 1}) +/// w(r_0,\cdots, r_{i - 1}, X_i, x_{i + 1}, \cdots, x_{\nu - 1}). +/// +/// b. The Verifier checks that s_{i - 1}(r_{i - 1}) = s_{i}(0) + s_{i}(1) rejecting if not. +/// +/// c. The Verifier samples a random challenge `r_i ∈ 𝔽` and sends it to the Prover. +/// +/// 5. The Verifier now queries each of the oracles behind the commitments i.e., `f_i` at +/// `(r_0, \cdots , r_{\nu - 1})` to get u_i = f_i(r_0, \cdots , r_{\nu - 1}). +/// The Verifier then accepts if and only if: +/// +/// s_{\nu - 1}(r_{\nu - 1}) = g(u_0, \cdots , u_{\nu - 1}) +/// +/// A few remarks: +/// +/// 1. The degree bound on `g` implies that each of the `s_i` polynomials is a univariate polynomial +/// of degree at most `d`. Thus, the Prover in each round sends `d + 1` values, either +/// the coefficients or the evaluations of `s_i`. +/// +/// 2. The Prover has each `f_i` in its evaluation form over the hyper-cube \{0 , 1\}^{\nu}. +/// +/// 3. An optimization is for the Prover to not send `s_i(0)` as it can be recovered from the +/// current +/// reduced claim s_{i - 1}(r_{i - 1}) using the relation s_{i}(0) = s_{i}(1) - s_{i - 1}(r_{i - +/// 1}). This also means that the Verifier can skip point 4.b. +pub struct SumCheckProver +where + E: FieldElement, + C: RandomCoin, + H: ElementHasher, + V: FinalClaimBuilder, +{ + composition_poly: P, + final_claim_builder: V, + + _challenger: PhantomData, +} + +impl SumCheckProver +where + E: FieldElement, + P: CompositionPolynomial, + C: RandomCoin, + H: ElementHasher, + V: FinalClaimBuilder, +{ + /// Constructs a new [SumCheckProver] given a multivariate composition polynomial. + /// The multivariate composition polynomial corresponds to the `g` polynomial in the + /// description of the [SumCheckProver] struct. + pub fn new(composition_poly: P, final_claim_builder: V) -> Self { + Self { + composition_poly, + final_claim_builder, + _challenger: PhantomData, + } + } + + /// Given an initial claim `claim`, a mutable vector of multi-linear polynomials `mls` and + /// a number of rounds `num_rounds`, computes `num_rounds` iterations of the sum-check protocol + /// starting from claim `claim`. + /// + /// More specifically, executes the sum-check protocol for the following relation + /// + /// v = \sum_{(x_0,\cdots, x_{\nu - 1}) \in \{0 , 1\}^{\nu}} + /// g(f_0((x_0,\cdots, x_{\nu - 1})), \cdots , f_c((x_0,\cdots, x_{\nu - + /// 1}))) + /// + /// where: + /// + /// 1. `claim` is v. + /// 2. `mls` is [f_0, ..., f_c]. + /// 3. `self.composition_poly` is g. + /// + /// # Errors + /// Returns an error if: + /// - No multi-linears were provided. + /// - Number of rounds is zero or is greater than the number of variables of the multilinears. + /// - The provided multi-linears have different arities. + pub fn prove( + &self, + claim: E, + mut mls: Vec>, + coin: &mut C, + ) -> Result, Error> { + let num_rounds = mls[0].num_variables(); + let ( + RoundClaim { + eval_point, + claim: _claim, + }, + round_proofs, + ) = self.prove_rounds(claim, &mut mls, num_rounds, coin)?; + + let openings = mls.iter_mut().map(|ml| ml.evaluations()[0]).collect(); + let openings_claim = self.final_claim_builder.build_claim(openings, &eval_point); + + Ok(Proof { + openings_claim, + round_proofs, + }) + } + + /// Proves a round of the sum-check protocol. + pub fn prove_rounds( + &self, + claim: E, + mls: &mut [MultiLinearPoly], + num_rounds: usize, + coin: &mut C, + ) -> Result<(RoundClaim, Vec>), Error> { + // there should be at least one multi-linear polynomial provided + if mls.is_empty() { + return Err(Error::NoMlsProvided); + } + + // there should be at least one round to prove + if num_rounds == 0 { + return Err(Error::NumRoundsZero); + } + + // there can not be more rounds than variables of the multi-linears + let ml_variables = mls[0].num_variables(); + if num_rounds > ml_variables { + return Err(Error::TooManyRounds); + } + + // there should at least be one variable for the protocol to be non-trivial + if ml_variables < 1 { + return Err(Error::AtLeastOneVariable); + } + + // all multi-linears should have the same arity + if !mls.iter().all(|ml| ml.num_variables() == ml_variables) { + return Err(Error::MlesDifferentArities); + } + + let mut round_proofs = vec![]; + + // setup first round claim + let mut current_round_claim = RoundClaim { + eval_point: vec![], + claim, + }; + + // run the first round of the protocol + let round_poly_evals = sumcheck_round(&self.composition_poly, mls); + let round_poly_coefs = round_poly_evals.to_poly(current_round_claim.claim); + + // reseed with the s_0 polynomial + coin.reseed(H::hash_elements(&round_poly_coefs.coefficients)); + round_proofs.push(RoundProof { round_poly_coefs }); + + for i in 1..num_rounds { + // generate random challenge r_i for the i-th round + let round_challenge = coin.draw().map_err(|_| Error::FailedToGenerateChallenge)?; + + // compute the new reduced round claim + let new_round_claim = + reduce_claim(&round_proofs[i - 1], current_round_claim, round_challenge); + + // fold each multi-linear using the round challenge + mls.iter_mut() + .for_each(|ml| ml.bind_least_significant_variable(round_challenge)); + + // run the i-th round of the protocol using the folded multi-linears for the new reduced + // claim. This basically computes the s_i polynomial. + let round_poly_evals = sumcheck_round(&self.composition_poly, mls); + + // update the claim + current_round_claim = new_round_claim; + + let round_poly_coefs = round_poly_evals.to_poly(current_round_claim.claim); + + // reseed with the s_i polynomial + coin.reseed(H::hash_elements(&round_poly_coefs.coefficients)); + let round_proof = RoundProof { round_poly_coefs }; + round_proofs.push(round_proof); + } + + // generate the last random challenge + let round_challenge = coin.draw().map_err(|_| Error::FailedToGenerateChallenge)?; + // fold each multi-linear using the last random challenge + mls.iter_mut() + .for_each(|ml| ml.bind_least_significant_variable(round_challenge)); + + let round_claim = + reduce_claim(&round_proofs[num_rounds - 1], current_round_claim, round_challenge); + Ok((round_claim, round_proofs)) + } +} + +/// Computes the polynomial +/// +/// s_i(X_i) := \sum_{(x_{i + 1},\cdots, x_{\nu - 1}) +/// w(r_0,\cdots, r_{i - 1}, X_i, x_{i + 1}, \cdots, x_{\nu - 1}). +/// where +/// +/// w(x_0,\cdots, x_{\nu - 1}) := g(f_0((x_0,\cdots, x_{\nu - 1})), +/// \cdots , f_c((x_0,\cdots, x_{\nu - 1}))). +/// +/// Given a degree bound `d_max` for all variables, it suffices to compute the evaluations of `s_i` +/// at `d_max + 1` points. Given that `s_{i}(0) = s_{i}(1) - s_{i - 1}(r_{i - 1})` it is sufficient +/// to compute the evaluations on only `d_max` points. +/// +/// The algorithm works by iterating over the variables (x_{i + 1}, \cdots, x_{\nu - 1}) in +/// {0, 1}^{\nu - 1 - i}. For each such tuple, we store the evaluations of the (folded) +/// multi-linears at (0, x_{i + 1}, \cdots, x_{\nu - 1}) and +/// (1, x_{i + 1}, \cdots, x_{\nu - 1}) in two arrays, `evals_zero` and `evals_one`. +/// Using `evals_one`, remember that we optimize evaluating at 0 away, we get the first evaluation +/// i.e., `s_i(1)`. +/// +/// For the remaining evaluations, we use the fact that the folded `f_i` is multi-linear and hence +/// we can write +/// +/// f_i(X_i, x_{i + 1}, \cdots, x_{\nu - 1}) = +/// (1 - X_i) . f_i(0, x_{i + 1}, \cdots, x_{\nu - 1}) + X_i . f_i(1, x_{i + 1}, \cdots, +/// x_{\nu - 1}) +/// +/// Note that we omitted writing the folding randomness for readability. +/// Since the evaluation domain is {0, 1, ... , d_max}, we can compute the evaluations based on +/// the previous one using only additions. This is the purpose of `deltas`, to hold the increments +/// added to each multi-linear to compute the evaluation at the next point, and `evals_x` to hold +/// the current evaluation at `x` in {2, ... , d_max}. +fn sumcheck_round>( + composition_poly: &P, + mls: &[MultiLinearPoly], +) -> UnivariatePolyEvals { + let num_ml = mls.len(); + let num_vars = mls[0].num_variables(); + let num_rounds = num_vars - 1; + + let mut evals_zero = vec![E::ZERO; num_ml]; + let mut evals_one = vec![E::ZERO; num_ml]; + let mut deltas = vec![E::ZERO; num_ml]; + let mut evals_x = vec![E::ZERO; num_ml]; + + let total_evals = (0..1 << num_rounds).map(|i| { + for (j, ml) in mls.iter().enumerate() { + // Holds `f(x_{i+1}, ..., x_v, 0)` + evals_zero[j] = ml.evaluations()[2 * i]; + + // Holds `f(x_{i+1}, ..., x_v, 1)` + evals_one[j] = ml.evaluations()[2 * i + 1]; + } + + let mut total_evals = vec![E::ZERO; composition_poly.max_degree() as usize]; + total_evals[0] = composition_poly.evaluate(&evals_one); + + evals_zero + .iter() + .zip(evals_one.iter().zip(deltas.iter_mut().zip(evals_x.iter_mut()))) + .for_each(|(a0, (a1, (delta, evx)))| { + *delta = *a1 - *a0; + *evx = *a1; + }); + + total_evals.iter_mut().skip(1).for_each(|e| { + evals_x.iter_mut().zip(deltas.iter()).for_each(|(evx, delta)| { + *evx += *delta; + }); + *e = composition_poly.evaluate(&evals_x); + }); + total_evals + }); + + let evaluations = total_evals.fold( + vec![E::ZERO; composition_poly.max_degree() as usize], + |mut acc, evals| { + acc.iter_mut().zip(evals.iter()).for_each(|(a, ev)| *a += *ev); + acc + }, + ); + + UnivariatePolyEvals { + partial_evaluations: evaluations, + } +} + +pub trait FinalClaimBuilder { + type Field: FieldElement; + + fn build_claim( + &self, + openings: alloc::vec::Vec, + evaluation_point: &[Self::Field], + ) -> FinalOpeningClaim; +} diff --git a/processor/src/trace/virtual_bus/sum_check/verifier/error.rs b/processor/src/trace/virtual_bus/sum_check/verifier/error.rs new file mode 100644 index 0000000000..19163283cd --- /dev/null +++ b/processor/src/trace/virtual_bus/sum_check/verifier/error.rs @@ -0,0 +1,9 @@ +#[derive(Debug, thiserror::Error)] +pub enum Error { + #[error("the final evaluation check of sum-check failed")] + FinalEvaluationCheckFailed, + #[error("failed to generate round challenge")] + FailedToGenerateChallenge, + #[error("wrong opening point for the oracles")] + WrongOpeningPoint, +} diff --git a/processor/src/trace/virtual_bus/sum_check/verifier/mod.rs b/processor/src/trace/virtual_bus/sum_check/verifier/mod.rs new file mode 100644 index 0000000000..50c045056b --- /dev/null +++ b/processor/src/trace/virtual_bus/sum_check/verifier/mod.rs @@ -0,0 +1,152 @@ +use super::{CompositionPolynomial, FinalOpeningClaim, Proof, RoundClaim, RoundProof}; +use alloc::vec::Vec; +use core::marker::PhantomData; +use vm_core::FieldElement; +use winter_prover::crypto::{ElementHasher, RandomCoin}; + +mod error; +pub use self::error::Error; + +// SUM-CHECK VERIFIER +// ================================================================================================ + +/// A struct that contains relevant information for the execution of the multivariate sum-check +/// protocol verifier. The protocol is described in [`super::SumCheckProver`]. +/// +/// The sum-check Verifier is composed of two parts: +/// +/// 1. A multi-round interaction where it sends challenges and receives polynomials. For each +/// polynomial received it uses the sent randomness to reduce the current claim to a new one. +/// 2. A final round where the Verifier queries the multi-linear oracles it received at the outset +/// of the protocol (i.e., commitments) for their evaluations at the random point +/// `(r_0, ... , r_{\nu - 1})` where $\nu$ is the number of rounds of the sum-check protocol and +/// `r_i` is the randomness sent by the Verifier at each round. +pub struct SumCheckVerifier +where + E: FieldElement, + C: RandomCoin, + P: CompositionPolynomial, + H: ElementHasher, + V: CompositionPolyQueryBuilder, +{ + composition_poly: P, + final_query_builder: V, + _field: PhantomData, + _challenger: PhantomData, +} + +impl SumCheckVerifier +where + E: FieldElement, + C: RandomCoin, + P: CompositionPolynomial, + H: ElementHasher, + V: CompositionPolyQueryBuilder, +{ + /// Create a new [SumCheckVerifier] from a composition polynomial and final query builder. + pub fn new(composition_poly: P, final_query_builder: V) -> Self { + Self { + composition_poly, + final_query_builder, + _field: PhantomData, + _challenger: PhantomData, + } + } + + /// Verifies a sum-check proof [Proof] and returns a claim on the openings of the multi-linear + /// oracles that are part of the statement being proven. + /// + /// More precisely, the method: + /// + /// 1. Generates a `claimed_evaluation` from the round proof polynomials and the round challenge + /// randomness. + /// 2. Computes a query that is built using the [FinalQueryBuilder] from the multi-linear + /// openings and the round challenges. + /// 3. Evaluates the composition polynomial at the query and checks that it is equal + /// `claimed_evaluation`. + /// 4. Outputs a `FinalOpeningClaim` on the multi-linear oracles. + /// + /// Thus, the proof is correct if the method outputs a [FinalOpeningClaim] and this latter is + /// a valid claim on the multi-linear oracles i.e., each multi-linear oracle opens to the + /// claimed value at the specified opening point. + /// + /// # Errors + /// Returns an error if: + /// - No openings were provided. + /// - Draw the round challenge fails. + /// - The final evaluation check fails. + pub fn verify( + &self, + claim: E, + proof: Proof, + coin: &mut C, + ) -> Result, Error> { + let Proof { + openings_claim, + round_proofs, + } = proof; + + let RoundClaim { + eval_point: evaluation_point, + claim: claimed_evaluation, + } = self.verify_rounds(claim, round_proofs, coin)?; + + if openings_claim.eval_point != evaluation_point { + return Err(Error::WrongOpeningPoint); + } + let query = self.final_query_builder.build_query(&openings_claim, &evaluation_point); + + if self.composition_poly.evaluate(&query) != claimed_evaluation { + Err(Error::FinalEvaluationCheckFailed) + } else { + Ok(openings_claim) + } + } + + /// Verifies a round of the sum-check protocol. + pub fn verify_rounds( + &self, + claim: E, + round_proofs: Vec>, + coin: &mut C, + ) -> Result, Error> { + let mut round_claim = claim; + let mut evaluation_point = vec![]; + for round_proof in round_proofs { + let round_poly_coefs = round_proof.round_poly_coefs.clone(); + coin.reseed(H::hash_elements(&round_poly_coefs.coefficients)); + + let r = coin.draw().map_err(|_| Error::FailedToGenerateChallenge)?; + + round_claim = round_proof.round_poly_coefs.evaluate_using_claim(&round_claim, &r); + evaluation_point.push(r); + } + + Ok(RoundClaim { + eval_point: evaluation_point, + claim: round_claim, + }) + } +} + +/// Contains the logic for building the final query made to the virtual polynomial. +/// +/// During the last step of the sum-check protocol, the Verifier must evaluate the composed +/// multi-linear polynomials at a random point `(r_0, ... ,r_{\nu - 1})`. To do this, the Verifier +/// asks the Prover for the openings of the multi-linear oracles at `(r_0, ... ,r_{\nu - 1})` i.e., +/// `v_i = f_i(r_0, ... ,r_{\nu - 1})`. The Verifier then evaluates `g(v_0, ... , v_{\nu - 1})` and +/// compares it to the reduced claim resulting from the round proofs and challenges. +/// +/// At this point, for the Verifier to accept the proof, it needs to check that indeed +/// `v_i = f_i(r_0, ... ,r_{\nu - 1})`, this is the exact content of [`FinalOpeningClaim`], which +/// can be either answered by a direct query to the oracles (i.e., in the compiled protocol this +/// would be answered with an opening proof against the commitment) or through further interaction +/// (as in the case of the GKR protocol). +/// +/// The purpose of [`CompositionPolyQueryBuilder`] is to abstract the logic for evaluating the +/// multi-linear polynomials that the Verifier can do by herself. For example, this is the case +/// for periodic columns where given `(r_0, ... ,r_{\nu - 1})` the Verifier can evaluate +/// it at `(r_0, ... ,r_{\nu - 1})` unassisted. +pub trait CompositionPolyQueryBuilder { + fn build_query(&self, openings_claim: &FinalOpeningClaim, evaluation_point: &[E]) -> Vec; +} diff --git a/processor/src/trace/virtual_bus/tests.rs b/processor/src/trace/virtual_bus/tests.rs new file mode 100644 index 0000000000..6492c9a946 --- /dev/null +++ b/processor/src/trace/virtual_bus/tests.rs @@ -0,0 +1,81 @@ +use crate::{prove_virtual_bus, verify_virtual_bus, DefaultHost, ExecutionTrace, Process}; +use alloc::vec::Vec; +use miden_air::{ + trace::{main_trace::MainTrace, range::M_COL_IDX}, + ExecutionOptions, +}; +use vm_core::crypto::random::RpoRandomCoin; +use vm_core::{ + code_blocks::CodeBlock, CodeBlockTable, Felt, FieldElement, Kernel, Operation, StackInputs, +}; + +#[test] +fn test_vb_prover_verifier() { + let s = 6; + let o = 6; + let stack: Vec<_> = (0..(1 << s)).into_iter().collect(); + let operations: Vec<_> = (0..(1 << o)) + .flat_map(|_| { + vec![Operation::U32split, Operation::U32add, Operation::U32xor, Operation::MStoreW] + }) + .collect(); + + let trace = build_full_trace(&stack, operations, Kernel::default()); + + // this should be generated using the transcript up to when the prover sends the commitment + // to the main trace. + let log_up_randomness: Vec = vec![test_utils::rand::rand_value()]; + + let seed = [Felt::ZERO; 4]; // should be initialized with the appropriate transcript + let mut transcript = RpoRandomCoin::new(seed.into()); + let proof = prove_virtual_bus(&trace, log_up_randomness.clone(), &mut transcript).unwrap(); + + let seed = [Felt::ZERO; 4]; // should be initialized with the appropriate transcript + let mut transcript = RpoRandomCoin::new(seed.into()); + let final_opening_claim = + verify_virtual_bus(Felt::ZERO, proof, log_up_randomness, &mut transcript); + assert!(final_opening_claim.is_ok()) +} + +#[test] +fn test_vb_prover_verifier_failure() { + let s = 6; + let o = 6; + let stack: Vec<_> = (0..(1 << s)).into_iter().collect(); + let operations: Vec<_> = (0..(1 << o)) + .flat_map(|_| { + vec![Operation::U32split, Operation::U32add, Operation::U32xor, Operation::MStoreW] + }) + .collect(); + + // modifying the multiplicity + let mut trace = build_full_trace(&stack, operations, Kernel::default()); + let index = trace.get_column(M_COL_IDX).iter().position(|v| *v != Felt::ZERO).unwrap(); + trace.get_column_mut(M_COL_IDX)[index] = Felt::ONE; + + // this should be generated using the transcript up to when the prover sends the commitment + // to the main trace. + let log_up_randomness: Vec = vec![test_utils::rand::rand_value()]; + + let seed = [Felt::ZERO; 4]; // should be initialized with the appropriate transcript + let mut transcript = RpoRandomCoin::new(seed.into()); + let proof = prove_virtual_bus(&trace, log_up_randomness.clone(), &mut transcript).unwrap(); + + let seed = [Felt::ZERO; 4]; // should be initialized with the appropriate transcript + let mut transcript = RpoRandomCoin::new(seed.into()); + let final_opening_claim = + verify_virtual_bus(Felt::ZERO, proof, log_up_randomness, &mut transcript); + assert!(final_opening_claim.is_err()) +} + +fn build_full_trace(stack_inputs: &[u64], operations: Vec, kernel: Kernel) -> MainTrace { + let stack_inputs: Vec = stack_inputs.iter().map(|a| Felt::new(*a)).collect(); + let stack_inputs = StackInputs::new(stack_inputs).unwrap(); + let host = DefaultHost::default(); + let mut process = Process::new(kernel, stack_inputs, host, ExecutionOptions::default()); + let program = CodeBlock::new_span(operations); + process.execute_code_block(&program, &CodeBlockTable::default()).unwrap(); + let (trace, _, _): (MainTrace, _, _) = ExecutionTrace::test_finalize_trace(process); + + trace +} diff --git a/processor/src/trace/virtual_bus/univariate.rs b/processor/src/trace/virtual_bus/univariate.rs new file mode 100644 index 0000000000..66a3545a91 --- /dev/null +++ b/processor/src/trace/virtual_bus/univariate.rs @@ -0,0 +1,253 @@ +use alloc::vec::Vec; +use vm_core::{polynom, FieldElement}; +use winter_prover::math::batch_inversion; + +// UNIVARIATE POLYNOMIAL (EVALUATION FORM) +// ================================================================================================ + +/// The evaluations of a univariate polynomial of degree n at 0, 1, ..., n with the evaluation at 0 +/// omitted. +#[derive(Clone, Debug)] +pub struct UnivariatePolyEvals { + pub(crate) partial_evaluations: Vec, +} + +impl UnivariatePolyEvals { + /// Gives the coefficient representation of a polynomial represented in evaluation form. + /// + /// Since the evaluation at 0 is omitted, we need to use the round claim to recover + /// the evaluation at 0 using the identity p(0) + p(1) = claim. + /// Now, we have that for any polynomial p(x) = c0 + c1 * x + ... + c_{n-1} * x^{n - 1}: + /// + /// 1. p(0) = c0. + /// 2. p(x) = c0 + x * q(x) where q(x) = c1 + ... + c_{n-1} * x^{n - 2}. + /// + /// This means that we can compute the evaluations of q at 1, ..., n - 1 using the evaluations + /// of p and thus reduce by 1 the size of the interpolation problem. + /// Once the coefficient of q are recovered, the c0 coefficient is appended to these and this + /// is precisely the coefficient representation of the original polynomial q. + /// Note that the coefficient of the linear term is removed as this coefficient can be recovered + /// from the remaining coefficients, again, using the round claim using the relation + /// 2 * c0 + c1 + ... c_{n - 1} = claim. + pub fn to_poly(&self, round_claim: E) -> UnivariatePolyCoef { + // construct the vector of interpolation points 1, ..., n + let n_minus_1 = self.partial_evaluations.len(); + let mut points = vec![E::BaseField::ZERO; n_minus_1]; + + // construct their inverses. These will be needed for computing the evaluations + // of the q polynomial as well as for doing the interpolation on q + points + .iter_mut() + .enumerate() + .for_each(|(i, node)| *node = E::BaseField::from(i as u32 + 1)); + let points_inv = batch_inversion(&points); + + // compute the zeroth coefficient + let c0 = round_claim - self.partial_evaluations[0]; + + // compute the evaluations of q + let mut q_evals: Vec = vec![E::ZERO; n_minus_1]; + q_evals.iter_mut().zip(self.partial_evaluations.iter()).enumerate().for_each( + |(i, (normalized, evals))| *normalized = (*evals - c0).mul_base(points_inv[i]), + ); + + // interpolate q + let q_coefs = multiply_by_inverse_vandermonde(&q_evals, &points_inv); + + // append c0 to the coefficients of q to get the coefficients of p. The linear term + // coefficient is removed as this can be recovered from the other coefficients using + // the reduced claim. + let mut coefficients = Vec::with_capacity(self.partial_evaluations.len() + 1); + coefficients.push(c0); + coefficients.extend_from_slice(&q_coefs[1..]); + + UnivariatePolyCoef { coefficients } + } +} + +// UNIVARIATE POLYNOMIAL (COEFFICIENT FORM) +// ================================================================================================ + +/// The coefficients of a univariate polynomial of degree n with the linear term coefficient +/// omitted. +#[derive(Clone, Debug)] +pub struct UnivariatePolyCoef { + pub(crate) coefficients: Vec, +} + +impl UnivariatePolyCoef { + /// Evaluates a polynomial at a challenge point using a round claim. + /// + /// The round claim is used to recover the coefficient of the linear term using the relation + /// 2 * c0 + c1 + ... c_{n - 1} = claim. Using the complete list of coefficients, the polynomial + /// is then evaluated using Horner's method. + pub fn evaluate_using_claim(&self, claim: &E, challenge: &E) -> E { + // recover the coefficient of the linear term + let c1 = *claim + - self.coefficients.iter().fold(E::ZERO, |acc, term| acc + *term) + - self.coefficients[0]; + + // construct the full coefficient list + let mut complete_coefficients = vec![self.coefficients[0], c1]; + complete_coefficients.extend_from_slice(&self.coefficients[1..]); + + // evaluate + polynom::eval(&complete_coefficients, *challenge) + } +} + +// HELPER FUNCTIONS +// ================================================================================================ + +/// Given a (row) vector `v`, computes the vector-matrix product `v * V^{-1}` where `V` is +/// the Vandermonde matrix over the points `1, ..., n` where `n` is the length of `v`. +/// The resulting vector will then be the coefficients of the minimal interpolating polynomial +/// through the points `(i+1, v[i])` for `i` in `0, ..., n - 1` +/// +/// The naive way would be to invert the matrix `V` and then compute the vector-matrix product +/// this will cost `O(n^3)` operations and `O(n^2)` memory. We can also try Gaussian elimination +/// but this is also worst case `O(n^3)` operations and `O(n^2)` memory. +/// In the following implementation, we use the fact that the points over which we are interpolating +/// is a set of equidistant points and thus both the Vandermonde matrix and its inverse can be +/// described by sparse linear recurrence equations. +/// More specifically, we use the representation given in [1], where `V^{-1}` is represented as +/// `U * M` where: +/// +/// 1. `M` is a lower triangular matrix where its entries are given by M(i, j) = M(i - 1, j) - M(i - +/// 1, j - 1) / (i - 1) +/// with boundary conditions M(i, 1) = 1 and M(i, j) = 0 when j > i. +/// +/// 2. `U` is an upper triangular (involutory) matrix where its entries are given by U(i, j) = U(i, +/// j - 1) - U(i - 1, j - 1) +/// with boundary condition U(1, j) = 1 and U(i, j) = 0 when i > j. +/// +/// Note that the matrix indexing in the formulas above matches the one in the reference and starts +/// from 1. +/// +/// The above implies that we can do the vector-matrix multiplication in `O(n^2)` and using only +/// `O(n)` space. +/// +/// [1]: https://link.springer.com/article/10.1007/s002110050360 +fn multiply_by_inverse_vandermonde( + vector: &[E], + nodes_inv: &[E::BaseField], +) -> Vec { + let res = multiply_by_u(vector); + multiply_by_m(&res, nodes_inv) +} + +/// Multiplies a (row) vector `v` by an upper triangular matrix `U` to compute `v * U`. +/// +/// `U` is an upper triangular (involutory) matrix with its entries given by +/// U(i, j) = U(i, j - 1) - U(i - 1, j - 1) +/// with boundary condition U(1, j) = 1 and U(i, j) = 0 when i > j. +fn multiply_by_u(vector: &[E]) -> Vec { + let n = vector.len(); + let mut previous_u_col = vec![E::BaseField::ZERO; n]; + previous_u_col[0] = E::BaseField::ONE; + let mut current_u_col = vec![E::BaseField::ZERO; n]; + current_u_col[0] = E::BaseField::ONE; + + let mut result: Vec = vec![E::ZERO; n]; + for (i, res) in result.iter_mut().enumerate() { + *res = vector[0]; + + for (j, v) in vector.iter().enumerate().take(i + 1).skip(1) { + let u_entry: E::BaseField = + compute_u_entry::(j, &mut previous_u_col, &mut current_u_col); + *res += v.mul_base(u_entry); + } + previous_u_col.clone_from(¤t_u_col); + } + + result +} + +/// Multiplies a (row) vector `v` by a lower triangular matrix `M` to compute `v * M`. +/// +/// `M` is a lower triangular matrix with its entries given by +/// M(i, j) = M(i - 1, j) - M(i - 1, j - 1) / (i - 1) +/// with boundary conditions M(i, 1) = 1 and M(i, j) = 0 when j > i. +fn multiply_by_m(vector: &[E], nodes_inv: &[E::BaseField]) -> Vec { + let n = vector.len(); + let mut previous_m_col = vec![E::BaseField::ONE; n]; + let mut current_m_col = vec![E::BaseField::ZERO; n]; + current_m_col[0] = E::BaseField::ONE; + + let mut result: Vec = vec![E::ZERO; n]; + result[0] = vector.iter().fold(E::ZERO, |acc, term| acc + *term); + for (i, res) in result.iter_mut().enumerate().skip(1) { + current_m_col = vec![E::BaseField::ZERO; n]; + + for (j, v) in vector.iter().enumerate().skip(i) { + let m_entry: E::BaseField = + compute_m_entry::(j, &mut previous_m_col, &mut current_m_col, nodes_inv[j - 1]); + *res += v.mul_base(m_entry); + } + previous_m_col.clone_from(¤t_m_col); + } + + result +} + +/// Returns the j-th entry of the i-th column of matrix `U` given the values of the (i - 1)-th +/// column. The i-th column is also updated with the just computed `U(i, j)` entry. +/// +/// `U` is an upper triangular (involutory) matrix with its entries given by +/// U(i, j) = U(i, j - 1) - U(i - 1, j - 1) +/// with boundary condition U(1, j) = 1 and U(i, j) = 0 when i > j. +fn compute_u_entry( + j: usize, + col_prev: &mut [E::BaseField], + col_cur: &mut [E::BaseField], +) -> E::BaseField { + let value = col_prev[j] - col_prev[j - 1]; + col_cur[j] = value; + value +} + +/// Returns the j-th entry of the i-th column of matrix `M` given the values of the (i - 1)-th +/// and the i-th columns. The i-th column is also updated with the just computed `M(i, j)` entry. +/// +/// `M` is a lower triangular matrix with its entries given by +/// M(i, j) = M(i - 1, j) - M(i - 1, j - 1) / (i - 1) +/// with boundary conditions M(i, 1) = 1 and M(i, j) = 0 when j > i. +fn compute_m_entry( + j: usize, + col_previous: &mut [E::BaseField], + col_current: &mut [E::BaseField], + node_inv: E::BaseField, +) -> E::BaseField { + let value = col_current[j - 1] - node_inv * col_previous[j - 1]; + col_current[j] = value; + value +} + +// TESTS +// ================================================================================================ + +#[test] +fn test_poly_partial() { + use vm_core::Felt; + + use test_utils::rand; + let degree = 1000; + let mut points: Vec = vec![Felt::ZERO; degree]; + points.iter_mut().enumerate().for_each(|(i, node)| *node = Felt::from(i as u32)); + + let p: Vec = rand::rand_vector(degree); + let evals = vm_core::polynom::eval_many(&p, &points); + + let mut partial_evals = evals.clone(); + partial_evals.remove(0); + + let partial_poly = UnivariatePolyEvals { + partial_evaluations: partial_evals, + }; + let claim = evals[0] + evals[1]; + let poly_coeff = partial_poly.to_poly(claim); + + let r = rand::rand_vector(1); + + assert_eq!(vm_core::polynom::eval(&p, r[0]), poly_coeff.evaluate_using_claim(&claim, &r[0])) +}