diff --git a/Cargo.toml b/Cargo.toml index eed7f1be3..81a003780 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [workspace] -members = ["math", "crypto", "gpu", "benches", "provers/plonk", "provers/stark", "provers/groth16", "provers/groth16/arkworks-adapter", "provers/groth16/circom-adapter", "examples/merkle-tree-cli", "examples/prove-miden", "provers/winterfell_adapter", "examples/shamir_secret_sharing","examples/pinocchio", "examples/prove-verify-circom", "examples/baby-snark"] +members = ["math", "crypto", "gpu", "benches", "provers/plonk", "provers/stark", "provers/groth16", "provers/groth16/arkworks-adapter", "provers/groth16/circom-adapter", "examples/merkle-tree-cli", "examples/prove-miden", "provers/winterfell_adapter", "examples/shamir_secret_sharing","examples/pinocchio", "examples/prove-verify-circom", "examples/baby-snark", "provers/circle_stark"] exclude = ["ensure-no_std"] resolver = "2" diff --git a/math/src/circle/cfft.rs b/math/src/circle/cfft.rs new file mode 100644 index 000000000..f060e02d6 --- /dev/null +++ b/math/src/circle/cfft.rs @@ -0,0 +1,220 @@ +extern crate alloc; +use crate::field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}; +use alloc::vec::Vec; + +#[cfg(feature = "alloc")] +/// fft in place algorithm used to evaluate a polynomial of degree 2^n - 1 in 2^n points. +/// Input must be of size 2^n for some n. +pub fn cfft( + input: &mut [FieldElement], + twiddles: Vec>>, +) { + // If the input size is 2^n, then log_2_size is n. + let log_2_size = input.len().trailing_zeros(); + + // The cfft has n layers. + (0..log_2_size).for_each(|i| { + // In each layer i we split the current input in chunks of size 2^{i+1}. + let chunk_size = 1 << (i + 1); + let half_chunk_size = 1 << i; + input.chunks_mut(chunk_size).for_each(|chunk| { + // We split each chunk in half, calling the first half hi_part and the second hal low_part. + let (hi_part, low_part) = chunk.split_at_mut(half_chunk_size); + + // We apply the corresponding butterfly for every element j of the high and low part. + hi_part + .iter_mut() + .zip(low_part) + .enumerate() + .for_each(|(j, (hi, low))| { + let temp = *low * twiddles[i as usize][j]; + *low = *hi - temp; + *hi += temp + }); + }); + }); +} + +#[cfg(feature = "alloc")] +/// The inverse fft algorithm used to interpolate 2^n points. +/// Input must be of size 2^n for some n. +pub fn icfft( + input: &mut [FieldElement], + twiddles: Vec>>, +) { + // If the input size is 2^n, then log_2_size is n. + let log_2_size = input.len().trailing_zeros(); + + // The icfft has n layers. + (0..log_2_size).for_each(|i| { + // In each layer i we split the current input in chunks of size 2^{n - i}. + let chunk_size = 1 << (log_2_size - i); + let half_chunk_size = chunk_size >> 1; + input.chunks_mut(chunk_size).for_each(|chunk| { + // We split each chunk in half, calling the first half hi_part and the second hal low_part. + let (hi_part, low_part) = chunk.split_at_mut(half_chunk_size); + + // We apply the corresponding butterfly for every element j of the high and low part. + hi_part + .iter_mut() + .zip(low_part) + .enumerate() + .for_each(|(j, (hi, low))| { + let temp = *hi + *low; + *low = (*hi - *low) * twiddles[i as usize][j]; + *hi = temp; + }); + }); + }); +} + +/// This function permutes a slice of field elements to order the result of the cfft in the natural way. +/// We call the natural order to [P(x0, y0), P(x1, y1), P(x2, y2), ...], +/// where (x0, y0) is the first point of the corresponding coset. +/// The cfft doesn't return the evaluations in the natural order. +/// For example, if we apply the cfft to 8 coefficients of a polynomial of degree 7 we'll get the evaluations in this order: +/// [P(x0, y0), P(x2, y2), P(x4, y4), P(x6, y6), P(x7, y7), P(x5, y5), P(x3, y3), P(x1, y1)], +/// where the even indices are found first in ascending order and then the odd indices in descending order. +/// This function permutes the slice [0, 2, 4, 6, 7, 5, 3, 1] into [0, 1, 2, 3, 4, 5, 6, 7]. +/// TODO: This can be optimized by performing in-place value swapping (WIP). +pub fn order_cfft_result_naive( + input: &[FieldElement], +) -> Vec> { + let mut result = Vec::new(); + let length = input.len(); + for i in 0..length / 2 { + result.push(input[i]); // We push the left index. + result.push(input[length - i - 1]); // We push the right index. + } + result +} + +/// This function permutes a slice of field elements to order the input of the icfft in a specific way. +/// For example, if we want to interpolate 8 points we should input them in the icfft in this order: +/// [(x0, y0), (x2, y2), (x4, y4), (x6, y6), (x7, y7), (x5, y5), (x3, y3), (x1, y1)], +/// where the even indices are found first in ascending order and then the odd indices in descending order. +/// This function permutes the slice [0, 1, 2, 3, 4, 5, 6, 7] into [0, 2, 4, 6, 7, 5, 3, 1]. +/// TODO: This can be optimized by performing in-place value swapping (WIP). +pub fn order_icfft_input_naive( + input: &mut [FieldElement], +) -> Vec> { + let mut result = Vec::new(); + + // We push the even indices. + (0..input.len()).step_by(2).for_each(|i| { + result.push(input[i]); + }); + + // We push the odd indices. + (1..input.len()).step_by(2).rev().for_each(|i| { + result.push(input[i]); + }); + result +} + +#[cfg(test)] +mod tests { + use super::*; + type FE = FieldElement; + + #[test] + fn ordering_cfft_result_works_for_4_points() { + let expected_slice = [FE::from(0), FE::from(1), FE::from(2), FE::from(3)]; + + let slice = [FE::from(0), FE::from(2), FE::from(3), FE::from(1)]; + + let res = order_cfft_result_naive(&slice); + + assert_eq!(res, expected_slice) + } + + #[test] + fn ordering_cfft_result_works_for_16_points() { + let expected_slice = [ + FE::from(0), + FE::from(1), + FE::from(2), + FE::from(3), + FE::from(4), + FE::from(5), + FE::from(6), + FE::from(7), + FE::from(8), + FE::from(9), + FE::from(10), + FE::from(11), + FE::from(12), + FE::from(13), + FE::from(14), + FE::from(15), + ]; + + let slice = [ + FE::from(0), + FE::from(2), + FE::from(4), + FE::from(6), + FE::from(8), + FE::from(10), + FE::from(12), + FE::from(14), + FE::from(15), + FE::from(13), + FE::from(11), + FE::from(9), + FE::from(7), + FE::from(5), + FE::from(3), + FE::from(1), + ]; + + let res = order_cfft_result_naive(&slice); + + assert_eq!(res, expected_slice) + } + + #[test] + fn from_natural_to_icfft_input_order_works() { + let mut slice = [ + FE::from(0), + FE::from(1), + FE::from(2), + FE::from(3), + FE::from(4), + FE::from(5), + FE::from(6), + FE::from(7), + FE::from(8), + FE::from(9), + FE::from(10), + FE::from(11), + FE::from(12), + FE::from(13), + FE::from(14), + FE::from(15), + ]; + + let expected_slice = [ + FE::from(0), + FE::from(2), + FE::from(4), + FE::from(6), + FE::from(8), + FE::from(10), + FE::from(12), + FE::from(14), + FE::from(15), + FE::from(13), + FE::from(11), + FE::from(9), + FE::from(7), + FE::from(5), + FE::from(3), + FE::from(1), + ]; + + let res = order_icfft_input_naive(&mut slice); + + assert_eq!(res, expected_slice) + } +} diff --git a/math/src/circle/cosets.rs b/math/src/circle/cosets.rs new file mode 100644 index 000000000..b9fcdafcc --- /dev/null +++ b/math/src/circle/cosets.rs @@ -0,0 +1,95 @@ +extern crate alloc; +use crate::circle::point::CirclePoint; +use crate::field::fields::mersenne31::field::Mersenne31Field; +use alloc::vec::Vec; + +/// Given g_n, a generator of the subgroup of size n of the circle, i.e. , +/// and given a shift, that is a another point of the circle, +/// we define the coset shift + which is the set of all the points in +/// plus the shift. +/// For example, if = {p1, p2, p3, p4}, then g_8 + = {g_8 + p1, g_8 + p2, g_8 + p3, g_8 + p4}. + +#[derive(Debug, Clone)] +pub struct Coset { + // Coset: shift + where n = 2^{log_2_size}. + // Example: g_16 + , n = 8, log_2_size = 3, shift = g_16. + pub log_2_size: u32, //TODO: Change log_2_size to u8 because log_2_size < 31. + pub shift: CirclePoint, +} + +impl Coset { + pub fn new(log_2_size: u32, shift: CirclePoint) -> Self { + Coset { log_2_size, shift } + } + + /// Returns the coset g_2n + + pub fn new_standard(log_2_size: u32) -> Self { + // shift is a generator of the subgroup of order 2n = 2^{log_2_size + 1}. + let shift = CirclePoint::get_generator_of_subgroup(log_2_size + 1); + Coset { log_2_size, shift } + } + + /// Returns g_n, the generator of the subgroup of order n = 2^log_2_size. + pub fn get_generator(&self) -> CirclePoint { + CirclePoint::GENERATOR.repeated_double(31 - self.log_2_size) + } + + /// Given a standard coset g_2n + , returns the subcoset with half size g_2n + + pub fn half_coset(coset: Self) -> Self { + Coset { + log_2_size: coset.log_2_size - 1, + shift: coset.shift, + } + } + + /// Given a coset shift + G returns the coset -shift + G. + /// Note that (g_2n + ) U (-g_2n + ) = g_2n + . + pub fn conjugate(coset: Self) -> Self { + Coset { + log_2_size: coset.log_2_size, + shift: coset.shift.conjugate(), + } + } + + /// Returns the vector of shift + g for every g in . + /// where g = i * g_n for i = 0, ..., n-1. + #[cfg(feature = "alloc")] + pub fn get_coset_points(coset: &Self) -> Vec> { + // g_n the generator of the subgroup of order n. + let generator_n = CirclePoint::get_generator_of_subgroup(coset.log_2_size); + let size: usize = 1 << coset.log_2_size; + core::iter::successors(Some(coset.shift.clone()), move |prev| { + Some(prev + &generator_n) + }) + .take(size.into()) + .collect() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn coset_points_vector_has_right_size() { + let coset = Coset::new_standard(3); + let points = Coset::get_coset_points(&coset); + assert_eq!(1 << coset.log_2_size, points.len()) + } + + #[test] + fn antipode_of_coset_point_is_in_coset() { + let coset = Coset::new_standard(3); + let points = Coset::get_coset_points(&coset); + let point = points[2].clone(); + let anitpode_point = points[6].clone(); + assert_eq!(anitpode_point, point.antipode()) + } + + #[test] + fn coset_generator_has_right_order() { + let coset = Coset::new(2, CirclePoint::GENERATOR * 3); + let generator_n = coset.get_generator(); + assert_eq!(generator_n.repeated_double(2), CirclePoint::zero()); + } +} diff --git a/math/src/circle/errors.rs b/math/src/circle/errors.rs new file mode 100644 index 000000000..51dcb720b --- /dev/null +++ b/math/src/circle/errors.rs @@ -0,0 +1,4 @@ +#[derive(Debug)] +pub enum CircleError { + PointDoesntSatisfyCircleEquation, +} diff --git a/math/src/circle/mod.rs b/math/src/circle/mod.rs new file mode 100644 index 000000000..ac576194f --- /dev/null +++ b/math/src/circle/mod.rs @@ -0,0 +1,6 @@ +pub mod cfft; +pub mod cosets; +pub mod errors; +pub mod point; +pub mod polynomial; +pub mod twiddles; diff --git a/math/src/circle/point.rs b/math/src/circle/point.rs new file mode 100644 index 000000000..e0e7aa210 --- /dev/null +++ b/math/src/circle/point.rs @@ -0,0 +1,302 @@ +use super::errors::CircleError; +use crate::field::traits::IsField; +use crate::field::{ + element::FieldElement, + fields::mersenne31::{extensions::Degree4ExtensionField, field::Mersenne31Field}, +}; +use core::ops::{Add, AddAssign, Mul, MulAssign}; + +/// Given a Field F, we implement here the Group which consists of all the points (x, y) such as +/// x in F, y in F and x^2 + y^2 = 1, i.e. the Circle. The operation of the group will have +/// additive notation and is as follows: +/// (a, b) + (c, d) = (a * c - b * d, a * d + b * c) +#[derive(Debug, Clone)] +pub struct CirclePoint { + pub x: FieldElement, + pub y: FieldElement, +} + +impl> CirclePoint { + pub fn new(x: FieldElement, y: FieldElement) -> Result { + if x.square() + y.square() == FieldElement::one() { + Ok(Self { x, y }) + } else { + Err(CircleError::PointDoesntSatisfyCircleEquation) + } + } + + /// Neutral element of the Circle group (with additive notation). + pub fn zero() -> Self { + Self::new(FieldElement::one(), FieldElement::zero()).unwrap() + } + + /// Computes 2(x, y) = (2x^2 - 1, 2xy). + pub fn double(&self) -> Self { + Self::new( + self.x.square().double() - FieldElement::one(), + self.x.double() * self.y.clone(), + ) + .unwrap() + } + + /// Computes 2^n * (x, y). + pub fn repeated_double(self, n: u32) -> Self { + let mut res = self; + for _ in 0..n { + res = res.double(); + } + res + } + + /// Computes the inverse of the point. + /// We are using -(x, y) = (x, -y), i.e. the inverse of the group opertion is conjugation + /// because the norm of every point in the circle is one. + pub fn conjugate(self) -> Self { + Self { + x: self.x, + y: -self.y, + } + } + + pub fn antipode(self) -> Self { + Self { + x: -self.x, + y: -self.y, + } + } + + pub const GENERATOR: Self = Self { + x: F::CIRCLE_GENERATOR_X, + y: F::CIRCLE_GENERATOR_Y, + }; + + /// Returns the generator of the subgroup of order n = 2^log_2_size. + /// We are using that 2^k * g is a generator of the subgroup of order 2^{31 - k}. + pub fn get_generator_of_subgroup(log_2_size: u32) -> Self { + Self::GENERATOR.repeated_double(31 - log_2_size) + } + + pub const ORDER: u128 = F::ORDER; +} + +/// Parameters of the base field that we'll need to define its Circle. +pub trait HasCircleParams { + type FE; + + /// Coordinate x of the generator of the circle group. + const CIRCLE_GENERATOR_X: FieldElement; + + /// Coordinate y of the generator of the circle group. + const CIRCLE_GENERATOR_Y: FieldElement; + + const ORDER: u128; +} + +impl HasCircleParams for Mersenne31Field { + type FE = FieldElement; + + const CIRCLE_GENERATOR_X: Self::FE = Self::FE::const_from_raw(2); + + const CIRCLE_GENERATOR_Y: Self::FE = Self::FE::const_from_raw(1268011823); + + /// ORDER = 2^31 + const ORDER: u128 = 2147483648; +} + +impl HasCircleParams for Degree4ExtensionField { + type FE = FieldElement; + + // These parameters were taken from stwo's implementation: + // https://github.com/starkware-libs/stwo/blob/9cfd48af4e8ac5dd67643a92927c894066fa989c/crates/prover/src/core/circle.rs + const CIRCLE_GENERATOR_X: Self::FE = + Degree4ExtensionField::const_from_coefficients(1, 0, 478637715, 513582971); + + const CIRCLE_GENERATOR_Y: Self::FE = + Degree4ExtensionField::const_from_coefficients(992285211, 649143431, 740191619, 1186584352); + + /// ORDER = (2^31 - 1)^4 - 1 + const ORDER: u128 = 21267647892944572736998860269687930880; +} + +/// Equality between two cricle points. +impl> PartialEq for CirclePoint { + fn eq(&self, other: &Self) -> bool { + self.x == other.x && self.y == other.y + } +} + +/// Addition (i.e. group operation) between two points: +/// (a, b) + (c, d) = (a * c - b * d, a * d + b * c) +impl> Add for &CirclePoint { + type Output = CirclePoint; + fn add(self, other: Self) -> Self::Output { + let x = &self.x * &other.x - &self.y * &other.y; + let y = &self.x * &other.y + &self.y * &other.x; + CirclePoint { x, y } + } +} +impl> Add for CirclePoint { + type Output = CirclePoint; + fn add(self, rhs: CirclePoint) -> Self::Output { + &self + &rhs + } +} +impl> Add> for &CirclePoint { + type Output = CirclePoint; + fn add(self, rhs: CirclePoint) -> Self::Output { + self + &rhs + } +} +impl> Add<&CirclePoint> for CirclePoint { + type Output = CirclePoint; + fn add(self, rhs: &CirclePoint) -> Self::Output { + &self + rhs + } +} +impl> AddAssign<&CirclePoint> for CirclePoint { + fn add_assign(&mut self, rhs: &CirclePoint) { + *self = &*self + rhs; + } +} +impl> AddAssign> for CirclePoint { + fn add_assign(&mut self, rhs: CirclePoint) { + *self += &rhs; + } +} +/// Multiplication between a point and a scalar (i.e. group operation repeatedly): +/// (x, y) * n = (x ,y) + ... + (x, y) n-times. +impl> Mul for &CirclePoint { + type Output = CirclePoint; + fn mul(self, scalar: u128) -> Self::Output { + let mut scalar = scalar; + let mut res = CirclePoint::::zero(); + let mut cur = self.clone(); + loop { + if scalar == 0 { + return res; + } + if scalar & 1 == 1 { + res += &cur; + } + cur = cur.double(); + scalar >>= 1; + } + } +} +impl> Mul for CirclePoint { + type Output = CirclePoint; + fn mul(self, scalar: u128) -> Self::Output { + &self * scalar + } +} +impl> MulAssign for CirclePoint { + fn mul_assign(&mut self, scalar: u128) { + let mut scalar = scalar; + let mut res = CirclePoint::::zero(); + loop { + if scalar == 0 { + *self = res.clone(); + } + if scalar & 1 == 1 { + res += &*self; + } + *self = self.double(); + scalar >>= 1; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + type F = Mersenne31Field; + type FE = FieldElement; + type G = CirclePoint; + + type Fp4 = Degree4ExtensionField; + type Fp4E = FieldElement; + type G4 = CirclePoint; + + #[test] + fn create_new_valid_g_point() { + let valid_point = G::new(FE::one(), FE::zero()).unwrap(); + let expected = G { + x: FE::one(), + y: FE::zero(), + }; + assert_eq!(valid_point, expected) + } + + #[test] + fn create_new_valid_g4_point() { + let valid_point = G4::new(Fp4E::one(), Fp4E::zero()).unwrap(); + let expected = G4 { + x: Fp4E::one(), + y: Fp4E::zero(), + }; + assert_eq!(valid_point, expected) + } + + #[test] + fn create_new_invalid_circle_point() { + let invalid_point = G::new(FE::one(), FE::one()); + assert!(invalid_point.is_err()) + } + + #[test] + fn create_new_invalid_g4_circle_point() { + let invalid_point = G4::new(Fp4E::one(), Fp4E::one()); + assert!(invalid_point.is_err()) + } + + #[test] + fn zero_plus_zero_is_zero() { + let a = G::zero(); + let b = G::zero(); + assert_eq!(&a + &b, G::zero()) + } + + #[test] + fn generator_plus_zero_is_generator() { + let g = G::GENERATOR; + let zero = G::zero(); + assert_eq!(&g + &zero, g) + } + + #[test] + fn double_equals_mul_two() { + let g = G::GENERATOR; + assert_eq!(g.clone().double(), g * 2) + } + + #[test] + fn mul_eight_equals_double_three_times() { + let g = G::GENERATOR; + assert_eq!(g.clone().repeated_double(3), g * 8) + } + + #[test] + fn generator_g1_has_order_two_pow_31() { + let g = G::GENERATOR; + let n = 31; + assert_eq!(g.repeated_double(n), G::zero()) + } + + #[test] + fn generator_g4_has_the_order_of_the_group() { + let g = G4::GENERATOR; + assert_eq!(g * G4::ORDER, G4::zero()) + } + + #[test] + fn conjugation_is_inverse_operation() { + let g = G::GENERATOR; + assert_eq!(&g.clone() + &g.conjugate(), G::zero()) + } + + #[test] + fn subgroup_generator_has_correct_order() { + let generator_n = G::get_generator_of_subgroup(7); + assert_eq!(generator_n.repeated_double(7), G::zero()); + } +} diff --git a/math/src/circle/polynomial.rs b/math/src/circle/polynomial.rs new file mode 100644 index 000000000..4cab7c0c3 --- /dev/null +++ b/math/src/circle/polynomial.rs @@ -0,0 +1,353 @@ +extern crate alloc; +use super::point::CirclePoint; +#[cfg(feature = "alloc")] +use super::{ + cfft::{cfft, icfft, order_cfft_result_naive, order_icfft_input_naive}, + cosets::Coset, + twiddles::{get_twiddles, TwiddlesConfig}, +}; +use crate::{ + fft::cpu::bit_reversing::in_place_bit_reverse_permute, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; +use alloc::vec::Vec; + +/// Given the 2^n coefficients of a two-variables polynomial of degree 2^n - 1 in the basis {1, y, x, xy, 2xˆ2 -1, 2xˆ2y-y, 2xˆ3-x, 2xˆ3y-xy,...} +/// returns the evaluation of the polynomial on the points of the standard coset of size 2^n. +/// Note that coeff has to be a vector with length a power of two 2^n. +#[cfg(feature = "alloc")] +pub fn evaluate_cfft( + coeff: Vec>, +) -> Vec> { + let mut coeff = coeff; + + // We get the twiddles for the Evaluation. + let domain_log_2_size: u32 = coeff.len().trailing_zeros(); + let coset = Coset::new_standard(domain_log_2_size); + let config = TwiddlesConfig::Evaluation; + let twiddles = get_twiddles(coset, config); + + // For our algorithm to work, we must give as input the coefficients in bit reverse order. + in_place_bit_reverse_permute::>(&mut coeff); + cfft(&mut coeff, twiddles); + + // The cfft returns the evaluations in a certain order, so we permute them to get the natural order. + order_cfft_result_naive(&coeff) +} + +/// Interpolates the 2^n evaluations of a two-variables polynomial of degree 2^n - 1 on the points of the standard coset of size 2^n. +/// As a result we obtain the coefficients of the polynomial in the basis: {1, y, x, xy, 2xˆ2 -1, 2xˆ2y-y, 2xˆ3-x, 2xˆ3y-xy,...} +/// Note that eval has to be a vector of length a power of two 2^n. +/// If the vector of evaluations is empty, it returns an empty vector. +#[cfg(feature = "alloc")] +pub fn interpolate_cfft( + eval: Vec>, +) -> Vec> { + let mut eval = eval; + + if eval.is_empty() { + let poly: Vec> = Vec::new(); + return poly; + } + + // We get the twiddles for the interpolation. + let domain_log_2_size: u32 = eval.len().trailing_zeros(); + let coset = Coset::new_standard(domain_log_2_size); + let config = TwiddlesConfig::Interpolation; + let twiddles = get_twiddles(coset, config); + + // For our algorithm to work, we must give as input the evaluations ordered in a certain way. + let mut eval_ordered = order_icfft_input_naive(&mut eval); + icfft(&mut eval_ordered, twiddles); + + // The icfft returns the polynomial coefficients in bit reverse order. So we premute it to get the natural order. + in_place_bit_reverse_permute::>(&mut eval_ordered); + + // The icfft returns all the coefficients multiplied by 2^n, the length of the evaluations. + // So we multiply every element that outputs the icfft by the inverse of 2^n to get the actual coefficients. + // Note that this `unwrap` will never panic because eval.len() != 0. + let factor = (FieldElement::::from(eval.len() as u64)) + .inv() + .unwrap(); + eval_ordered.iter().map(|coef| coef * factor).collect() +} + +/// Note: This implementation uses a straightforward approach and is intended for testing purposes only. +pub fn evaluate_point( + coef: &Vec>, + point: &CirclePoint, +) -> FieldElement { + let order = coef.len(); + assert!( + order.is_power_of_two(), + "Coefficient length must be a power of 2" + ); + + let v_len = order.trailing_zeros() as usize; + + let mut a = point.x; + let mut v = Vec::with_capacity(v_len); + v.push(FieldElement::one()); + v.push(point.x); + for _ in 2..v_len { + a = a.square().double() - FieldElement::one(); + v.push(a); + } + + let mut result = FieldElement::zero(); + + for i in 0..order { + let mut term = coef[i]; + + if i % 2 == 1 { + term = term * point.y; + } + + let mut idx = i / 2; + let mut pos = 0; + while idx > 0 { + if idx % 2 == 1 { + term = term * v[pos + 1]; + } + idx /= 2; + pos += 1; + } + + result = result + term; + } + + result +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::circle::cosets::Coset; + type FE = FieldElement; + use alloc::vec; + + /// Naive evaluation of a polynomial of degree 3. + fn evaluate_poly_4(coef: &[FE; 4], x: FE, y: FE) -> FE { + coef[0] + coef[1] * y + coef[2] * x + coef[3] * x * y + } + + /// Naive evaluation of a polynomial of degree 7. + fn evaluate_poly_8(coef: &[FE; 8], x: FE, y: FE) -> FE { + coef[0] + + coef[1] * y + + coef[2] * x + + coef[3] * x * y + + coef[4] * (x.square().double() - FE::one()) + + coef[5] * (x.square().double() - FE::one()) * y + + coef[6] * ((x.square() * x).double() - x) + + coef[7] * ((x.square() * x).double() - x) * y + } + + /// Naive evaluation of a polynomial of degree 15. + fn evaluate_poly_16(coef: &[FE; 16], x: FE, y: FE) -> FE { + let mut a = x; + let mut v = Vec::new(); + v.push(FE::one()); + v.push(x); + for _ in 2..4 { + a = a.square().double() - FE::one(); + v.push(a); + } + + coef[0] * v[0] + + coef[1] * y * v[0] + + coef[2] * v[1] + + coef[3] * y * v[1] + + coef[4] * v[2] + + coef[5] * y * v[2] + + coef[6] * v[1] * v[2] + + coef[7] * y * v[1] * v[2] + + coef[8] * v[3] + + coef[9] * y * v[3] + + coef[10] * v[1] * v[3] + + coef[11] * y * v[1] * v[3] + + coef[12] * v[2] * v[3] + + coef[13] * y * v[2] * v[3] + + coef[14] * v[1] * v[2] * v[3] + + coef[15] * y * v[1] * v[2] * v[3] + } + + #[test] + /// cfft evaluation equals naive evaluation. + fn cfft_evaluation_4_points() { + // We define the coefficients of a polynomial of degree 3. + let input = [FE::from(1), FE::from(2), FE::from(3), FE::from(4)]; + + // We create the coset points and evaluate the polynomial with the naive function. + let coset = Coset::new_standard(2); + let points = Coset::get_coset_points(&coset); + let mut expected_result: Vec = Vec::new(); + for point in points { + let point_eval = evaluate_poly_4(&input, point.x, point.y); + expected_result.push(point_eval); + } + + let input_vec = input.to_vec(); + // We evaluate the polynomial using now the cfft. + let result = evaluate_cfft(input_vec); + let slice_result: &[FE] = &result; + + assert_eq!(slice_result, expected_result); + } + + #[test] + /// cfft evaluation equals naive evaluation. + fn cfft_evaluation_8_points() { + // We define the coefficients of a polynomial of degree 7. + let input = [ + FE::from(1), + FE::from(2), + FE::from(3), + FE::from(4), + FE::from(5), + FE::from(6), + FE::from(7), + FE::from(8), + ]; + + // We create the coset points and evaluate them without the fft. + let coset = Coset::new_standard(3); + let points = Coset::get_coset_points(&coset); + let mut expected_result: Vec = Vec::new(); + for point in points { + let point_eval = evaluate_poly_8(&input, point.x, point.y); + expected_result.push(point_eval); + } + + // We evaluate the polynomial using now the cfft. + let result = evaluate_cfft(input.to_vec()); + let slice_result: &[FE] = &result; + + assert_eq!(slice_result, expected_result); + } + + #[test] + /// cfft evaluation equals naive evaluation. + fn cfft_evaluation_16_points() { + // We define the coefficients of a polynomial of degree 15. + let input = [ + FE::from(1), + FE::from(2), + FE::from(3), + FE::from(4), + FE::from(5), + FE::from(6), + FE::from(7), + FE::from(8), + FE::from(9), + FE::from(10), + FE::from(11), + FE::from(12), + FE::from(13), + FE::from(14), + FE::from(15), + FE::from(16), + ]; + + // We create the coset points and evaluate them without the fft. + let coset = Coset::new_standard(4); + let points = Coset::get_coset_points(&coset); + let mut expected_result: Vec = Vec::new(); + for point in points { + let point_eval = evaluate_poly_16(&input, point.x, point.y); + expected_result.push(point_eval); + } + + // We evaluate the polynomial using now the cfft. + let result = evaluate_cfft(input.to_vec()); + let slice_result: &[FE] = &result; + + assert_eq!(slice_result, expected_result); + } + + #[test] + fn evaluate_and_interpolate_8_points_is_identity() { + // We define the 8 coefficients of a polynomial of degree 7. + let coeff = vec![ + FE::from(1), + FE::from(2), + FE::from(3), + FE::from(4), + FE::from(5), + FE::from(6), + FE::from(7), + FE::from(8), + ]; + let evals = evaluate_cfft(coeff.clone()); + let new_coeff = interpolate_cfft(evals); + + assert_eq!(coeff, new_coeff); + } + + #[test] + fn evaluate_and_interpolate_8_other_points() { + let coeff = vec![ + FE::from(2147483650), + FE::from(147483647), + FE::from(2147483700), + FE::from(2147483647), + FE::from(3147483647), + FE::from(4147483647), + FE::from(2147483640), + FE::from(5147483647), + ]; + let evals = evaluate_cfft(coeff.clone()); + let new_coeff = interpolate_cfft(evals); + + assert_eq!(coeff, new_coeff); + } + + #[test] + fn evaluate_and_interpolate_32_points() { + // We define 32 coefficients of a polynomial of degree 31. + let coeff = vec![ + FE::from(1), + FE::from(2), + FE::from(3), + FE::from(4), + FE::from(5), + FE::from(6), + FE::from(7), + FE::from(8), + FE::from(9), + FE::from(10), + FE::from(11), + FE::from(12), + FE::from(13), + FE::from(14), + FE::from(15), + FE::from(16), + FE::from(17), + FE::from(18), + FE::from(19), + FE::from(20), + FE::from(21), + FE::from(22), + FE::from(23), + FE::from(24), + FE::from(25), + FE::from(26), + FE::from(27), + FE::from(28), + FE::from(29), + FE::from(30), + FE::from(31), + FE::from(32), + ]; + let evals = evaluate_cfft(coeff.clone()); + + let coset = Coset::new_standard(5); + let coset_points = Coset::get_coset_points(&coset); + + assert_eq!(evals[4], evaluate_point(&coeff, &coset_points[4].clone())); + + let new_coeff = interpolate_cfft(evals); + + assert_eq!(coeff, new_coeff); + } +} diff --git a/math/src/circle/twiddles.rs b/math/src/circle/twiddles.rs new file mode 100644 index 000000000..6a07804c4 --- /dev/null +++ b/math/src/circle/twiddles.rs @@ -0,0 +1,83 @@ +extern crate alloc; +use crate::{ + circle::cosets::Coset, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; +use alloc::vec::Vec; + +#[derive(PartialEq)] +pub enum TwiddlesConfig { + Evaluation, + Interpolation, +} +#[cfg(feature = "alloc")] +pub fn get_twiddles( + domain: Coset, + config: TwiddlesConfig, +) -> Vec>> { + // We first take the half coset. + let half_domain_points = Coset::get_coset_points(&Coset::half_coset(domain.clone())); + + // The first set of twiddles are all the y coordinates of the half coset. + let mut twiddles: Vec>> = Vec::new(); + twiddles.push(half_domain_points.iter().map(|p| p.y).collect()); + + if domain.log_2_size >= 2 { + // The second set of twiddles are the x coordinates of the first half of the half coset. + twiddles.push( + half_domain_points + .iter() + .take(half_domain_points.len() / 2) + .map(|p| p.x) + .collect(), + ); + for _ in 0..(domain.log_2_size - 2) { + // The rest of the sets of twiddles are the "square" of the x coordinates of the first half of the previous set. + let prev = twiddles.last().unwrap(); + let cur = prev + .iter() + .take(prev.len() / 2) + .map(|x| x.square().double() - FieldElement::::one()) + .collect(); + twiddles.push(cur); + } + } + + if config == TwiddlesConfig::Interpolation { + // For the interpolation, we need to take the inverse element of each twiddle in the default order. + // We can take inverse being sure that the `unwrap` won't panic because the twiddles are coordinates + // of elements of the coset (or their squares) so they can't be zero. + twiddles.iter_mut().for_each(|x| { + FieldElement::::inplace_batch_inverse(x).unwrap(); + }); + } else { + // For the evaluation, we need reverse the order of the vector of twiddles. + twiddles.reverse(); + } + twiddles +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn evaluation_twiddles_vectors_length_is_correct() { + let domain = Coset::new_standard(3); + let config = TwiddlesConfig::Evaluation; + let twiddles = get_twiddles(domain, config); + for i in 0..twiddles.len() - 1 { + assert_eq!(2 * twiddles[i].len(), twiddles[i + 1].len()) + } + } + + #[test] + fn interpolation_twiddles_vectors_length_is_correct() { + let domain = Coset::new_standard(3); + let config = TwiddlesConfig::Interpolation; + let twiddles = get_twiddles(domain, config); + for i in 0..twiddles.len() - 1 { + assert_eq!(twiddles[i].len(), 2 * twiddles[i + 1].len()) + } + } +} diff --git a/math/src/field/fields/mersenne31/extensions.rs b/math/src/field/fields/mersenne31/extensions.rs index 27c2ab118..69c64f096 100644 --- a/math/src/field/fields/mersenne31/extensions.rs +++ b/math/src/field/fields/mersenne31/extensions.rs @@ -8,6 +8,8 @@ use crate::field::{ use alloc::vec::Vec; type FpE = FieldElement; +type Fp2E = FieldElement; +type Fp4E = FieldElement; #[derive(Clone, Debug)] pub struct Degree2ExtensionField; @@ -132,11 +134,18 @@ impl IsSubFieldOf for Mersenne31Field { } } -type Fp2E = FieldElement; - #[derive(Clone, Debug)] pub struct Degree4ExtensionField; +impl Degree4ExtensionField { + pub const fn const_from_coefficients(a: u32, b: u32, c: u32, d: u32) -> Fp4E { + Fp4E::const_from_raw([ + Fp2E::const_from_raw([FpE::const_from_raw(a), FpE::const_from_raw(b)]), + Fp2E::const_from_raw([FpE::const_from_raw(c), FpE::const_from_raw(d)]), + ]) + } +} + impl IsField for Degree4ExtensionField { type BaseType = [Fp2E; 2]; diff --git a/math/src/field/fields/mersenne31/field.rs b/math/src/field/fields/mersenne31/field.rs index 1c8b2dc58..7654dc29b 100644 --- a/math/src/field/fields/mersenne31/field.rs +++ b/math/src/field/fields/mersenne31/field.rs @@ -1,3 +1,4 @@ +use crate::traits::{AsBytes, ByteConversion}; use crate::{ errors::CreationError, field::{ @@ -203,6 +204,12 @@ impl Display for FieldElement { } } +impl AsBytes for FieldElement { + fn as_bytes(&self) -> alloc::vec::Vec { + self.to_bytes_be() + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/math/src/lib.rs b/math/src/lib.rs index 56c6e598e..1f5ae60d6 100644 --- a/math/src/lib.rs +++ b/math/src/lib.rs @@ -3,6 +3,7 @@ #[cfg(feature = "alloc")] extern crate alloc; +pub mod circle; pub mod cyclic_group; pub mod elliptic_curve; pub mod errors; diff --git a/provers/circle_stark/Cargo.toml b/provers/circle_stark/Cargo.toml new file mode 100644 index 000000000..27a7a54c0 --- /dev/null +++ b/provers/circle_stark/Cargo.toml @@ -0,0 +1,20 @@ +[package] +name = "circle_stark" +version.workspace = true +edition.workspace = true +license.workspace = true +repository.workspace = true + +[dependencies] +lambdaworks-math = { workspace = true, features = [ + "std", + "lambdaworks-serde-binary", +] } +lambdaworks-crypto = { workspace = true, features = ["std", "serde"] } + +thiserror = "1.0.38" +itertools = "0.11.0" +serde = { version = "1.0", features = ["derive"] } + +[dev-dependencies] +test-log = { version = "0.2.11", features = ["log"] } diff --git a/provers/circle_stark/src/air.rs b/provers/circle_stark/src/air.rs new file mode 100644 index 000000000..a59f9b9ba --- /dev/null +++ b/provers/circle_stark/src/air.rs @@ -0,0 +1,95 @@ +use super::{constraints::boundary::BoundaryConstraints, frame::Frame}; +use crate::{ + air_context::AirContext, constraints::transition::TransitionConstraint, domain::Domain, +}; +use lambdaworks_math::{ + circle::point::{CirclePoint, HasCircleParams}, + field::{ + element::FieldElement, + fields::mersenne31::field::Mersenne31Field, + traits::{IsFFTField, IsField, IsSubFieldOf}, + }, +}; +use std::collections::HashMap; +// type ZerofierGroupKey = (usize, usize, Option, Option, usize); +type ZerofierGroupKey = (usize); + +/// AIR is a representation of the Constraints +pub trait AIR { + type PublicInputs; + + fn new(trace_length: usize, pub_inputs: &Self::PublicInputs) -> Self; + + /// Returns the amount trace columns. + fn trace_layout(&self) -> usize; + + fn composition_poly_degree_bound(&self) -> usize; + + /// The method called by the prover to evaluate the transitions corresponding to an evaluation frame. + /// In the case of the prover, the main evaluation table of the frame takes values in + /// `Self::Field`, since they are the evaluations of the main trace at the LDE domain. + fn compute_transition_prover(&self, frame: &Frame) -> Vec> { + let mut evaluations = + vec![FieldElement::::zero(); self.num_transition_constraints()]; + self.transition_constraints() + .iter() + .for_each(|c| c.evaluate(frame, &mut evaluations)); + evaluations + } + + fn boundary_constraints(&self) -> BoundaryConstraints; + + /// The method called by the verifier to evaluate the transitions at the out of domain frame. + /// In the case of the verifier, both main and auxiliary tables of the evaluation frame take + /// values in `Self::FieldExtension`, since they are the evaluations of the trace polynomials + /// at the out of domain challenge. + /// In case `Self::Field` coincides with `Self::FieldExtension`, this method and + /// `compute_transition_prover` should return the same values. + fn compute_transition_verifier(&self, frame: &Frame) -> Vec>; + + fn context(&self) -> &AirContext; + + fn trace_length(&self) -> usize; + + fn blowup_factor(&self) -> u8 { + 2 + } + + fn trace_group_generator(&self) -> CirclePoint { + let trace_length = self.trace_length(); + let log_2_length = trace_length.trailing_zeros(); + CirclePoint::get_generator_of_subgroup(log_2_length) + } + + fn num_transition_constraints(&self) -> usize { + self.context().num_transition_constraints + } + + fn pub_inputs(&self) -> &Self::PublicInputs; + + fn transition_constraints(&self) -> &Vec>; + + fn transition_zerofier_evaluations( + &self, + domain: &Domain, + ) -> Vec>> { + let mut evals = vec![Vec::new(); self.num_transition_constraints()]; + let mut zerofier_groups: HashMap>> = + HashMap::new(); + + self.transition_constraints().iter().for_each(|c| { + let end_exemptions = c.end_exemptions(); + // This hashmap is used to avoid recomputing with an fft the same zerofier evaluation + // If there are multiple domain and subdomains it can be further optimized + // as to share computation between them + let zerofier_group_key = end_exemptions; + zerofier_groups + .entry(zerofier_group_key) + .or_insert_with(|| c.zerofier_evaluations_on_extended_domain(domain)); + let zerofier_evaluations = zerofier_groups.get(&zerofier_group_key).unwrap(); + // println!("ZEROFIER_EVALUATIONS: {:?}", zerofier_evaluations); + evals[c.constraint_idx()] = zerofier_evaluations.clone(); + }); + evals + } +} diff --git a/provers/circle_stark/src/air_context.rs b/provers/circle_stark/src/air_context.rs new file mode 100644 index 000000000..efee83d08 --- /dev/null +++ b/provers/circle_stark/src/air_context.rs @@ -0,0 +1,31 @@ +use std::collections::HashSet; + +#[derive(Clone, Debug)] +pub struct AirContext { + pub trace_columns: usize, + + /// This is a vector with the indices of all the rows that constitute + /// an evaluation frame. Note that, because of how we write all constraints + /// in one method (`compute_transitions`), this vector needs to include the + /// offsets that are needed to compute EVERY transition constraint, even if some + /// constraints don't use all of the indexes in said offsets. + pub transition_offsets: Vec, + pub transition_exemptions: Vec, + pub num_transition_constraints: usize, +} + +impl AirContext { + pub fn num_transition_constraints(&self) -> usize { + self.num_transition_constraints + } + + /// Returns the number of non-trivial different + /// transition exemptions. + pub fn num_transition_exemptions(&self) -> usize { + self.transition_exemptions + .iter() + .filter(|&x| *x != 0) + .collect::>() + .len() + } +} diff --git a/provers/circle_stark/src/composition_polynomial.rs b/provers/circle_stark/src/composition_polynomial.rs new file mode 100644 index 000000000..2704ee336 --- /dev/null +++ b/provers/circle_stark/src/composition_polynomial.rs @@ -0,0 +1,184 @@ +use crate::air::AIR; +use crate::{ + constraints::boundary::BoundaryConstraints, domain::Domain, frame::Frame, trace::LDETraceTable, +}; +use itertools::Itertools; +use lambdaworks_math::{ + circle::{ + point::CirclePoint, + polynomial::{evaluate_cfft, evaluate_point, interpolate_cfft}, + }, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; + +pub(crate) fn evaluate_cp( + air: &A, + lde_trace: &LDETraceTable, + domain: &Domain, + transition_coefficients: &[FieldElement], + boundary_coefficients: &[FieldElement], +) -> Vec> { + // >>> First, we compute the part of the Composition Polynomial related to the boundary constraints. + + let boundary_constraints = &air.boundary_constraints(); + let number_of_b_constraints = boundary_constraints.constraints.len(); + let trace_coset = &domain.trace_coset_points; + let lde_coset = &domain.lde_coset_points; + + // For each pair of boundary constraints, we calculate the denominator's evaluations. + let boundary_zerofiers_inverse_evaluations = + boundary_constraints.evaluate_zerofiers(&trace_coset, &lde_coset); + + // For each pair of boundary constraints, we calculate the numerator's evaluations. + let boundary_polys_evaluations = + boundary_constraints.evaluate_poly_constraints(&trace_coset, &lde_coset, lde_trace); + + // We begin to construct the cp by adding each numerator mulpitlied by the denominator and the beta coefficient. + let cp_boundary: Vec> = (0..lde_coset.len()) + .map(|domain_index| { + (0..number_of_b_constraints) + .step_by(2) + .zip(boundary_coefficients) + .fold(FieldElement::zero(), |acc, (constraint_index, beta)| { + acc + &boundary_zerofiers_inverse_evaluations[constraint_index][domain_index] + * beta + * &boundary_polys_evaluations[constraint_index][domain_index] + }) + }) + .collect(); + + // >>> Now we compute the part of the CP related to the transition constraints and add it to the already + // computed part of the boundary constraints. + + // For each transition constraint, we calulate its zerofier's evaluations. + let transition_zerofiers_inverse_evaluations = air.transition_zerofier_evaluations(domain); + + // + let cp_evaluations = (0..lde_coset.len()) + .zip(cp_boundary) + .map(|(eval_index, boundary_eval)| { + let frame = + Frame::read_from_lde(lde_trace, eval_index, &air.context().transition_offsets); + let transition_poly_evaluations = air.compute_transition_prover(&frame); + let transition_polys_accumulator = itertools::izip!( + transition_poly_evaluations, + &transition_zerofiers_inverse_evaluations, + transition_coefficients + ) + .fold( + FieldElement::zero(), + |acc, (transition_eval, zerof_eval, beta)| { + acc + &zerof_eval[eval_index] * transition_eval * beta + }, + ); + transition_polys_accumulator + boundary_eval + }) + .collect(); + + cp_evaluations +} + +#[cfg(test)] +mod test { + + use crate::examples::simple_fibonacci::{self, FibonacciAIR, FibonacciPublicInputs}; + + use super::*; + + type FE = FieldElement; + + fn build_fibonacci_example() {} + + #[test] + fn boundary_zerofiers_vanish_correctly() { + // Build Fibonacci Example + let trace = simple_fibonacci::fibonacci_trace([FE::one(), FE::one()], 32); + let pub_inputs = FibonacciPublicInputs { + a0: FE::one(), + a1: FE::one(), + }; + let air = FibonacciAIR::new(trace.n_rows(), &pub_inputs); + let boundary_constraints = air.boundary_constraints(); + let domain = Domain::new(&air); + + // Calculate the boundary zerofiers evaluations (L function). + let boundary_zerofiers_inverse_evaluations = boundary_constraints + .evaluate_zerofiers(&domain.trace_coset_points, &domain.lde_coset_points); + + // Interpolate the boundary zerofiers evaluations. + let boundary_zerofiers_coeff = boundary_zerofiers_inverse_evaluations + .iter() + .map(|evals| { + let mut inverse_evals = evals.clone(); + FieldElement::inplace_batch_inverse(&mut inverse_evals).unwrap(); + interpolate_cfft(inverse_evals.to_vec()) + }) + .collect::>>>(); + + // Since simple fibonacci only has one pair of boundary constraints we only check that + // the corresponding polynomial evaluates 0 in the first two coset points and different from 0 + // in the rest of the points. + assert_eq!( + evaluate_point(&boundary_zerofiers_coeff[0], &domain.trace_coset_points[0]), + FE::zero() + ); + + assert_eq!( + evaluate_point(&boundary_zerofiers_coeff[0], &domain.trace_coset_points[1]), + FE::zero() + ); + + for point in domain.trace_coset_points.iter().skip(2) { + assert_ne!( + evaluate_point(&boundary_zerofiers_coeff[0], &point), + FE::zero() + ); + } + } + + #[test] + fn boundary_polys_vanish_correctly() { + // Build Fibonacci Example + let trace = simple_fibonacci::fibonacci_trace([FE::one(), FE::one()], 32); + let pub_inputs = FibonacciPublicInputs { + a0: FE::one(), + a1: FE::one(), + }; + let air = FibonacciAIR::new(trace.n_rows(), &pub_inputs); + let boundary_constraints = air.boundary_constraints(); + let domain = Domain::new(&air); + + // Evaluate each polynomial in the lde domain. + let lde_trace = LDETraceTable::new(trace.table.data.clone(), 1, 1); + + // Calculate boundary polynomials evaluations (the polynomial f - I). + let boundary_polys_evaluations = boundary_constraints.evaluate_poly_constraints( + &domain.trace_coset_points, + &domain.lde_coset_points, + &lde_trace, + ); + + // Interpolate the boundary polynomials evaluations. + let boundary_poly_coeff = boundary_polys_evaluations + .iter() + .map(|evals| interpolate_cfft(evals.to_vec())) + .collect::>>>(); + + // Since simple fibonacci only has one pair of boundary constraints we only check that + // the corresponding polynomial evaluates 0 in the first two coset points and different from 0 + // in the rest of the points. + assert_eq!( + evaluate_point(&boundary_poly_coeff[0], &domain.trace_coset_points[0]), + FE::zero() + ); + + assert_eq!( + evaluate_point(&boundary_poly_coeff[0], &domain.trace_coset_points[1]), + FE::zero() + ); + + for point in domain.trace_coset_points.iter().skip(2) { + assert_ne!(evaluate_point(&boundary_poly_coeff[0], &point), FE::zero()); + } + } +} diff --git a/provers/circle_stark/src/config.rs b/provers/circle_stark/src/config.rs new file mode 100644 index 000000000..52df5fcef --- /dev/null +++ b/provers/circle_stark/src/config.rs @@ -0,0 +1,17 @@ +use lambdaworks_crypto::merkle_tree::{ + backends::types::{BatchKeccak256Backend, Keccak256Backend}, + merkle::MerkleTree, +}; + +// Merkle Trees configuration + +// Security of both hashes should match + +pub type FriMerkleTreeBackend = Keccak256Backend; +pub type FriMerkleTree = MerkleTree>; + +pub const COMMITMENT_SIZE: usize = 32; +pub type Commitment = [u8; COMMITMENT_SIZE]; + +pub type BatchedMerkleTreeBackend = BatchKeccak256Backend; +pub type BatchedMerkleTree = MerkleTree>; diff --git a/provers/circle_stark/src/constraints/boundary.rs b/provers/circle_stark/src/constraints/boundary.rs new file mode 100644 index 000000000..dab5a92e1 --- /dev/null +++ b/provers/circle_stark/src/constraints/boundary.rs @@ -0,0 +1,282 @@ +use itertools::Itertools; +use lambdaworks_math::{ + circle::point::CirclePoint, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field, traits::IsField}, + polynomial::Polynomial, +}; + +use crate::trace::LDETraceTable; + +use super::{evaluator::line, utils::interpolant}; + +#[derive(Debug)] +/// Represents a boundary constraint that must hold in an execution +/// trace: +/// * col: The column of the trace where the constraint must hold +/// * step: The step (or row) of the trace where the constraint must hold +/// * value: The value the constraint must have in that column and step +pub struct BoundaryConstraint { + pub col: usize, + pub step: usize, + pub value: FieldElement, +} + +impl BoundaryConstraint { + pub fn new(col: usize, step: usize, value: FieldElement) -> Self { + Self { col, step, value } + } + + /// Used for creating boundary constraints for a trace with only one column + pub fn new_simple(step: usize, value: FieldElement) -> Self { + Self { + col: 0, + step, + value, + } + } +} + +/// Data structure that stores all the boundary constraints that must +/// hold for the execution trace +#[derive(Default, Debug)] +pub struct BoundaryConstraints { + pub constraints: Vec, +} + +impl BoundaryConstraints { + #[allow(dead_code)] + pub fn new() -> Self { + Self { + constraints: Vec::::new(), + } + } + + /// To instantiate from a vector of BoundaryConstraint elements + pub fn from_constraints(constraints: Vec) -> Self { + Self { constraints } + } + + /// Returns all the steps where boundary conditions exist for the given column + pub fn steps(&self, col: usize) -> Vec { + self.constraints + .iter() + .filter(|v| v.col == col) + .map(|c| c.step) + .collect() + } + + /// Return all the steps where boundary constraints hold. + pub fn steps_for_boundary(&self) -> Vec { + self.constraints + .iter() + .unique_by(|elem| elem.step) + .map(|v| v.step) + .collect() + } + + /// Return all the columns where boundary constraints hold. + pub fn cols_for_boundary(&self) -> Vec { + self.constraints + .iter() + .unique_by(|elem| elem.col) + .map(|v| v.col) + .collect() + } + + /// Given the group generator of some domain, returns for each column the domain values corresponding + /// to the steps where the boundary conditions hold. This is useful when interpolating + /// the boundary conditions, since we must know the x values + pub fn generate_roots_of_unity( + &self, + group_generator: &CirclePoint, + cols_trace: &[usize], + ) -> Vec>> { + cols_trace + .iter() + .map(|i| { + self.steps(*i) + .into_iter() + .map(|s| group_generator * (s as u128)) + .collect::>>() + }) + .collect::>>>() + } + + /// For every trace column, give all the values the trace must be equal to in + /// the steps where the boundary constraints hold + pub fn values(&self, cols_trace: &[usize]) -> Vec>> { + cols_trace + .iter() + .map(|i| { + self.constraints + .iter() + .filter(|c| c.col == *i) + .map(|c| c.value.clone()) + .collect() + }) + .collect() + } + + /// For every column, it returns for each pair of boundary constraints in that column, the corresponding evaluation + /// in all the `eval_points`. + // We assume that there are an even number of boundary contrainsts in the column `col`. + // TODO: If two columns have a pair of boundary constraints that hold in the same rows, we can optimize + // this function so that we don't calculate the same line evaluations for the two of them. Maybe a hashmap? + pub fn evaluate_zerofiers( + &self, + trace_coset: &Vec>, + eval_points: &Vec>, + ) -> Vec>> { + let mut zerofiers_evaluations = Vec::new(); + for col in self.cols_for_boundary() { + self.constraints + .iter() + .filter(|constraint| constraint.col == col) + .chunks(2) + .into_iter() + .map(|chunk| { + let chunk: Vec<_> = chunk.collect(); + let first_constraint = &chunk[0]; + let second_constraint = &chunk[1]; + let first_vanish_point = &trace_coset[first_constraint.step]; + let second_vanish_point = &trace_coset[second_constraint.step]; + + let mut boundary_evaluations: Vec> = eval_points + .iter() + .map(|eval_point| { + line(eval_point, &first_vanish_point, &second_vanish_point) + }) + .collect(); + FieldElement::inplace_batch_inverse(&mut boundary_evaluations).unwrap(); + zerofiers_evaluations.push(boundary_evaluations); + }) + .collect() // TODO: We don't use this collect. + } + zerofiers_evaluations + } + + /// For every columnm, and every constrain, returns the evaluation on the entire lde domain of the constrain function F(x) - I(x). + // Note this is the numerator of each pair of constrains for the composition polynomial. + pub fn evaluate_poly_constraints( + &self, + trace_coset: &Vec>, + eval_points: &Vec>, + lde_trace: &LDETraceTable, + ) -> Vec>> { + let mut poly_evaluations = Vec::new(); + for col in self.cols_for_boundary() { + self.constraints + .iter() + .filter(|constraint| constraint.col == col) + .chunks(2) + .into_iter() + .map(|chunk| { + let chunk: Vec<_> = chunk.collect(); + let first_constraint = &chunk[0]; + let second_constraint = &chunk[1]; + let first_vanish_point = &trace_coset[first_constraint.step]; + let first_value = first_constraint.value; + let second_vanish_point = &trace_coset[second_constraint.step]; + let second_value = second_constraint.value; + let boundary_evaluations = eval_points + .iter() + .zip(&lde_trace.table.columns()[col]) + .map(|(eval_point, lde_eval)| { + lde_eval + - interpolant( + &first_vanish_point, + &second_vanish_point, + first_value, + second_value, + eval_point, + ) + }) + .collect::>>(); + poly_evaluations.push(boundary_evaluations); + }) + .collect() // TODO: We are not using this collect. + } + poly_evaluations + } + + // /// For every columnm, and every constrain, returns the evaluation on the entire lde domain of the constrain function F(x) - I(x). + // // Note this is the numerator of each pair of constrains for the composition polynomial. + // pub fn evaluate_poly_constraints( + // &self, + // trace_coset: &Vec>, + // eval_points: &Vec>, + // lde_trace: &LDETraceTable, + // ) -> Vec>>> { + // self.cols_for_boundary() + // .iter() + // .map(|col| { + // self.constraints + // .iter() + // .filter(|constraint| constraint.col == *col) + // .chunks(2) + // .into_iter() + // .map(|chunk| { + // let chunk: Vec<_> = chunk.collect(); + // let first_constraint = &chunk[0]; + // let second_constraint = &chunk[1]; + // let first_vanish_point = &trace_coset[first_constraint.step]; + // let first_value = first_constraint.value; + // let second_vanish_point = &trace_coset[second_constraint.step]; + // let second_value = second_constraint.value; + // eval_points + // .iter() + // .zip(&lde_trace.table.columns()[*col]) + // .map(|(eval_point, lde_eval)| { + // lde_eval + // - interpolant( + // &first_vanish_point, + // &second_vanish_point, + // first_value, + // second_value, + // eval_point, + // ) + // }) + // .collect::>>() + // }) + // .collect() + // }) + // .collect() + // } +} + +#[cfg(test)] +mod test { + use lambdaworks_math::{ + circle::cosets::Coset, + field::{ + fields::fft_friendly::stark_252_prime_field::Stark252PrimeField, traits::IsFFTField, + }, + }; + type PrimeField = Stark252PrimeField; + + use super::*; + + #[test] + fn simple_fibonacci_boundary_zerofiers() { + let one = FieldElement::::one(); + + // Fibonacci constraints: + // a0 = 1 + // a1 = 1 + let a0 = BoundaryConstraint::new_simple(0, one); + let a1 = BoundaryConstraint::new_simple(1, one); + + let trace_coset = Coset::get_coset_points(&Coset::new_standard(3)); + let eval_points = Coset::get_coset_points(&Coset::new_standard(4)); + + let expected_zerofier: Vec>> = vec![eval_points + .iter() + .map(|point| line(point, &trace_coset[0], &trace_coset[1])) + .collect()]; + + let constraints = BoundaryConstraints::from_constraints(vec![a0, a1]); + let zerofier = constraints.evaluate_zerofiers(&trace_coset, &eval_points); + + assert_eq!(expected_zerofier, zerofier); + } +} diff --git a/provers/circle_stark/src/constraints/evaluator.rs b/provers/circle_stark/src/constraints/evaluator.rs new file mode 100644 index 000000000..956071461 --- /dev/null +++ b/provers/circle_stark/src/constraints/evaluator.rs @@ -0,0 +1,258 @@ +use super::boundary::BoundaryConstraints; +use crate::air::AIR; +use std::marker::PhantomData; + +use crate::{domain::Domain, frame::Frame, trace::LDETraceTable}; +use itertools::Itertools; +use lambdaworks_math::circle::point::CirclePoint; +use lambdaworks_math::circle::polynomial::{evaluate_point, interpolate_cfft}; +use lambdaworks_math::field::element::FieldElement; +use lambdaworks_math::field::fields::mersenne31::field::Mersenne31Field; + +pub struct ConstraintEvaluator { + boundary_constraints: BoundaryConstraints, + phantom: PhantomData, +} + +// See https://vitalik.eth.limo/general/2024/07/23/circlestarks.html (Section: Quetienting). +// https://github.com/ethereum/research/blob/master/circlestark/line_functions.py#L10 +pub fn line( + eval_point: &CirclePoint, + first_vanish_point: &CirclePoint, + second_vanish_point: &CirclePoint, +) -> FieldElement { + (first_vanish_point.y - second_vanish_point.y) * eval_point.x + + (second_vanish_point.x - first_vanish_point.x) * eval_point.y + + (first_vanish_point.x * second_vanish_point.y + - first_vanish_point.y * second_vanish_point.x) +} + +// See https://vitalik.eth.limo/general/2024/07/23/circlestarks.html (Section: Quetienting). +// https://github.com/ethereum/research/blob/master/circlestark/line_functions.py#L16 +// Evaluates the polybomial I at eval_point. I is the polynomial such that I(point_1) = value_1 and +// I(point_2) = value_2. +pub fn interpolant( + point_1: &CirclePoint, + point_2: &CirclePoint, + value_1: FieldElement, + value_2: FieldElement, + eval_point: &CirclePoint, +) -> FieldElement { + let dx = point_2.x - point_1.x; + let dy = point_2.y - point_1.y; + // CHECK: can dx^2 + dy^2 = 0 even if dx!=0 and dy!=0 ? (using that they are FE of Mersenne31). + let invdist = (dx * dx + dy * dy).inv().unwrap(); + let dot = (eval_point.x - point_1.x) * dx + (eval_point.y - point_1.y) * dy; + value_1 + (value_2 - value_1) * dot * invdist +} + +impl ConstraintEvaluator { + pub fn new(air: &A) -> Self { + let boundary_constraints = air.boundary_constraints(); + + Self { + boundary_constraints, + phantom: PhantomData, + } + } + + pub(crate) fn evaluate( + &self, + air: &A, + lde_trace: &LDETraceTable, + domain: &Domain, + transition_coefficients: &[FieldElement], + boundary_coefficients: &[FieldElement], + ) -> Vec> { + let boundary_constraints = &self.boundary_constraints; + let number_of_b_constraints = boundary_constraints.constraints.len(); + + let boundary_zerofiers_inverse_evaluations: Vec>> = + boundary_constraints + .constraints + .chunks(2) + .map(|chunk| { + let first_constraint = &chunk[0]; + let second_constraint = &chunk[1]; + let first_vanish_point = &domain.trace_coset_points[first_constraint.step]; + let second_vanish_point = &domain.trace_coset_points[second_constraint.step]; + let mut evals = domain + .lde_coset_points + .iter() + .map(|eval_point| line(eval_point, first_vanish_point, second_vanish_point)) + .collect::>>(); + FieldElement::inplace_batch_inverse(&mut evals).unwrap(); + evals + }) + .collect::>>>(); + + let boundary_polys_evaluations: Vec>> = + boundary_constraints + .constraints + .chunks(2) + .map(|chunk| { + let first_constraint = &chunk[0]; + let second_constraint = &chunk[1]; + let first_vanish_point = &domain.trace_coset_points[first_constraint.step]; + let first_value = first_constraint.value; + let second_vanish_point = &domain.trace_coset_points[second_constraint.step]; + let second_value = second_constraint.value; + let evals = domain + .lde_coset_points + .iter() + .zip(&lde_trace.table.data) + .map(|(eval_point, lde_eval)| { + lde_eval - interpolant( + first_vanish_point, + second_vanish_point, + first_value, + second_value, + eval_point, + ) + }) + .collect::>>(); + evals + }) + .collect::>>>(); + + // let boundary_polys_evaluations = boundary_constraints + // .constraints + // .iter() + // .map(|constraint| { + // (0..lde_trace.num_rows()) + // .map(|row| { + // let v = lde_trace.table.get(row, constraint.col); + // v - &constraint.value + // }) + // .collect_vec() + // }) + // .collect_vec(); + + // --------------- BEGIN TESTING ---------------------------- + // Interpolate lde trace evaluations. + let l_poly_coefficients = boundary_zerofiers_inverse_evaluations + .iter() + .map(|evals| { + let mut inverse_evals = evals.clone(); + FieldElement::inplace_batch_inverse(&mut inverse_evals).unwrap(); + interpolate_cfft(inverse_evals.to_vec()) + }) + .collect::>>>(); + + let fi_poly_coefficients = boundary_polys_evaluations + .iter() + .map(|evals| { + interpolate_cfft(evals.to_vec()) + }) + .collect::>>>(); + + // Evaluate lde trace interpolating polynomial in trace domain. + for point in &domain.trace_coset_points { + println!("-----------------------"); + println!( + "L evaluation: {:?}", + evaluate_point(&l_poly_coefficients[0], &point) + ); + println!( + "F-I evaluation: {:?}", + evaluate_point(&fi_poly_coefficients[0], &point) + ); + } + + // --------------- END TESTING ---------------------------- + + let boundary_eval_iter = 0..domain.lde_coset_points.len(); + + let boundary_evaluation: Vec<_> = boundary_eval_iter + .map(|domain_index| { + (0..number_of_b_constraints) + .step_by(2) + .zip(boundary_coefficients) + .fold(FieldElement::zero(), |acc, (constraint_index, beta)| { + acc + &boundary_zerofiers_inverse_evaluations[constraint_index] + [domain_index] + * beta + * &boundary_polys_evaluations[constraint_index][domain_index] + }) + }) + .collect(); + + // Iterate over all LDE domain and compute the part of the composition polynomial + // related to the transition constraints and add it to the already computed part of the + // boundary constraints. + + let zerofiers_evals = air.transition_zerofier_evaluations(domain); + + // --------------- BEGIN TESTING ---------------------------- + // Interpolate lde trace evaluations. + // let zerofier_poly_coefficients = zerofiers_evals + // .iter() + // .map(|evals| interpolate_cfft(evals.to_vec())) + // .collect::>>>(); + + // // Evaluate lde trace interpolating polynomial in trace domain. + // // This should print all zeroes except in the end exceptions points. + // for point in &domain.trace_coset_points { + // println!( + // "{:?}", + // evaluate_point(&zerofier_poly_coefficients[0], &point) + // ); + // } + // --------------- END TESTING ---------------------------- + + let evaluations_t_iter = 0..domain.lde_coset_points.len(); + + let evaluations_t = evaluations_t_iter + .zip(boundary_evaluation) + .map(|(i, boundary)| { + let frame = Frame::read_from_lde(lde_trace, i, &air.context().transition_offsets); + + // Compute all the transition constraints at this point of the LDE domain. + let evaluations_transition = air.compute_transition_prover(&frame); + + // Add each term of the transition constraints to the composition polynomial, including the zerofier, + // the challenge and the exemption polynomial if it is necessary. + let acc_transition = itertools::izip!( + evaluations_transition, + &zerofiers_evals, + transition_coefficients + ) + .fold(FieldElement::zero(), |acc, (eval, zerof_eval, beta)| { + acc + &zerof_eval[i] * eval * beta + }); + + acc_transition + boundary + }) + .collect(); + + evaluations_t + } +} + +#[cfg(test)] +mod tests { + use super::*; + use lambdaworks_math::circle::cosets::Coset; + + type FE = FieldElement; + + #[test] + fn line_vanishes_in_vanishing_points(){ + let first_vanish_point = CirclePoint::GENERATOR * 3; + let second_vanish_point = CirclePoint::GENERATOR * 5; + assert_eq!(line(&first_vanish_point, &first_vanish_point, &second_vanish_point), FE::zero()); + assert_eq!(line(&second_vanish_point, &first_vanish_point, &second_vanish_point), FE::zero()); + } + + #[test] + fn interpolant_takes_the_corresponding_values(){ + let point_1 = CirclePoint::GENERATOR * 4; + let point_2 = CirclePoint::GENERATOR * 9; + let value_1 = FE::from(5); + let value_2 = FE::from(3); + assert_eq!(interpolant(&point_1, &point_2, value_1, value_2, &point_1), value_1); + assert_eq!(interpolant(&point_1, &point_2, value_1, value_2, &point_2), value_2); + + } + +} diff --git a/provers/circle_stark/src/constraints/mod.rs b/provers/circle_stark/src/constraints/mod.rs new file mode 100644 index 000000000..f6f474a44 --- /dev/null +++ b/provers/circle_stark/src/constraints/mod.rs @@ -0,0 +1,4 @@ +pub mod boundary; +pub mod evaluator; +pub mod transition; +pub mod utils; diff --git a/provers/circle_stark/src/constraints/transition.rs b/provers/circle_stark/src/constraints/transition.rs new file mode 100644 index 000000000..08ba957ff --- /dev/null +++ b/provers/circle_stark/src/constraints/transition.rs @@ -0,0 +1,188 @@ +use crate::domain::Domain; +use crate::frame::Frame; +use crate::constraints::evaluator::line; +use lambdaworks_math::circle::point::CirclePoint; +use lambdaworks_math::circle::polynomial::{evaluate_point, interpolate_cfft}; +use lambdaworks_math::field::element::FieldElement; +use lambdaworks_math::field::fields::mersenne31::field::Mersenne31Field; +/// TransitionConstraint represents the behaviour that a transition constraint +/// over the computation that wants to be proven must comply with. +pub trait TransitionConstraint { + /// The degree of the constraint interpreting it as a multivariate polynomial. + fn degree(&self) -> usize; + + /// The index of the constraint. + /// Each transition constraint should have one index in the range [0, N), + /// where N is the total number of transition constraints. + fn constraint_idx(&self) -> usize; + + /// The function representing the evaluation of the constraint over elements + /// of the trace table. + /// + /// Elements of the trace table are found in the `frame` input, and depending on the + /// constraint, elements of `periodic_values` and `rap_challenges` may be used in + /// the evaluation. + /// Once computed, the evaluation should be inserted in the `transition_evaluations` + /// vector, in the index corresponding to the constraint as given by `constraint_idx()`. + fn evaluate(&self, frame: &Frame, transition_evaluations: &mut [FieldElement]); + + /// The periodicity the constraint is applied over the trace. + /// + /// Default value is 1, meaning that the constraint is applied to every + /// step of the trace. + fn period(&self) -> usize { + 1 + } + + /// The offset with respect to the first trace row, where the constraint + /// is applied. + /// For example, if the constraint has periodicity 2 and offset 1, this means + /// the constraint will be applied over trace rows of index 1, 3, 5, etc. + /// + /// Default value is 0, meaning that the constraint is applied from the first + /// element of the trace on. + fn offset(&self) -> usize { + 0 + } + + /// For a more fine-grained description of where the constraint should apply, + /// an exemptions period can be defined. + /// This specifies the periodicity of the row indexes where the constraint should + /// NOT apply, within the row indexes where the constraint applies, as specified by + /// `period()` and `offset()`. + /// + /// Default value is None. + fn exemptions_period(&self) -> Option { + None + } + + /// The offset value for periodic exemptions. Check documentation of `period()`, + /// `offset()` and `exemptions_period` for a better understanding. + fn periodic_exemptions_offset(&self) -> Option { + None + } + + /// The number of exemptions at the end of the trace. + /// + /// This method's output defines what trace elements should not be considered for + /// the constraint evaluation at the end of the trace. For example, for a fibonacci + /// computation that has to use the result 2 following steps, this method is defined + /// to return the value 2. + fn end_exemptions(&self) -> usize; + + /// Evaluate the `eval_point` in the polynomial that vanishes in all the exemptions points. + fn evaluate_end_exemptions_poly( + &self, + eval_point: &CirclePoint, + // `trace_group_generator` can be calculated with `trace_length` but it is better to precompute it + trace_group_generator: &CirclePoint, + trace_length: usize, + ) -> FieldElement { + + let one = FieldElement::::one(); + + if self.end_exemptions() == 0 { + return one; + } + + let double_group_generator = CirclePoint::::get_generator_of_subgroup( + trace_length.trailing_zeros() + 1, + ); + + (1..=self.end_exemptions()) + .step_by(2) + .map(|exemption| { + println!("EXEMPTION: {:?}", exemption); + let first_vanish_point = &double_group_generator + (trace_group_generator * ((trace_length - exemption) as u128)); + + let second_vanish_point = &double_group_generator + (trace_group_generator * ((trace_length - (exemption + 1)) as u128)); + + line(eval_point, &first_vanish_point, &second_vanish_point) + }) + .fold(one, |acc, eval| acc * eval) + } + + /// Compute evaluations of the constraints zerofier over a LDE domain. + /// TODO: See if we can evaluate using cfft. + /// TODO: See if we can optimize computing only some evaluations and cycle them as in regular stark. + #[allow(unstable_name_collisions)] + fn zerofier_evaluations_on_extended_domain( + &self, + domain: &Domain, + ) -> Vec> { + let blowup_factor = domain.blowup_factor; + let trace_length = domain.trace_length; + let trace_log_2_size = trace_length.trailing_zeros(); + let lde_log_2_size = (blowup_factor * trace_length).trailing_zeros(); + let trace_group_generator = &domain.trace_group_generator; + + // if let Some(exemptions_period) = self.exemptions_period() { + + // } else { + + let lde_points = &domain.lde_coset_points; + + let mut zerofier_evaluations: Vec<_> = lde_points + .iter() + .map(|point| { + // TODO: Is there a way to avoid this clone()? + let mut x = point.x.clone(); + for _ in 1..trace_log_2_size { + x = x.square().double() - FieldElement::::one(); + } + x + }) + .collect(); + FieldElement::inplace_batch_inverse(&mut zerofier_evaluations).unwrap(); + + let end_exemptions_evaluations: Vec<_> = lde_points + .iter() + .map(|point| { + self.evaluate_end_exemptions_poly(point, trace_group_generator, trace_length) + }) + .collect(); + + // --------------- BEGIN TESTING ---------------------------- + // Interpolate lde trace evaluations. + let end_exemptions_coeff = interpolate_cfft(end_exemptions_evaluations.clone()); + + // Evaluate lde trace interpolating polynomial in trace domain. + // This should print zeroes only in the end exceptions points. + for point in &domain.trace_coset_points { + println!( + "EXEMPTIONS POLYS EVALUATED ON TRACE DOMAIN {:?}", + evaluate_point(&end_exemptions_coeff, &point) + ); + } + // --------------- END TESTING ---------------------------- + + std::iter::zip(zerofier_evaluations, end_exemptions_evaluations) + .map(|(eval, exemptions_eval)| eval * exemptions_eval) + .collect() + } + + ///// Returns the evaluation of the zerofier corresponding to this constraint in some point + ///// `eval_point`, (which is in the circle over the extension field). + // #[allow(unstable_name_collisions)] + // fn evaluate_zerofier( + // &self, + // eval_point: &CirclePoint, + // trace_group_generator: &CirclePoint, + // trace_length: usize, + // ) -> FieldElement { + // // if let Some(exemptions_period) = self.exemptions_period() { + + // // } else { + + // let end_exemptions_evaluation = + // self.evaluate_end_exemptions_poly(eval_point, trace_group_generator, trace_length); + + // let trace_log_2_size = trace_length.trailing_zeros(); + // let mut x = eval_point.x.clone(); + // for _ in 1..trace_log_2_size { + // x = x.square().double() - FieldElement::::one(); + // } + + // x.inv().unwrap() * end_exemptions_evaluation + // } +} diff --git a/provers/circle_stark/src/constraints/utils.rs b/provers/circle_stark/src/constraints/utils.rs new file mode 100644 index 000000000..b92ddde6b --- /dev/null +++ b/provers/circle_stark/src/constraints/utils.rs @@ -0,0 +1,36 @@ +use lambdaworks_math::{ + circle::point::CirclePoint, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; + +// See https://vitalik.eth.limo/general/2024/07/23/circlestarks.html (Section: Quetienting). +// https://github.com/ethereum/research/blob/master/circlestark/line_functions.py#L10 +pub fn line( + eval_point: &CirclePoint, + first_vanish_point: &CirclePoint, + second_vanish_point: &CirclePoint, +) -> FieldElement { + (first_vanish_point.y - second_vanish_point.y) * eval_point.x + + (second_vanish_point.x - first_vanish_point.x) * eval_point.y + + (first_vanish_point.x * second_vanish_point.y + - first_vanish_point.y * second_vanish_point.x) +} + +// See https://vitalik.eth.limo/general/2024/07/23/circlestarks.html (Section: Quetienting). +// https://github.com/ethereum/research/blob/master/circlestark/line_functions.py#L16 +// Evaluates the polybomial I at eval_point. I is the polynomial such that I(point_1) = value_1 and +// I(point_2) = value_2. +pub fn interpolant( + point_1: &CirclePoint, + point_2: &CirclePoint, + value_1: FieldElement, + value_2: FieldElement, + eval_point: &CirclePoint, +) -> FieldElement { + let dx = point_2.x - point_1.x; + let dy = point_2.y - point_1.y; + // CHECK: can dx^2 + dy^2 = 0 even if dx!=0 and dy!=0 ? (using that they are FE of Mersenne31). + let invdist = (dx * dx + dy * dy).inv().unwrap(); + let dot = (eval_point.x - point_1.x) * dx + (eval_point.y - point_1.y) * dy; + value_1 + (value_2 - value_1) * dot * invdist +} diff --git a/provers/circle_stark/src/domain.rs b/provers/circle_stark/src/domain.rs new file mode 100644 index 000000000..a684fddb7 --- /dev/null +++ b/provers/circle_stark/src/domain.rs @@ -0,0 +1,46 @@ +use lambdaworks_math::{ + circle::{ + cosets::Coset, + point::CirclePoint, + }, + field::fields::mersenne31::field::Mersenne31Field, +}; + +use super::air::AIR; + +pub struct Domain { + pub(crate) trace_length: usize, + pub(crate) trace_log_2_length: u32, + pub(crate) blowup_factor: usize, + pub(crate) trace_coset_points: Vec>, + pub(crate) lde_coset_points: Vec>, + pub(crate) trace_group_generator: CirclePoint, +} + +impl Domain { + pub fn new(air: &A) -> Self + where + A: AIR, + { + // Initial definitions + let trace_length = air.trace_length(); + let trace_log_2_length = trace_length.trailing_zeros(); + let blowup_factor = air.blowup_factor() as usize; + + // * Generate Coset + let trace_coset_points = Coset::get_coset_points(&Coset::new_standard(trace_log_2_length)); + let lde_coset_points = Coset::get_coset_points(&Coset::new_standard( + (blowup_factor * trace_length).trailing_zeros(), + )); + let trace_group_generator = CirclePoint::get_generator_of_subgroup(trace_log_2_length); + + Self { + trace_length, + trace_log_2_length, + blowup_factor, + trace_coset_points, + lde_coset_points, + trace_group_generator, + } + } +} diff --git a/provers/circle_stark/src/examples/mod.rs b/provers/circle_stark/src/examples/mod.rs new file mode 100644 index 000000000..c8bec26b5 --- /dev/null +++ b/provers/circle_stark/src/examples/mod.rs @@ -0,0 +1 @@ +pub mod simple_fibonacci; \ No newline at end of file diff --git a/provers/circle_stark/src/examples/simple_fibonacci.rs b/provers/circle_stark/src/examples/simple_fibonacci.rs new file mode 100644 index 000000000..927bf13a4 --- /dev/null +++ b/provers/circle_stark/src/examples/simple_fibonacci.rs @@ -0,0 +1,117 @@ +use crate::{ + constraints::{ + boundary::{BoundaryConstraint, BoundaryConstraints}, + transition::TransitionConstraint, + }, + air_context::AirContext, + frame::Frame, + trace::TraceTable, + air::AIR, +}; +use lambdaworks_math::field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field, traits::IsFFTField}; +// use std::marker::PhantomData; +#[derive(Clone)] +struct FibConstraint; + +impl FibConstraint { + pub fn new() -> Self { + Self {} + } +} +impl TransitionConstraint for FibConstraint { + fn degree(&self) -> usize { + 1 + } + fn constraint_idx(&self) -> usize { + 0 + } + fn end_exemptions(&self) -> usize { + 4 + } + fn evaluate( + &self, + frame: &Frame, + transition_evaluations: &mut [FieldElement], + ) { + let first_step = frame.get_evaluation_step(0); + let second_step = frame.get_evaluation_step(1); + let third_step = frame.get_evaluation_step(2); + let a0 = first_step[0]; + let a1 = second_step[0]; + let a2 = third_step[0]; + let res = a2 - a1 - a0; + transition_evaluations[self.constraint_idx()] = res; + } +} +pub struct FibonacciAIR { + context: AirContext, + trace_length: usize, + pub_inputs: FibonacciPublicInputs, + constraints: Vec>, +} +#[derive(Clone, Debug)] +pub struct FibonacciPublicInputs{ + pub a0: FieldElement, + pub a1: FieldElement, +} +impl AIR for FibonacciAIR +{ + type PublicInputs = FibonacciPublicInputs; + fn new(trace_length: usize, pub_inputs: &Self::PublicInputs) -> Self { + let constraints: Vec> = + vec![Box::new(FibConstraint::new())]; + let context = AirContext { + trace_columns: 1, + transition_exemptions: vec![2], + transition_offsets: vec![0, 1, 2], + num_transition_constraints: constraints.len(), + }; + Self { + pub_inputs: pub_inputs.clone(), + context, + trace_length, + constraints, + } + } + fn composition_poly_degree_bound(&self) -> usize { + self.trace_length() + } + fn transition_constraints(&self) -> &Vec> { + &self.constraints + } + fn boundary_constraints(&self) -> BoundaryConstraints { + let a0 = BoundaryConstraint::new_simple(0, self.pub_inputs.a0.clone()); + let a1 = BoundaryConstraint::new_simple(1, self.pub_inputs.a1.clone()); + BoundaryConstraints::from_constraints(vec![a0, a1]) + } + fn context(&self) -> &AirContext { + &self.context + } + fn trace_length(&self) -> usize { + self.trace_length + } + fn trace_layout(&self) -> usize { + 1 + } + fn pub_inputs(&self) -> &Self::PublicInputs { + &self.pub_inputs + } + fn compute_transition_verifier( + &self, + frame: &Frame, + ) -> Vec> { + self.compute_transition_prover(frame) + } +} +pub fn fibonacci_trace( + initial_values: [FieldElement; 2], + trace_length: usize, +) -> TraceTable { + let mut ret: Vec> = vec![]; + ret.push(initial_values[0].clone()); + ret.push(initial_values[1].clone()); + for i in 2..(trace_length) { + ret.push(ret[i - 1].clone() + ret[i - 2].clone()); + } + TraceTable::from_columns(vec![ret]) +} diff --git a/provers/circle_stark/src/frame.rs b/provers/circle_stark/src/frame.rs new file mode 100644 index 000000000..315190f5f --- /dev/null +++ b/provers/circle_stark/src/frame.rs @@ -0,0 +1,32 @@ +use crate::trace::LDETraceTable; +use itertools::Itertools; +use lambdaworks_math::field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}; + +/// A frame represents a collection of trace steps. +/// The collected steps are all the necessary steps for +/// all transition costraints over a trace to be evaluated. +#[derive(Clone, Debug, PartialEq)] +pub struct Frame { + steps: Vec>>, +} + +impl Frame { + pub fn new(steps: Vec>>) -> Self { + Self { steps } + } + + pub fn get_evaluation_step(&self, step: usize) -> &Vec> { + &self.steps[step] + } + + pub fn read_from_lde(lde_trace: &LDETraceTable, row: usize, offsets: &[usize]) -> Self { + let num_rows = lde_trace.num_rows(); + + let lde_steps = offsets + .iter() + .map(|offset| lde_trace.get_row((row + offset) % num_rows).to_vec()) + .collect_vec(); + + Frame::new(lde_steps) + } +} diff --git a/provers/circle_stark/src/lib.rs b/provers/circle_stark/src/lib.rs new file mode 100644 index 000000000..999534a5e --- /dev/null +++ b/provers/circle_stark/src/lib.rs @@ -0,0 +1,15 @@ +pub mod air; +pub mod air_context; +pub mod composition_polynomial; +pub mod config; +pub mod constraints; +pub mod domain; +pub mod examples; +pub mod frame; +pub mod prover; +pub mod table; +pub mod trace; +pub mod vanishing_poly; + +#[cfg(test)] +pub mod tests; diff --git a/provers/circle_stark/src/prover.rs b/provers/circle_stark/src/prover.rs new file mode 100644 index 000000000..6032aaab2 --- /dev/null +++ b/provers/circle_stark/src/prover.rs @@ -0,0 +1,143 @@ +use crate::{air::AIR, config::FriMerkleTree, constraints::evaluator::ConstraintEvaluator, domain::Domain, frame::Frame, trace::{LDETraceTable, TraceTable}}; +use std::marker::PhantomData; + +use super::config::Commitment; +use lambdaworks_math::{ + circle::polynomial::{evaluate_cfft, evaluate_point, interpolate_cfft}, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; + +/// A default STARK prover implementing `IsStarkProver`. +pub struct Prover { + phantom: PhantomData, +} + +impl IsStarkProver for Prover {} + +#[derive(Debug)] +pub enum ProvingError { + WrongParameter(String), +} + +pub trait IsStarkProver { + fn prove(trace: &TraceTable, pub_inputs: &A::PublicInputs) { + let air = A::new(trace.n_rows(), pub_inputs); + let domain = Domain::new(&air); + let lde_domain_length = domain.blowup_factor * domain.trace_length; + + // For each column, calculate the coefficients of the trace interpolating polynomial. + let mut trace_polys = trace.compute_trace_polys(); + + // Evaluate each polynomial in the lde domain. + let lde_trace_evaluations = trace_polys + .iter_mut() + .map(|poly| { + // Padding with zeros the coefficients of the polynomial, so we can evaluate it in the lde domain. + poly.resize(lde_domain_length, FieldElement::zero()); + evaluate_cfft(poly.to_vec()) + }) + .collect::>>>(); + + // TODO: Commit on lde trace evaluations. + + // --------- VALIDATE LDE TRACE EVALUATIONS ------------ + + // Interpolate lde trace evaluations. + let lde_coefficients = lde_trace_evaluations + .iter() + .map(|evals| interpolate_cfft(evals.to_vec())) + .collect::>>>(); + + // Evaluate lde trace interpolating polynomial in trace domain. + for point in &domain.trace_coset_points { + // println!("{:?}", evaluate_point(&lde_coefficients[0], &point)); + } + + // Crate a LDE_TRACE with a blow up factor of one, so its the same values as the trace. + let not_lde_trace = LDETraceTable::new(trace.table.data.clone(), 1, 1); + + // --------- VALIDATE BOUNDARY CONSTRAINTS ------------ + air.boundary_constraints() + .constraints + .iter() + .for_each(|constraint| { + let col = constraint.col; + let step = constraint.step; + let boundary_value = constraint.value.clone(); + + let trace_value = trace.table.get(step, col).clone(); + + if boundary_value.clone() != trace_value { + // println!("Boundary constraint inconsistency - Expected value {:?} in step {} and column {}, found: {:?}", boundary_value, step, col, trace_value); + } else { + // println!("Consistent Boundary constraint - Expected value {:?} in step {} and column {}, found: {:?}", boundary_value, step, col, trace_value) + } + }); + + // --------- VALIDATE TRANSITION CONSTRAINTS ----------- + for row_index in 0..not_lde_trace.table.height - 2 { + let frame = Frame::read_from_lde(¬_lde_trace, row_index, &air.context().transition_offsets); + let evaluations = air.compute_transition_prover(&frame); + // println!("Transition constraints evaluations: {:?}", evaluations); + } + + + let transition_coefficients: Vec> = vec![FieldElement::::one(); air.context().num_transition_constraints()]; + let boundary_coefficients: Vec> = vec![FieldElement::::one(); air.boundary_constraints().constraints.len() / 2]; + + // Compute the evaluations of the composition polynomial on the LDE domain. + let lde_trace = LDETraceTable::from_columns(lde_trace_evaluations, domain.blowup_factor); + let evaluator = ConstraintEvaluator::new(&air); + let constraint_evaluations = evaluator.evaluate( + &air, + &lde_trace, + &domain, + &transition_coefficients, + &boundary_coefficients, + ); + } +} + +// const BLOW_UP_FACTOR: usize = 2; + +// pub fn prove(trace: Vec>) -> Commitment { + +// let lde_domain_size = trace.len() * BLOW_UP_FACTOR; + +// // Returns the coef of the interpolating polinomial of the trace on a natural domain. +// let mut trace_poly = interpolate_cfft(trace); + +// // Padding with zeros the coefficients of the polynomial, so we can evaluate it in the lde domain. +// trace_poly.resize(lde_domain_size, FieldElement::zero()); +// let lde_eval = evaluate_cfft(trace_poly); + +// let tree = FriMerkleTree::::build(&lde_eval).unwrap(); +// let commitment = tree.root; + +// commitment +// } + +#[cfg(test)] +mod tests { + + use super::*; + + type FE = FieldElement; + + // #[test] + // fn basic_test() { + // let trace = vec![ + // FE::from(1), + // FE::from(2), + // FE::from(3), + // FE::from(4), + // FE::from(5), + // FE::from(6), + // FE::from(7), + // FE::from(8), + // ]; + + // let commitmet = prove(trace); + // println!("{:?}", commitmet); + // } +} diff --git a/provers/circle_stark/src/table.rs b/provers/circle_stark/src/table.rs new file mode 100644 index 000000000..18133d112 --- /dev/null +++ b/provers/circle_stark/src/table.rs @@ -0,0 +1,111 @@ +use lambdaworks_math::field::{ + element::FieldElement, fields::mersenne31::field::Mersenne31Field, +}; + +/// A two-dimensional Table holding field elements, arranged in a row-major order. +/// This is the basic underlying data structure used for any two-dimensional component in the +/// the STARK protocol implementation, such as the `TraceTable` and the `EvaluationFrame`. +/// Since this struct is a representation of a two-dimensional table, all rows should have the same +/// length. +#[derive(Clone, Default, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)] +pub struct Table { + pub data: Vec>, + pub width: usize, + pub height: usize, +} + +impl Table { + /// Crates a new Table instance from a one-dimensional array in row major order + /// and the intended width of the table. + pub fn new(data: Vec>, width: usize) -> Self { + // Check if the intented width is 0, used for creating an empty table. + if width == 0 { + return Self { + data: Vec::new(), + width, + height: 0, + }; + } + + // Check that the one-dimensional data makes sense to be interpreted as a 2D one. + // debug_assert!(crate::debug::validate_2d_structure(&data, width)); + let height = data.len() / width; + + Self { + data, + width, + height, + } + } + + /// Creates a Table instance from a vector of the intended columns. + pub fn from_columns(columns: Vec>>) -> Self { + if columns.is_empty() { + return Self::new(Vec::new(), 0); + } + let height = columns[0].len(); + + // Check that all columns have the same length for integrity + // debug_assert!(columns.iter().all(|c| c.len() == height)); + + let width = columns.len(); + let mut data = Vec::with_capacity(width * height); + + for row_idx in 0..height { + for column in columns.iter() { + data.push(column[row_idx].clone()); + } + } + + Self::new(data, width) + } + + /// Returns a vector of vectors of field elements representing the table rows + pub fn rows(&self) -> Vec>> { + self.data.chunks(self.width).map(|r| r.to_vec()).collect() + } + + /// Given a row index, returns a reference to that row as a slice of field elements. + pub fn get_row(&self, row_idx: usize) -> &[FieldElement] { + let row_offset = row_idx * self.width; + &self.data[row_offset..row_offset + self.width] + } + + /// Given a row index, returns a mutable reference to that row as a slice of field elements. + pub fn get_row_mut(&mut self, row_idx: usize) -> &mut [FieldElement] { + let n_cols = self.width; + let row_offset = row_idx * n_cols; + &mut self.data[row_offset..row_offset + n_cols] + } + + /// Given a slice of field elements representing a row, appends it to + /// the end of the table. + pub fn append_row(&mut self, row: &[FieldElement]) { + debug_assert_eq!(row.len(), self.width); + self.data.extend_from_slice(row); + self.height += 1 + } + + /// Returns a reference to the last row of the table + pub fn last_row(&self) -> &[FieldElement] { + self.get_row(self.height - 1) + } + + /// Returns a vector of vectors of field elements representing the table + /// columns + pub fn columns(&self) -> Vec>> { + (0..self.width) + .map(|col_idx| { + (0..self.height) + .map(|row_idx| self.data[row_idx * self.width + col_idx].clone()) + .collect() + }) + .collect() + } + + /// Given row and column indexes, returns the stored field element in that position of the table. + pub fn get(&self, row: usize, col: usize) -> &FieldElement { + let idx = row * self.width + col; + &self.data[idx] + } +} diff --git a/provers/circle_stark/src/tests/integration.rs b/provers/circle_stark/src/tests/integration.rs new file mode 100644 index 000000000..9fe78a89b --- /dev/null +++ b/provers/circle_stark/src/tests/integration.rs @@ -0,0 +1,34 @@ +use lambdaworks_math::field::{ + element::FieldElement, + fields::{ + fft_friendly::stark_252_prime_field::Stark252PrimeField, mersenne31::field::Mersenne31Field, + }, +}; + +use crate::{ + examples::simple_fibonacci::{self, FibonacciAIR, FibonacciPublicInputs}, + prover::{IsStarkProver, Prover}, + // proof::options::ProofOptions, + // prover::{IsStarkProver, Prover}, +}; + +#[test_log::test] +fn test_prove_fib() { + type FE = FieldElement; + + let trace = simple_fibonacci::fibonacci_trace([FE::one(), FE::one()], 32); + + let pub_inputs = FibonacciPublicInputs { + a0: FE::one(), + a1: FE::one(), + }; + + let proof = Prover::::prove(&trace, &pub_inputs); + // .unwrap(); + // assert!(Verifier::>::verify( + // &proof, + // &pub_inputs, + // &proof_options, + // StoneProverTranscript::new(&[]), + // )); +} diff --git a/provers/circle_stark/src/tests/mod.rs b/provers/circle_stark/src/tests/mod.rs new file mode 100644 index 000000000..605134bbc --- /dev/null +++ b/provers/circle_stark/src/tests/mod.rs @@ -0,0 +1 @@ +mod integration; \ No newline at end of file diff --git a/provers/circle_stark/src/trace.rs b/provers/circle_stark/src/trace.rs new file mode 100644 index 000000000..4ebb5535e --- /dev/null +++ b/provers/circle_stark/src/trace.rs @@ -0,0 +1,251 @@ +use crate::table::Table; +use itertools::Itertools; +use lambdaworks_math::{ + circle::{ + point::CirclePoint, + polynomial::{evaluate_point, interpolate_cfft}, + }, + fft::errors::FFTError, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; +#[cfg(feature = "parallel")] +use rayon::prelude::{IntoParallelRefIterator, ParallelIterator}; + +/// A two-dimensional representation of an execution trace of the STARK +/// protocol. +/// +/// For the moment it is mostly a wrapper around the `Table` struct. It is a +/// layer above the raw two-dimensional table, with functionality relevant to the +/// STARK protocol, such as the step size (number of consecutive rows of the table) +/// of the computation being proven. +#[derive(Clone, Default, Debug, PartialEq, Eq)] +pub struct TraceTable { + pub table: Table, + pub num_columns: usize, +} + +impl TraceTable { + pub fn new(data: Vec>, num_columns: usize) -> Self { + let table = Table::new(data, num_columns); + Self { table, num_columns } + } + + pub fn from_columns(columns: Vec>>) -> Self { + println!("COLUMNS LEN: {}", columns.len()); + let num_columns = columns.len(); + let table = Table::from_columns(columns); + Self { table, num_columns } + } + + pub fn empty() -> Self { + Self::new(Vec::new(), 0) + } + + pub fn is_empty(&self) -> bool { + self.table.width == 0 + } + + pub fn n_rows(&self) -> usize { + self.table.height + } + + pub fn n_cols(&self) -> usize { + self.table.width + } + + pub fn rows(&self) -> Vec>> { + self.table.rows() + } + + pub fn get_row(&self, row_idx: usize) -> &[FieldElement] { + self.table.get_row(row_idx) + } + + pub fn get_row_mut(&mut self, row_idx: usize) -> &mut [FieldElement] { + self.table.get_row_mut(row_idx) + } + + pub fn last_row(&self) -> &[FieldElement] { + self.get_row(self.n_rows() - 1) + } + + pub fn columns(&self) -> Vec>> { + self.table.columns() + } + + /// Given a slice of integer numbers representing column indexes, merge these columns into + /// a one-dimensional vector. + /// + /// The particular way they are merged is not really important since this function is used to + /// aggreagate values distributed across various columns with no importance on their ordering, + /// such as to sort them. + pub fn merge_columns(&self, column_indexes: &[usize]) -> Vec> { + let mut data = Vec::with_capacity(self.n_rows() * column_indexes.len()); + for row_index in 0..self.n_rows() { + for column in column_indexes { + data.push(self.table.data[row_index * self.n_cols() + column].clone()); + } + } + data + } + + pub fn compute_trace_polys(&self) -> Vec>> { + let columns = self.columns(); + #[cfg(feature = "parallel")] + let iter = columns.par_iter(); + #[cfg(not(feature = "parallel"))] + let iter = columns.iter(); + // FIX: Replace the .to_vec() + iter.map(|col| interpolate_cfft(col.to_vec())) + .collect::>>>() + } + + /// Given the padding length, appends the last row of the trace table + /// that many times. + /// This is useful for example when the desired trace length should be power + /// of two, and only the last row is the one that can be appended without affecting + /// the integrity of the constraints. + pub fn pad_with_last_row(&mut self, padding_len: usize) { + let last_row = self.last_row().to_vec(); + (0..padding_len).for_each(|_| { + self.table.append_row(&last_row); + }) + } + + /// Given a row index, a column index and a value, tries to set that location + /// of the trace with the given value. + /// The row_idx passed as argument may be greater than the max row index by 1. In this case, + /// last row of the trace is cloned, and the value is set in that cloned row. Then, the row is + /// appended to the end of the trace. + pub fn set_or_extend( + &mut self, + row_idx: usize, + col_idx: usize, + value: &FieldElement, + ) { + debug_assert!(col_idx < self.n_cols()); + // NOTE: This is not very nice, but for how the function is being used at the moment, + // the passed `row_idx` should never be greater than `self.n_rows() + 1`. This is just + // an integrity check for ease in the developing process, we should think a better alternative + // in the future. + debug_assert!(row_idx <= self.n_rows() + 1); + if row_idx >= self.n_rows() { + let mut last_row = self.last_row().to_vec(); + last_row[col_idx] = value.clone(); + self.table.append_row(&last_row) + } else { + let row = self.get_row_mut(row_idx); + row[col_idx] = value.clone(); + } + } +} +pub struct LDETraceTable { + pub(crate) table: Table, + pub(crate) blowup_factor: usize, +} + +impl LDETraceTable { + pub fn new( + data: Vec>, + n_columns: usize, + blowup_factor: usize, + ) -> Self { + let table = Table::new(data, n_columns); + + Self { + table, + blowup_factor, + } + } + + pub fn from_columns( + columns: Vec>>, + blowup_factor: usize, + ) -> Self { + let table = Table::from_columns(columns); + + Self { + table, + blowup_factor, + } + } + + pub fn num_cols(&self) -> usize { + self.table.width + } + + pub fn num_rows(&self) -> usize { + self.table.height + } + + pub fn get_row(&self, row_idx: usize) -> &[FieldElement] { + self.table.get_row(row_idx) + } + + pub fn get_table(&self, row: usize, col: usize) -> &FieldElement { + self.table.get(row, col) + } +} + +/// Given a slice of trace polynomials, an evaluation point `x`, the frame offsets +/// corresponding to the computation of the transitions, and a primitive root, +/// outputs the trace evaluations of each trace polynomial over the values used to +/// compute a transition. +/// Example: For a simple Fibonacci computation, if t(x) is the trace polynomial of +/// the computation, this will output evaluations t(x), t(g * x), t(g^2 * z). +pub fn get_trace_evaluations( + trace_polys: &[Vec>], + point: &CirclePoint, + frame_offsets: &[usize], + group_generator: &CirclePoint, +) -> Table { + let evaluation_points = frame_offsets + .iter() + .map(|offset| (group_generator * (*offset as u128)) + point) + .collect_vec(); + + let evaluations: Vec<_> = evaluation_points + .iter() + .flat_map(|eval_point| { + trace_polys + .iter() + .map(|poly| evaluate_point(poly, eval_point)) + }) + .collect(); + + let table_width = trace_polys.len(); + Table::new(evaluations, table_width) +} + +pub fn columns2rows( + columns: Vec>>, +) -> Vec>> { + let num_rows = columns[0].len(); + let num_cols = columns.len(); + + (0..num_rows) + .map(|row_index| { + (0..num_cols) + .map(|col_index| columns[col_index][row_index].clone()) + .collect() + }) + .collect() +} + +// #[cfg(test)] +// mod test { +// use super::TraceTable; +// use lambdaworks_math::field::{element::FieldElement, fields::u64_prime_field::F17}; +// type FE = FieldElement; + +// #[test] +// fn test_cols() { +// let col_1 = vec![FE::from(1), FE::from(2), FE::from(5), FE::from(13)]; +// let col_2 = vec![FE::from(1), FE::from(3), FE::from(8), FE::from(21)]; + +// let trace_table = TraceTable::from_columns(vec![col_1.clone(), col_2.clone()]); +// let res_cols = trace_table.columns(); + +// assert_eq!(res_cols, vec![col_1, col_2]); +// } +// } diff --git a/provers/circle_stark/src/vanishing_poly.rs b/provers/circle_stark/src/vanishing_poly.rs new file mode 100644 index 000000000..302f8d003 --- /dev/null +++ b/provers/circle_stark/src/vanishing_poly.rs @@ -0,0 +1,71 @@ +use lambdaworks_math::{ + circle::point::CirclePoint, + field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field}, +}; + +/// Evaluate the vanishing polynomial of the standard coset of size 2^log_2_size in a point. +/// The vanishing polynomial of a coset is the polynomial that takes the value zero when evaluated +/// in all the points of the coset. +/// We are using that if we take a point in g_{2n} + and double it n-1 times, then +/// we'll get the point (0, 1) or (0, -1); so its coordinate x is always 0. +pub fn evaluate_vanishing_poly( + log_2_size: u32, + point: &CirclePoint, +) -> FieldElement { + let mut x = point.x; + for _ in 1..log_2_size { + x = x.square().double() - FieldElement::one(); + } + x +} + +// Evaluate the polynomial that vanishes at a specific point in the domain at an arbitrary point. +// This use the "tangent" line to the domain in the vanish point. +// Check: https://vitalik.eth.limo/general/2024/07/23/circlestarks.html for details. +pub fn evaluate_single_point_zerofier( + vanish_point: CirclePoint, + eval_point: &CirclePoint, +) -> FieldElement { + (eval_point + vanish_point.conjugate()).x - FieldElement::::one() +} + +#[cfg(test)] +mod tests { + use super::*; + use lambdaworks_math::circle::cosets::Coset; + + type FE = FieldElement; + + #[test] + fn vanishing_poly_vanishes_in_coset() { + let log_2_size = 3; + let coset = Coset::new_standard(log_2_size); + let points = Coset::get_coset_points(&coset); + for point in points { + assert_eq!(evaluate_vanishing_poly(log_2_size, &point), FE::zero()); + } + } + #[test] + fn vanishing_poly_doesnt_vanishes_outside_coset() { + let log_2_size = 3; + let coset = Coset::new_standard(log_2_size + 1); + let points = Coset::get_coset_points(&coset); + for point in points { + assert_ne!(evaluate_vanishing_poly(log_2_size, &point), FE::zero()); + } + } + + #[test] + fn single_point_zerofier_vanishes_only_in_vanish_point() { + let vanish_point = CirclePoint::GENERATOR; + let eval_point = &vanish_point * 3; + assert_eq!( + evaluate_single_point_zerofier(vanish_point.clone(), &vanish_point), + FE::zero() + ); + assert_ne!( + evaluate_single_point_zerofier(vanish_point.clone(), &eval_point), + FE::zero() + ); + } +}