diff --git a/src/error.rs b/src/error.rs index de7091eb..37087058 100644 --- a/src/error.rs +++ b/src/error.rs @@ -93,6 +93,9 @@ pub enum Error { /// Index of the offending polynomial. label: String, }, + + /// Attempt to `open_amortized` on too large a domain + AmortizedOpeningTooLarge(usize), } impl core::fmt::Display for Error { @@ -179,6 +182,8 @@ impl core::fmt::Display for Error { support up to degree ({:?})", label, poly_degree, supported_degree ), Error::IncorrectInputLength(err) => write!(f, "{}", err), + Error::AmortizedOpeningTooLarge(s) => write!(f, "tried to open_amortized on too large domain of size {:?}", s), + } } } diff --git a/src/kzg10/data_structures.rs b/src/kzg10/data_structures.rs index d14cbc23..405ef689 100644 --- a/src/kzg10/data_structures.rs +++ b/src/kzg10/data_structures.rs @@ -26,6 +26,8 @@ pub struct UniversalParams { pub h: E::G2Affine, /// \beta times the above generator of G2. pub beta_h: E::G2Affine, + /// Group elements of the form `{ \beta^i H }`, where `i` ranges from 0 to `degree`. + pub powers_of_h: Vec, /// Group elements of the form `{ \beta^i G2 }`, where `i` ranges from `0` to `-degree`. pub neg_powers_of_h: BTreeMap, /// The generator of G2, prepared for use in pairings. @@ -100,6 +102,7 @@ impl CanonicalDeserialize for UniversalParams { powers_of_gamma_g, h, beta_h, + powers_of_h: vec![], neg_powers_of_h, prepared_h, prepared_beta_h, @@ -123,6 +126,7 @@ impl CanonicalDeserialize for UniversalParams { powers_of_gamma_g, h, beta_h, + powers_of_h: vec![], neg_powers_of_h, prepared_h, prepared_beta_h, @@ -144,6 +148,7 @@ impl CanonicalDeserialize for UniversalParams { powers_of_gamma_g, h, beta_h, + powers_of_h: vec![], neg_powers_of_h, prepared_h, prepared_beta_h, @@ -618,3 +623,36 @@ impl ToBytes for Proof { .write(&mut writer) } } + +/// Opening proofs of a commitment on a large domain +pub struct DomainProof { + /// This is a vector of commitments to the witness polynomials + /// over a domain 1, omega, omega^2, ..., omega^{n-1} + /// where omega is a primitive n'th root of unity + pub w: Vec, + /// Scale factor whose multiplication is deferred until the proofs are combined + pub scale: E::Fr, +} + +impl DomainProof { + /// Combine opening proofs onto a subset of the domain + /// represented by the SubproductDomain s + pub fn combine_at_domain( + &self, + start: usize, // Domain is omega^{start}, ..., omega^{end-1} + end: usize, + s: &super::subproductdomain::SubproductDomain, // SubproductDomain of the domain + ) -> Proof { + let lagrange_coeff = s.inverse_lagrange_coefficients(); + let mut total = E::G1Projective::zero(); + for (c_i, point) in lagrange_coeff.iter().zip(self.w[start..end].iter()) { + total += point.into_affine().mul::(c_i.inverse().unwrap()); + } + Proof { + w: total.into(), + // NOTE: if the ifft had not multiplied by domain_size_inv + //w: total.into_affine().mul::(self.scale).into(), + random_v: None, + } + } +} diff --git a/src/kzg10/mod.rs b/src/kzg10/mod.rs index e4f1c49c..84f56666 100644 --- a/src/kzg10/mod.rs +++ b/src/kzg10/mod.rs @@ -9,7 +9,9 @@ use crate::{BTreeMap, Error, LabeledPolynomial, PCRandomness, ToString, Vec}; use ark_ec::msm::{FixedBaseMSM, VariableBaseMSM}; use ark_ec::{group::Group, AffineCurve, PairingEngine, ProjectiveCurve}; use ark_ff::{One, PrimeField, UniformRand, Zero}; -use ark_poly::UVPolynomial; +use ark_poly::polynomial::univariate::DensePolynomial; +use ark_poly::EvaluationDomain; +use ark_poly::{Polynomial, UVPolynomial}; use ark_std::{format, marker::PhantomData, ops::Div, vec}; use ark_std::rand::RngCore; @@ -18,6 +20,8 @@ use rayon::prelude::*; mod data_structures; pub use data_structures::*; +mod subproductdomain; +pub use subproductdomain::*; /// `KZG10` is an implementation of the polynomial commitment scheme of /// [Kate, Zaverucha and Goldbgerg][kzg10] @@ -39,6 +43,7 @@ where pub fn setup( max_degree: usize, produce_g2_powers: bool, + produce_h_powers: bool, rng: &mut R, ) -> Result, Error> { if max_degree < 1 { @@ -119,16 +124,37 @@ where end_timer!(neg_powers_of_h_time); + let powers_of_h_time = start_timer!(|| "Generating powers of h in G2"); + let powers_of_h = if produce_h_powers { + let h_table = FixedBaseMSM::get_window_table(scalar_bits, window_size, h); + + let powers_of_h = FixedBaseMSM::multi_scalar_mul::( + scalar_bits, + window_size, + &h_table, + &powers_of_beta, + ); + + E::G2Projective::batch_normalization_into_affine(&powers_of_h) + } else { + let h = h.into_affine(); + let beta_h = h.mul(beta).into_affine(); + vec![h, beta_h] + }; + end_timer!(powers_of_h_time); + let h = h.into_affine(); let beta_h = h.mul(beta).into_affine(); + let prepared_h = h.into(); - let prepared_beta_h = beta_h.into(); + let prepared_beta_h = powers_of_h[1].into(); let pp = UniversalParams { powers_of_g, powers_of_gamma_g, h, beta_h, + powers_of_h, neg_powers_of_h, prepared_h, prepared_beta_h, @@ -138,7 +164,7 @@ where } /// Outputs a commitment to `polynomial`. - pub fn commit( + pub fn commit_g1( powers: &Powers, polynomial: &P, hiding_bound: Option, @@ -191,6 +217,30 @@ where Ok((Commitment(commitment.into()), randomness)) } + /// Outputs a commitment to `polynomial` in G2 + pub fn commit_g2( + powers: &UniversalParams, + polynomial: &DensePolynomial, + ) -> Result { + Self::check_degree_is_too_large(polynomial.degree(), powers.powers_of_h.len())?; + + let commit_time = start_timer!(|| format!( + "Committing to polynomial of degree {} in G2", + polynomial.degree(), + )); + + let (num_leading_zeros, plain_coeffs) = + skip_leading_zeros_and_convert_to_bigints(polynomial); + + let msm_time = start_timer!(|| "MSM to compute commitment to plaintext poly in G2"); + let commitment = VariableBaseMSM::multi_scalar_mul( + &powers.powers_of_h[num_leading_zeros..], + &plain_coeffs, + ); + end_timer!(msm_time); + end_timer!(commit_time); + Ok(commitment) + } /// Compute witness polynomial. /// /// The witness polynomial w(x) the quotient of the division (p(x) - p(z)) / (x - z) @@ -265,7 +315,7 @@ where } /// On input a polynomial `p` and a point `point`, outputs a proof for the same. - pub(crate) fn open<'a>( + pub fn open<'a>( powers: &Powers, p: &P, point: P::Point, @@ -290,6 +340,32 @@ where proof } + /// Open a polynomial at an entire evaluation domain simultaneously + pub fn open_at_domain<'a>( + pp: &UniversalParams, + p: &DensePolynomial, + domain: &ark_poly::domain::Radix2EvaluationDomain, + ) -> Result, Error> { + Self::check_degree_is_too_large(p.degree(), pp.powers_of_g.len())?; + let open_time = + start_timer!(|| format!("Amortized opening polynomial of degree {}", p.degree())); + + let m = p.coeffs().len() - 1; + let powers_of_g = pp.powers_of_g[0..m].to_vec(); + + let toeplitz_time = start_timer!(|| "Computing Toeplitz product"); + let (mut h, scale) = toeplitz_mul::(p, &powers_of_g, domain.size())?; + + end_timer!(toeplitz_time); + + domain.fft_in_place(&mut h); + + let proofs = DomainProof { w: h, scale }; + + end_timer!(open_time); + Ok(proofs) + } + /// Verifies that `value` is the evaluation at `point` of the polynomial /// committed inside `comm`. pub fn check( @@ -427,6 +503,60 @@ where Ok(()) } } + + /// Verifies that `evals` are the evaluation at `s` of the polynomial + /// committed inside `comm` where `s` is a constructed SubproductDomain + pub fn check_at_domain( + pp: &UniversalParams, + ck: &Powers, + comm: &Commitment, + proof: &Proof, // A multi-reveal KZG proof over the SubproductDomain s + evals: &[E::Fr], // Evaluations at the SubproductDomain + s: &SubproductDomain, // SubproductDomain of the evaluation domain + ) -> Result { + let evaluation_interpolation_time = start_timer!(|| "Constructing evaluation polynomial"); + + let evaluation_polynomial = s.interpolate(evals); + end_timer!(evaluation_interpolation_time); + let evaluation_commit_time = + start_timer!(|| "Constructing commitment to evaluation polynomial in G1"); + + let evaluation_polynomial_commitment = + KZG10::>::commit_g1(&ck, &evaluation_polynomial, None, None)?; + + end_timer!(evaluation_commit_time); + + let s_commitment_time = + start_timer!(|| "Constructing commitment to SubproductDomain in G2"); + let s_commitment = Self::commit_g2(pp, &s.t.m)?; + end_timer!(s_commitment_time); + Self::check_at_domain_with_commitments( + pp, + comm, + proof, + &s_commitment, + &evaluation_polynomial_commitment.0, + ) + } + + /// Check combined domain proof with precomputed commitments + pub fn check_at_domain_with_commitments( + ck: &UniversalParams, + comm: &Commitment, + proof: &Proof, + s_commitment: &E::G2Projective, // G2 Commitment to the SubproductDomain + evaluation_polynomial_commitment: &Commitment, + ) -> Result { + let check_time = start_timer!(|| "Checking evaluation"); + let inner = comm.0.into_projective() - evaluation_polynomial_commitment.0.into_projective(); + + let lhs = E::pairing(inner, ck.powers_of_h[0]); + + let rhs = E::pairing(proof.w, *s_commitment); + + end_timer!(check_time, || format!("Result: {}", lhs == rhs)); + Ok(lhs == rhs) + } } fn skip_leading_zeros_and_convert_to_bigints>( @@ -512,12 +642,12 @@ mod tests { f_p += (f, &p); let degree = 4; - let pp = KZG_Bls12_381::setup(degree, false, rng).unwrap(); + let pp = KZG_Bls12_381::setup(degree, false, false, rng).unwrap(); let (powers, _) = KZG_Bls12_381::trim(&pp, degree).unwrap(); let hiding_bound = None; - let (comm, _) = KZG10::commit(&powers, &p, hiding_bound, Some(rng)).unwrap(); - let (f_comm, _) = KZG10::commit(&powers, &f_p, hiding_bound, Some(rng)).unwrap(); + let (comm, _) = KZG10::commit_g1(&powers, &p, hiding_bound, Some(rng)).unwrap(); + let (f_comm, _) = KZG10::commit_g1(&powers, &f_p, hiding_bound, Some(rng)).unwrap(); let mut f_comm_2 = Commitment::empty(); f_comm_2 += (f, &comm); @@ -536,11 +666,11 @@ mod tests { while degree <= 1 { degree = usize::rand(rng) % 20; } - let pp = KZG10::::setup(degree, false, rng)?; + let pp = KZG10::::setup(degree, false, false, rng)?; let (ck, vk) = KZG10::::trim(&pp, degree)?; let p = P::rand(degree, rng); let hiding_bound = Some(1); - let (comm, rand) = KZG10::::commit(&ck, &p, hiding_bound, Some(rng))?; + let (comm, rand) = KZG10::::commit_g1(&ck, &p, hiding_bound, Some(rng))?; let point = E::Fr::rand(rng); let value = p.evaluate(&point); let proof = KZG10::::open(&ck, &p, point, &rand)?; @@ -564,11 +694,11 @@ mod tests { let rng = &mut test_rng(); for _ in 0..100 { let degree = 50; - let pp = KZG10::::setup(degree, false, rng)?; + let pp = KZG10::::setup(degree, false, false, rng)?; let (ck, vk) = KZG10::::trim(&pp, 2)?; let p = P::rand(1, rng); let hiding_bound = Some(1); - let (comm, rand) = KZG10::::commit(&ck, &p, hiding_bound, Some(rng))?; + let (comm, rand) = KZG10::::commit_g1(&ck, &p, hiding_bound, Some(rng))?; let point = E::Fr::rand(rng); let value = p.evaluate(&point); let proof = KZG10::::open(&ck, &p, point, &rand)?; @@ -595,7 +725,7 @@ mod tests { while degree <= 1 { degree = usize::rand(rng) % 20; } - let pp = KZG10::::setup(degree, false, rng)?; + let pp = KZG10::::setup(degree, false, false, rng)?; let (ck, vk) = KZG10::::trim(&pp, degree)?; let mut comms = Vec::new(); let mut values = Vec::new(); @@ -604,7 +734,7 @@ mod tests { for _ in 0..10 { let p = P::rand(degree, rng); let hiding_bound = Some(1); - let (comm, rand) = KZG10::::commit(&ck, &p, hiding_bound, Some(rng))?; + let (comm, rand) = KZG10::::commit_g1(&ck, &p, hiding_bound, Some(rng))?; let point = E::Fr::rand(rng); let value = p.evaluate(&point); let proof = KZG10::::open(&ck, &p, point, &rand)?; @@ -646,11 +776,132 @@ mod tests { let rng = &mut test_rng(); let max_degree = 123; - let pp = KZG_Bls12_381::setup(max_degree, false, rng).unwrap(); + let pp = KZG_Bls12_381::setup(max_degree, false, false, rng).unwrap(); let (powers, _) = KZG_Bls12_381::trim(&pp, max_degree).unwrap(); let p = DensePoly::::rand(max_degree + 1, rng); assert!(p.degree() > max_degree); assert!(KZG_Bls12_381::check_degree_is_too_large(p.degree(), powers.size()).is_err()); } + #[test] + fn test_toeplitz() { + type G1Projective = ::G1Projective; + type G1Affine = ::G1Affine; + + let rng = &mut test_rng(); + for total_weight in 15..30 { + let degree = 2 * total_weight / 3; + + let pp = KZG_Bls12_381::setup(degree, false, false, rng).unwrap(); + let (ck, _) = KZG_Bls12_381::trim(&pp, degree).unwrap(); + + for _ in 0..4 { + let polynomial = DensePolynomial::::rand(degree, rng); + + let coeffs = polynomial.coeffs(); + let m = polynomial.coeffs.len() - 1; + + let domain = ark_poly::Radix2EvaluationDomain::::new(total_weight).unwrap(); + let p = ck.powers_of_g[0..m].to_vec(); + let (h, _scale): (Vec, Fr) = + toeplitz_mul::(&polynomial, &p, domain.size()).unwrap(); + + let h = h + .iter() + .map( + |p: &G1Projective| (*p).into_affine(), /*. mul(_scale).into_affine()*/ + ) + .collect::>(); + + for i in 1..=m { + let mut total = G1Projective::zero(); + for j in 0..=m - i { + total += p[j].mul(coeffs[i + j]); + } + assert_eq!(G1Affine::from(total), h[i - 1]); + } + } + } + } + #[test] + fn test_check_at_domain() { + let rng = &mut test_rng(); + let n = 10; + let eval_domain_size = 1 << n; + + let degree = 1 << (n - 1); + + let pp = KZG_Bls12_381::setup(eval_domain_size, false, true, rng).unwrap(); + let (ck, vk) = KZG_Bls12_381::trim(&pp, degree).unwrap(); + + let p = DensePoly::::rand(degree, rng); + + let (comm, _) = KZG_Bls12_381::commit_g1(&ck, &p, None, None).unwrap(); + + let domain = ark_poly::Radix2EvaluationDomain::::new(eval_domain_size).unwrap(); + let openings = KZG_Bls12_381::open_at_domain(&pp, &p, &domain).unwrap(); + + let evals = p.evaluate_over_domain_by_ref(domain).evals; + + assert!(KZG_Bls12_381::check( + &vk, + &comm, + Fr::one(), + evals[0], + &Proof { + w: openings.w[0].into_affine(), + /*.mul::(openings.scale) // If the ifft does not multiply by domain_size_inv + .into_affine()*/ + random_v: None, + } + ) + .unwrap()); + + let mut whole_domain = Vec::with_capacity(domain.size()); + let mut t = Fr::one(); + for _ in 0..domain.size() { + whole_domain.push(t); + t *= domain.group_gen; + } + + let subproduct_domains = [ + SubproductDomain::new(whole_domain[0..100].to_vec()), + SubproductDomain::new(whole_domain[100..400].to_vec()), + SubproductDomain::new(whole_domain[400..].to_vec()), + ]; + + let combined = [ + openings.combine_at_domain(0, 100, &subproduct_domains[0]), + openings.combine_at_domain(100, 400, &subproduct_domains[1]), + openings.combine_at_domain(400, eval_domain_size, &subproduct_domains[2]), + ]; + + assert!(KZG_Bls12_381::check_at_domain( + &pp, + &ck, + &comm, + &combined[0], + &evals[0..100], + &subproduct_domains[0] + ) + .unwrap()); + assert!(KZG_Bls12_381::check_at_domain( + &pp, + &ck, + &comm, + &combined[1], + &evals[100..400], + &subproduct_domains[1] + ) + .unwrap()); + assert!(KZG_Bls12_381::check_at_domain( + &pp, + &ck, + &comm, + &combined[2], + &evals[400..], + &subproduct_domains[2] + ) + .unwrap()); + } } diff --git a/src/kzg10/subproductdomain.rs b/src/kzg10/subproductdomain.rs new file mode 100644 index 00000000..95761da7 --- /dev/null +++ b/src/kzg10/subproductdomain.rs @@ -0,0 +1,422 @@ +use crate::Error; +use ark_ec::PairingEngine; +use ark_ff::{FftField, Field, Zero}; +use ark_poly::{Polynomial, UVPolynomial}; + +use ark_poly::polynomial::univariate::DensePolynomial as Poly; + +//! Where indicated, algorithms are from Modern Computer Algebra, 3rd edition, by Gathen and Gerhard +//! Abbreviated as GG +//! Let M(n) denote the time to multiply. + +/// GG Algorithm 9.3 +/// Computes the inverse of f mod x^l +/// Takes O(M(l)) field arithmetic operations +pub fn inverse_mod_xl(f: &Poly, l: usize) -> Option> { + let r = ark_std::log2(l); // Compute ceil(log_2(l)) + if let Some(f_0_inv) = f.coeffs[0].inverse() { + // Invert f(0)^-1 if possible + let mut g = Poly:: { + coeffs: vec![f_0_inv], // Constant polynomial f(0)^-1 + }; + + let mut i = 2usize; + // Use Newton iteration which converges to the inverse mod x^l + for _ in 0..r { + g = &(&g + &g) - &(f * &(&g * &g)); //TODO: is g*2 better than g+g? + g.coeffs + .resize(ark_std::cmp::min(g.coeffs.len(), i), F::zero()); // Take g remainder mod x^{2^i} + i *= 2; + } + Some(g) + } else { + // No inverse exists because f(0)^-1 does not exist + None + } +} + +/// GG chapter 9.1 +/// Computes the rev_m(f) function in place +/// rev_m(f) = x^m f(1/x) +pub fn rev(f: &mut Poly, m: usize) { + assert!(f.coeffs.len() - 1 <= m); + for _ in 0..(m - (f.coeffs.len() - 1)) { + f.coeffs.push(F::zero()); + } + f.reverse(); +} + +/// GG Algorithm 9.5 +/// Divide f by g in nearly linear time +pub fn fast_divide_monic(f: &Poly, g: &Poly) -> (Poly, Poly) { + assert_eq!(g.coeffs.last(), F::one()); // check that `g` is monic + + if f.coeffs().len() < g.coeffs().len() { + return ( + Poly::zero(), + f.clone(), + ); + } + let m = f.coeffs().len() - g.coeffs().len(); + + let mut rev_f = f.clone(); + let mut rev_g = g.clone(); + rev_f.reverse(); + rev_g.reverse(); + + let mut q = &rev_f * &inverse_mod_xl::(&rev_g, m + 1).unwrap(); + q.coeffs.resize(m + 1, F::zero()); + rev::(&mut q, m); + let r = f - &(g * &q); + (q, r) +} + +/// A subproduct domain is a domain { u_0, ..., u_{n-1} } of scalar values +/// accompanied by a subproduct tree of the polynomial: +/// m = (x - u_0)*...*(x-u_{n-1}) +/// Once the subproduct tree is constructed, operations over the +/// entire subproduct domain can be done in a fast, recursive way. +/// Unlike other fast algorithms, the subproduct domain may be +/// arbitrary points instead of a multiplicative subgroup +/// of roots of unity of the field +#[derive(Debug, Clone)] +pub struct SubproductDomain { + /// Domain values u = { u_0, ..., u_{n-1} } + pub u: Vec, + /// Subproduct tree over domain u + pub t: SubproductTree, + /// Derivative of the subproduct polynomial + pub prime: Poly, // Derivative of the polynomial m +} + +impl SubproductDomain { + /// Create a new subproduct tree domain over the domain { u_0, ..., u_{n-1} } + pub fn new(u: Vec) -> SubproductDomain { + let t = SubproductTree::new(&u); + let prime = derivative::(&t.m); + SubproductDomain { u, t, prime } + } + /// Evaluate a polynomial f over the subproduct domain u + pub fn evaluate(&self, f: &Poly) -> Vec { + let mut evals = vec![F::zero(); self.u.len()]; + self.t.evaluate(f, &self.u, &mut evals); + evals + } + /// Interpolate a polynomial f over the domain, such that f(u_i) = v_i + pub fn interpolate(&self, v: &[F]) -> Poly { + self.t.interpolate(&self.u, v) + } + /// Compute the inverse of the lagrange coefficients necessary to interpolate over u + pub fn inverse_lagrange_coefficients(&self) -> Vec { + self.t.inverse_lagrange_coefficients(&self.u) + } + /// Compute a linear combination of lagrange factors times c_i + pub fn linear_combine(&self, c: &[F]) -> Poly { + self.t.linear_combine(&self.u, &c) + } +} + +/// A subproduct tree of the subproduct domain +/// This type is defined separately from SubproductDomain +/// because the domain u is owned by SubproductDomain, whereas +/// m = (x - u_i)*...*(x-u_j) is owned by the SubproductTree +/// The subdomain { u_i, ..., u_j } is borrowed by SubproductTree for each operation +#[derive(Debug, Clone)] +pub struct SubproductTree { + /// The left child SubproductTree + pub left: Option>>, + /// The right child SubproductTree + pub right: Option>>, + /// The polynomial m = (x - u_i)*...*(x-u_j) for this subdomain + pub m: Poly, +} + +impl SubproductTree { + /// GG Algorithm 10.3 + /// Compute the subproduct tree of m = (x - u_0)*...*(x-u_{n-1}) + /// Takes O(M(r) log r) field operations + /// Specialized to assume the leaves are of the form m_i = x-u_i + /// Generalized to arbitrary r, not just powers of 2 + pub fn new(u: &[F]) -> SubproductTree { + // A degree 1 polynomial is a leaf of the tree + if u.len() == 1 { + SubproductTree { + left: None, + right: None, + m: Poly:: { + coeffs: vec![-u[0], F::one()], // m_0 = x - u_0 + }, + } + } else { + // Split as evenly as possible + let n = u.len() / 2; + let (u_0, u_1) = u.split_at(n); + let left = Box::new(SubproductTree::new(u_0)); + let right = Box::new(SubproductTree::new(u_1)); + let m = &left.m * &right.m; + SubproductTree { + left: Some(left), + right: Some(right), + m, + } + } + } + /// GG algorithm 9.5 + /// Evaluate f over this subproduct tree + /// deg(f) must be less than the size of the domain + /// self must be the subproduct tree of the slice u + /// The output is stored in the slice t: + /// t_i = f(u_i) + /// Takes O(M(n) log n) field operations + pub fn evaluate(&self, f: &Poly, u: &[F], t: &mut [F]) { + assert!(f.degree() < u.len()); + + if u.len() == 1 { + // By the assertion above, f must be a constant polynomial, so evaluating + t[0] = f.coeffs[0]; + return; + } + + // if u.len() > 1, then this SubproductTree must not be a leaf, and it has both children + let left = self.left.as_ref().unwrap(); + let right = self.right.as_ref().unwrap(); + + // if f = q*m + r, then on the domain u where m(u_i) = 0, f(u_i) = r(u_i) + let (_, r_0) = fast_divide_monic::(f, &left.m); + let (_, r_1) = fast_divide_monic::(f, &right.m); + + // divide the domain in the same way that the SubproductTree was constructed + let n = u.len() / 2; + let (u_0, u_1) = u.split_at(n); + let (t_0, t_1) = t.split_at_mut(n); + + left.evaluate(&r_0, u_0, t_0); + right.evaluate(&r_1, u_1, t_1); + } + /// Fast interpolate over this subproduct tree + pub fn interpolate(&self, u: &[F], v: &[F]) -> Poly { + let mut lagrange_coeff = self.inverse_lagrange_coefficients(u); + + for (s_i, v_i) in lagrange_coeff.iter_mut().zip(v.iter()) { + *s_i = s_i.inverse().unwrap() * *v_i; + } + + self.linear_combine(u, &lagrange_coeff) + } + /// Fast compute lagrange coefficients over this subproduct tree + pub fn inverse_lagrange_coefficients(&self, u: &[F]) -> Vec { + //assert u.len() == degree of s.m + if u.len() == 1 { + return vec![F::one()]; + } + let mut evals = vec![F::zero(); u.len()]; + let m_prime = derivative::(&self.m); + self.evaluate(&m_prime, u, &mut evals); + evals + } + /// GG Algorithm 10.9 + /// Fast linear combination of moduli over this subproduct tree + /// On input c = { c_0, ..., c_{n-1} } + /// output sum_i (c_i * m) / (x- u_i) + /// Takes O(M(n) log n) field operations + pub fn linear_combine(&self, u: &[F], c: &[F]) -> Poly { + if u.len() == 1 { + // Output c_0 * (x-u_0) / (x-u_0) = c_0 + return Poly:: { coeffs: vec![c[0]] }; + } + let n = u.len() / 2; + let (u_0, u_1) = u.split_at(n); + let (c_0, c_1) = c.split_at(n); + + // Linear combination of moduli over both halves of the subproduct tree + let left = self.left.as_ref().unwrap(); + let right = self.right.as_ref().unwrap(); + let r_0 = left.linear_combine(u_0, c_0); + let r_1 = right.linear_combine(u_1, c_1); + + // r_0 = sum_{i in left} c_i m_left / (x-u_i) + // so m_right * r_0 = sum_{i in left} c_i (m_left*m_right) / (x-u_i) = sum_{i in left} c_i m / (x-u_i) + &(&right.m * &r_0) + &(&left.m * &r_1) + } +} +/// compute the derivative of polynomial f +pub fn derivative(f: &Poly) -> Poly { + let mut coeffs = Vec::with_capacity(f.coeffs().len() - 1); + for (i, c) in f.coeffs.iter().enumerate().skip(1) { + coeffs.push(F::from(i as u64) * c); + } + Poly:: { coeffs } +} + +/// Build a vector representation of the (n x n) circulant matrix of polynomial f +/// Based on the blog post: +/// https://alinush.github.io/2020/03/19/multiplying-a-vector-by-a-toeplitz-matrix.html +pub fn build_circulant(f: &Poly, size: usize) -> Vec { + let mut circulant = vec![F::zero(); 2 * size]; + let coeffs = f.coeffs(); + if size == coeffs.len() - 1 { + circulant[0] = *coeffs.last().unwrap(); + circulant[size] = *coeffs.last().unwrap(); + circulant[size + 1..size + 1 + coeffs.len() - 2] + .copy_from_slice(&coeffs[1..coeffs.len() - 1]); + } else { + circulant[size + 1..size + 1 + coeffs.len() - 1].copy_from_slice(&coeffs[1..]); + } + circulant +} +/// Computes the Toeplitz matrix of polynomial times the vector v +pub fn toeplitz_mul( + polynomial: &Poly, + v: &[E::G1Affine], + size: usize, +) -> Result<(Vec, E::Fr), Error> { + use ark_ec::AffineCurve; + use ark_poly::EvaluationDomain; + let m = polynomial.coeffs().len() - 1; + let size = ark_std::cmp::max(size, m); + + let domain = ark_poly::Radix2EvaluationDomain::::new(2 * size) + .ok_or(Error::AmortizedOpeningTooLarge(size))?; + + let size = domain.size() / 2; + let mut circulant = build_circulant::(polynomial, size); + + let mut tmp: Vec = Vec::with_capacity(domain.size()); + + for _ in 0..(size - v.len()) { + tmp.push(E::G1Projective::zero()); + } + + for i in v.iter().rev() { + tmp.push(i.into_projective()); + } + + tmp.resize(domain.size(), E::G1Projective::zero()); + domain.fft_in_place(&mut tmp); + domain.fft_in_place(&mut circulant); + + for (i, j) in tmp.iter_mut().zip(circulant.iter()) { + *i *= *j; + } + + // NOTE: it is possible to avoid scaling by domain_size_inv, and defer this to later + domain.ifft_in_place(&mut tmp); + + Ok(( + tmp[..size].to_vec(), + E::Fr::from(domain.size() as u64).inverse().unwrap(), + )) +} + +#[cfg(test)] +mod tests { + use super::*; + use ark_ec::PairingEngine; + use ark_ff::{One, Zero}; + use ark_poly::polynomial::univariate::DensePolynomial; + use ark_poly::Polynomial; + use ark_poly::UVPolynomial; + use ark_std::UniformRand; + + type Fr = ::Fr; + #[test] + fn test_inverse() { + let rng = &mut ark_std::test_rng(); + + let l = 101; + for l in [1, 2, 3, 5, 19, 25, 101].iter() { + for degree in 0..*l { + for _ in 0..10 { + let p = DensePolynomial::::rand(degree, rng); + let p_inv = inverse_mod_xl::(&p, *l).unwrap(); + let mut t = &p * &p_inv; + t.coeffs.resize(*l, Fr::zero()); + assert_eq!(t.coeffs[0], Fr::one()); + for i in t.iter().skip(1) { + assert_eq!(*i, Fr::zero()); + } + } + } + } + } + + #[test] + fn test_divide() { + let rng = &mut ark_std::test_rng(); + + let degree = 100; + for g_deg in 1..100 { + let f = DensePolynomial::::rand(degree, rng); + let mut g = DensePolynomial::::rand(g_deg, rng); + *g.last_mut().unwrap() = Fr::one(); //monic + + let (q, r) = fast_divide_monic::(&f, &g); + + let t = &(&q * &g) + &r; + + for (i, j) in t.coeffs.iter().zip(f.coeffs.iter()) { + assert_eq!(*i, *j); + } + } + } + + #[test] + fn test_interpolate() { + let rng = &mut ark_std::test_rng(); + for d in 1..100 { + let mut points = vec![]; + let mut evals = vec![]; + for _ in 0..d { + points.push(Fr::rand(rng)); + evals.push(Fr::rand(rng)); + } + + let s = SubproductDomain::::new(points); + let p = s.interpolate(&evals); + + for (x, y) in s.u.iter().zip(evals.iter()) { + assert_eq!(p.evaluate(x), *y) + } + } + } + + #[test] + fn test_linear_combine() { + let rng = &mut ark_std::test_rng(); + for d in 1..100 { + let mut u = vec![]; + let mut c = vec![]; + for _ in 0..d { + u.push(Fr::rand(rng)); + c.push(Fr::rand(rng)); + } + let s = SubproductDomain::::new(u); + let f = s.linear_combine(&c); + + let r = Fr::rand(rng); + let m = s.t.m.evaluate(&r); + let mut total = Fr::zero(); + for (u_i, c_i) in s.u.iter().zip(c.iter()) { + total += m * *c_i / (r - u_i); + } + assert_eq!(f.evaluate(&r), total); + } + } + + #[test] + fn test_inv_lagrange() { + let rng = &mut ark_std::test_rng(); + for d in 1..100 { + let mut u = vec![]; + for _ in 0..d { + u.push(Fr::rand(rng)); + } + let s = SubproductDomain::::new(u); + let f = s.inverse_lagrange_coefficients(); + + for (i, j) in s.u.iter().zip(f.iter()) { + assert_eq!(s.prime.evaluate(i), *j); + } + } + } +} diff --git a/src/marlin/marlin_pc/mod.rs b/src/marlin/marlin_pc/mod.rs index 918d5c4c..3f3d9619 100644 --- a/src/marlin/marlin_pc/mod.rs +++ b/src/marlin/marlin_pc/mod.rs @@ -78,7 +78,7 @@ where _num_vars: Option, rng: &mut R, ) -> Result { - kzg10::KZG10::setup(max_degree, false, rng).map_err(Into::into) + kzg10::KZG10::setup(max_degree, false, false, rng).map_err(Into::into) } fn trim( @@ -219,13 +219,13 @@ where )); let (comm, rand) = - kzg10::KZG10::commit(&ck.powers(), polynomial, hiding_bound, Some(rng))?; + kzg10::KZG10::commit_g1(&ck.powers(), polynomial, hiding_bound, Some(rng))?; let (shifted_comm, shifted_rand) = if let Some(degree_bound) = degree_bound { let shifted_powers = ck .shifted_powers(degree_bound) .ok_or(Error::UnsupportedDegreeBound(degree_bound))?; let (shifted_comm, shifted_rand) = - kzg10::KZG10::commit(&shifted_powers, &polynomial, hiding_bound, Some(rng))?; + kzg10::KZG10::commit_g1(&shifted_powers, &polynomial, hiding_bound, Some(rng))?; (Some(shifted_comm), Some(shifted_rand)) } else { (None, None) diff --git a/src/sonic_pc/mod.rs b/src/sonic_pc/mod.rs index e02d635f..7e2890ba 100644 --- a/src/sonic_pc/mod.rs +++ b/src/sonic_pc/mod.rs @@ -157,7 +157,7 @@ where _: Option, rng: &mut R, ) -> Result { - kzg10::KZG10::::setup(max_degree, true, rng).map_err(Into::into) + kzg10::KZG10::::setup(max_degree, true, false, rng).map_err(Into::into) } fn trim( @@ -326,7 +326,8 @@ where ck.powers() }; - let (comm, rand) = kzg10::KZG10::commit(&powers, polynomial, hiding_bound, Some(rng))?; + let (comm, rand) = + kzg10::KZG10::commit_g1(&powers, polynomial, hiding_bound, Some(rng))?; labeled_comms.push(LabeledCommitment::new( label.to_string(),