From 6d03dd06f5806dc7dbd46dc95461a62cb273236e Mon Sep 17 00:00:00 2001 From: Daniel Lubarov Date: Wed, 12 May 2021 11:20:47 -0700 Subject: [PATCH 1/2] Finish merging in old_polynomial --- src/polynomial/division.rs | 169 ++++++--- src/polynomial/mod.rs | 1 - src/polynomial/old_polynomial.rs | 580 ------------------------------- src/polynomial/polynomial.rs | 290 +++++++++++++++- src/prover.rs | 3 +- 5 files changed, 399 insertions(+), 644 deletions(-) delete mode 100644 src/polynomial/old_polynomial.rs diff --git a/src/polynomial/division.rs b/src/polynomial/division.rs index 0b5055ef..c72b798c 100644 --- a/src/polynomial/division.rs +++ b/src/polynomial/division.rs @@ -3,73 +3,140 @@ use crate::field::field::Field; use crate::polynomial::polynomial::PolynomialCoeffs; use crate::util::log2_strict; -/// Takes a polynomial `a` in coefficient form, and divides it by `Z_H = X^n - 1`. -/// -/// This assumes `Z_H | a`, otherwise result is meaningless. -pub(crate) fn divide_by_z_h(mut a: PolynomialCoeffs, n: usize) -> PolynomialCoeffs { - // TODO: Is this special case needed? - if a.coeffs.iter().all(|p| *p == F::ZERO) { - return a; +impl PolynomialCoeffs { + /// Polynomial division. + /// Returns `(q, r)`, the quotient and remainder of the polynomial division of `a` by `b`. + pub(crate) fn div_rem(&self, b: &Self) -> (Self, Self) { + let (a_degree_plug_1, b_degree_plus_1) = (self.degree_plus_one(), b.degree_plus_one()); + if self.is_zero() { + (Self::zero(1), Self::empty()) + } else if b.is_zero() { + panic!("Division by zero polynomial"); + } else if a_degree_plug_1 < b_degree_plus_1 { + (Self::zero(1), self.clone()) + } else if b_degree_plus_1 == 1 { + (self * b.coeffs[0].inverse(), Self::empty()) + } else { + let rev_b = b.rev(); + let rev_b_inv = rev_b.inv_mod_xn(a_degree_plug_1 - b_degree_plus_1 + 1); + let rhs: Self = self.rev().coeffs[..=a_degree_plug_1 - b_degree_plus_1] + .to_vec() + .into(); + let rev_q: Self = (&rev_b_inv * &rhs).coeffs[..=a_degree_plug_1 - b_degree_plus_1] + .to_vec() + .into(); + let mut q = rev_q.rev(); + let mut qb = &q * b; + let mut r = self - &qb; + q.trim(); + r.trim(); + (q, r) + } } - let g = F::MULTIPLICATIVE_GROUP_GENERATOR; - let mut g_pow = F::ONE; - // Multiply the i-th coefficient of `a` by `g^i`. Then `new_a(w^j) = old_a(g.w^j)`. - a.coeffs.iter_mut().for_each(|x| { - *x *= g_pow; - g_pow *= g; - }); + /// Polynomial long division. + /// Returns `(q, r)`, the quotient and remainder of the polynomial division of `a` by `b`. + /// Generally slower that the equivalent function `Polynomial::polynomial_division`. + pub fn div_rem_long_division(&self, b: &Self) -> (Self, Self) { + let b = b.trimmed(); - let root = F::primitive_root_of_unity(log2_strict(a.len())); - // Equals to the evaluation of `a` on `{g.w^i}`. - let mut a_eval = fft(a); - // Compute the denominators `1/(g^n.w^(n*i) - 1)` using batch inversion. - let denominator_g = g.exp_usize(n); - let root_n = root.exp_usize(n); - let mut root_pow = F::ONE; - let denominators = (0..a_eval.len()) - .map(|i| { - if i != 0 { - root_pow = root_pow * root_n; + let (a_degree_plus_1, b_degree_plus_1) = (self.degree_plus_one(), b.degree_plus_one()); + if self.is_zero() { + (Self::zero(1), Self::empty()) + } else if b.is_zero() { + panic!("Division by zero polynomial"); + } else if a_degree_plus_1 < b_degree_plus_1 { + (Self::zero(1), self.clone()) + } else { + // Now we know that self.degree() >= divisor.degree(); + let mut quotient = Self::zero(a_degree_plus_1 - b_degree_plus_1 + 1); + let mut remainder = self.clone(); + // Can unwrap here because we know self is not zero. + let divisor_leading_inv = b.lead().inverse(); + while !remainder.is_zero() && remainder.degree_plus_one() >= b_degree_plus_1 { + let cur_q_coeff = remainder.lead() * divisor_leading_inv; + let cur_q_degree = remainder.degree_plus_one() - b_degree_plus_1; + quotient.coeffs[cur_q_degree] = cur_q_coeff; + + for (i, &div_coeff) in b.coeffs.iter().enumerate() { + remainder.coeffs[cur_q_degree + i] = + remainder.coeffs[cur_q_degree + i] - (cur_q_coeff * div_coeff); + } + remainder.trim(); } - denominator_g * root_pow - F::ONE - }) - .collect::>(); - let denominators_inv = F::batch_multiplicative_inverse(&denominators); - // Divide every element of `a_eval` by the corresponding denominator. - // Then, `a_eval` is the evaluation of `a/Z_H` on `{g.w^i}`. - a_eval - .values - .iter_mut() - .zip(denominators_inv.iter()) - .for_each(|(x, &d)| { - *x *= d; + (quotient, remainder) + } + } + + /// Takes a polynomial `a` in coefficient form, and divides it by `Z_H = X^n - 1`. + /// + /// This assumes `Z_H | a`, otherwise result is meaningless. + pub(crate) fn divide_by_z_h(&self, n: usize) -> PolynomialCoeffs { + let mut a = self.clone(); + + // TODO: Is this special case needed? + if a.coeffs.iter().all(|p| *p == F::ZERO) { + return a; + } + + let g = F::MULTIPLICATIVE_GROUP_GENERATOR; + let mut g_pow = F::ONE; + // Multiply the i-th coefficient of `a` by `g^i`. Then `new_a(w^j) = old_a(g.w^j)`. + a.coeffs.iter_mut().for_each(|x| { + *x *= g_pow; + g_pow *= g; }); - // `p` is the interpolating polynomial of `a_eval` on `{w^i}`. - let mut p = ifft(a_eval); - // We need to scale it by `g^(-i)` to get the interpolating polynomial of `a_eval` on `{g.w^i}`, - // a.k.a `a/Z_H`. - let g_inv = g.inverse(); - let mut g_inv_pow = F::ONE; - p.coeffs.iter_mut().for_each(|x| { - *x *= g_inv_pow; - g_inv_pow *= g_inv; - }); - p + + let root = F::primitive_root_of_unity(log2_strict(a.len())); + // Equals to the evaluation of `a` on `{g.w^i}`. + let mut a_eval = fft(a); + // Compute the denominators `1/(g^n.w^(n*i) - 1)` using batch inversion. + let denominator_g = g.exp_usize(n); + let root_n = root.exp_usize(n); + let mut root_pow = F::ONE; + let denominators = (0..a_eval.len()) + .map(|i| { + if i != 0 { + root_pow = root_pow * root_n; + } + denominator_g * root_pow - F::ONE + }) + .collect::>(); + let denominators_inv = F::batch_multiplicative_inverse(&denominators); + // Divide every element of `a_eval` by the corresponding denominator. + // Then, `a_eval` is the evaluation of `a/Z_H` on `{g.w^i}`. + a_eval + .values + .iter_mut() + .zip(denominators_inv.iter()) + .for_each(|(x, &d)| { + *x *= d; + }); + // `p` is the interpolating polynomial of `a_eval` on `{w^i}`. + let mut p = ifft(a_eval); + // We need to scale it by `g^(-i)` to get the interpolating polynomial of `a_eval` on `{g.w^i}`, + // a.k.a `a/Z_H`. + let g_inv = g.inverse(); + let mut g_inv_pow = F::ONE; + p.coeffs.iter_mut().for_each(|x| { + *x *= g_inv_pow; + g_inv_pow *= g_inv; + }); + p + } } #[cfg(test)] mod tests { use crate::field::crandall_field::CrandallField; use crate::field::field::Field; - use crate::polynomial::division::divide_by_z_h; use crate::polynomial::polynomial::PolynomialCoeffs; #[test] fn zero_div_z_h() { type F = CrandallField; let zero = PolynomialCoeffs::::zero(16); - let quotient = divide_by_z_h(zero.clone(), 4); + let quotient = zero.divide_by_z_h(4); assert_eq!(quotient, zero); } @@ -91,7 +158,7 @@ mod tests { let a = PolynomialCoeffs::new(vec![-six, -five, -four, -three, six, five, four, three]); let q = PolynomialCoeffs::new(vec![six, five, four, three, zero, zero, zero, zero]); - let computed_q = divide_by_z_h(a, 4); + let computed_q = a.divide_by_z_h(4); assert_eq!(computed_q, q); } } diff --git a/src/polynomial/mod.rs b/src/polynomial/mod.rs index 314ed8e6..4360d3d4 100644 --- a/src/polynomial/mod.rs +++ b/src/polynomial/mod.rs @@ -1,4 +1,3 @@ pub mod commitment; pub(crate) mod division; -mod old_polynomial; pub mod polynomial; diff --git a/src/polynomial/old_polynomial.rs b/src/polynomial/old_polynomial.rs deleted file mode 100644 index d011fff7..00000000 --- a/src/polynomial/old_polynomial.rs +++ /dev/null @@ -1,580 +0,0 @@ -#![allow(clippy::many_single_char_names)] -use crate::field::fft::{ - fft_precompute, fft_with_precomputation_power_of_2, ifft_with_precomputation_power_of_2, - FftPrecomputation, -}; -use crate::field::field::Field; -use crate::polynomial::polynomial::PolynomialCoeffs; -use crate::util::log2_ceil; -use std::cmp::Ordering; -use std::ops::{Index, IndexMut, RangeBounds}; -use std::slice::{Iter, IterMut, SliceIndex}; - -/// Polynomial struct holding a polynomial in coefficient form. -// TODO: Finish merging this with `PolynomialCoeffs`. -#[deprecated] -#[derive(Debug, Clone)] -pub struct Polynomial(Vec); - -impl PartialEq for Polynomial { - fn eq(&self, other: &Self) -> bool { - let max_terms = self.0.len().max(other.0.len()); - for i in 0..max_terms { - let self_i = self.0.get(i).cloned().unwrap_or(F::ZERO); - let other_i = other.0.get(i).cloned().unwrap_or(F::ZERO); - if self_i != other_i { - return false; - } - } - true - } -} - -impl Eq for Polynomial {} - -impl From> for Polynomial { - /// Takes a vector of coefficients and returns the corresponding polynomial. - fn from(coeffs: Vec) -> Self { - Self(coeffs) - } -} - -impl From> for Polynomial { - fn from(coeffs: PolynomialCoeffs) -> Self { - Self(coeffs.coeffs) - } -} -impl From> for PolynomialCoeffs { - fn from(poly: Polynomial) -> Self { - Self::new(poly.0) - } -} - -impl Index for Polynomial -where - F: Field, - I: SliceIndex<[F]>, -{ - type Output = I::Output; - - /// Indexing on the coefficients. - fn index(&self, index: I) -> &Self::Output { - &self.0[index] - } -} - -impl IndexMut for Polynomial -where - F: Field, - I: SliceIndex<[F]>, -{ - fn index_mut(&mut self, index: I) -> &mut >::Output { - &mut self.0[index] - } -} - -impl Polynomial { - /// Takes a slice of coefficients and returns the corresponding polynomial. - pub fn from_coeffs(coeffs: &[F]) -> Self { - Self(coeffs.to_vec()) - } - - /// Returns the coefficient vector. - pub fn coeffs(&self) -> &[F] { - &self.0 - } - - /// Empty polynomial; - pub fn empty() -> Self { - Self(Vec::new()) - } - - /// Zero polynomial with length `len`. - /// `len = 1` is the standard representation, but sometimes it's useful to set `len > 1` - /// to have polynomials with uniform length. - pub fn zero(len: usize) -> Self { - Self(vec![F::ZERO; len]) - } - - pub fn iter(&self) -> Iter { - self.0.iter() - } - - pub fn iter_mut(&mut self) -> IterMut { - self.0.iter_mut() - } - - pub fn is_zero(&self) -> bool { - self.0.iter().all(|x| x.is_zero()) - } - - /// Number of coefficients held by the polynomial. Is NOT equal to the degree in general. - pub fn len(&self) -> usize { - self.0.len() - } - - /// Degree of the polynomial. - /// Panics on zero polynomial. - pub fn degree(&self) -> usize { - (0usize..self.len()) - .rev() - .find(|&i| self[i].is_nonzero()) - .expect("Zero polynomial") - } - - /// Degree of the polynomial + 1. - fn degree_plus_one(&self) -> usize { - (0usize..self.len()) - .rev() - .find(|&i| self[i].is_nonzero()) - .map_or(0, |i| i + 1) - } - - pub fn is_empty(&self) -> bool { - self.0.is_empty() - } - - /// Removes the coefficients in the range `range`. - fn drain>(&mut self, range: R) { - self.0.drain(range); - } - - /// Evaluates the polynomial at a point `x`. - pub fn eval(&self, x: F) -> F { - self.iter().rev().fold(F::ZERO, |acc, &c| acc * x + c) - } - - /// Evaluates the polynomial at a point `x`, given the list of powers of `x`. - /// Assumes that `self.len() == x_pow.len()`. - pub fn eval_from_power(&self, x_pow: &[F]) -> F { - self.iter() - .zip(x_pow) - .fold(F::ZERO, |acc, (&c, &p)| acc + c * p) - } - - /// Evaluates the polynomial on subgroup of `F^*` with a given FFT precomputation. - pub(crate) fn eval_domain(&self, fft_precomputation: &FftPrecomputation) -> Vec { - let domain_size = fft_precomputation.size(); - if self.len() < domain_size { - // Need to pad the polynomial to have the same length as the domain. - fft_with_precomputation_power_of_2( - self.padded(domain_size).coeffs().to_vec().into(), - fft_precomputation, - ) - .values - } else { - fft_with_precomputation_power_of_2(self.coeffs().to_vec().into(), fft_precomputation) - .values - } - } - - /// Computes the interpolating polynomial of a list of `values` on a subgroup of `F^*`. - pub(crate) fn from_evaluations( - values: &[F], - fft_precomputation: &FftPrecomputation, - ) -> Self { - Self(ifft_with_precomputation_power_of_2(values.to_vec().into(), fft_precomputation).coeffs) - } - - /// Leading coefficient. - pub fn lead(&self) -> F { - self.iter() - .rev() - .find(|x| x.is_nonzero()) - .map_or(F::ZERO, |x| *x) - } - - /// Reverse the order of the coefficients, not taking into account the leading zero coefficients. - fn rev(&self) -> Self { - let d = self.degree(); - Self(self.0[..=d].iter().rev().copied().collect()) - } - - /// Negates the polynomial's coefficients. - pub fn neg(&self) -> Self { - Self(self.iter().map(|&x| -x).collect()) - } - - /// Multiply the polynomial's coefficients by a scalar. - pub(crate) fn scalar_mul(&self, c: F) -> Self { - Self(self.iter().map(|&x| c * x).collect()) - } - - /// Removes leading zero coefficients. - pub fn trim(&mut self) { - self.0.drain(self.degree_plus_one()..); - } - - /// Polynomial addition. - pub fn add(&self, other: &Self) -> Self { - let (mut a, mut b) = (self.clone(), other.clone()); - match a.len().cmp(&b.len()) { - Ordering::Less => a.pad(b.len()), - Ordering::Greater => b.pad(a.len()), - _ => (), - } - Self(a.iter().zip(b.iter()).map(|(&x, &y)| x + y).collect()) - } - - /// Zero-pad the coefficients to have a given length. - pub fn pad(&mut self, len: usize) { - self.trim(); - assert!(self.len() <= len); - self.0.extend((self.len()..len).map(|_| F::ZERO)); - } - - /// Returns the zero-padded polynomial. - pub fn padded(&self, len: usize) -> Self { - let mut a = self.clone(); - a.pad(len); - a - } - - /// Polynomial multiplication. - pub fn mul(&self, b: &Self) -> Self { - if self.is_zero() || b.is_zero() { - return Self::zero(1); - } - let a_deg = self.degree(); - let b_deg = b.degree(); - let new_deg = (a_deg + b_deg + 1).next_power_of_two(); - let a_pad = self.padded(new_deg); - let b_pad = b.padded(new_deg); - - let precomputation = fft_precompute(new_deg); - let a_evals = fft_with_precomputation_power_of_2(a_pad.0.to_vec().into(), &precomputation); - let b_evals = fft_with_precomputation_power_of_2(b_pad.0.to_vec().into(), &precomputation); - - let mul_evals: Vec = a_evals - .values - .iter() - .zip(b_evals.values.iter()) - .map(|(&pa, &pb)| pa * pb) - .collect(); - ifft_with_precomputation_power_of_2(mul_evals.to_vec().into(), &precomputation) - .coeffs - .into() - } - - /// Polynomial long division. - /// Returns `(q,r)` the quotient and remainder of the polynomial division of `a` by `b`. - /// Generally slower that the equivalent function `Polynomial::polynomial_division`. - pub fn polynomial_long_division(&self, b: &Self) -> (Self, Self) { - // Trim b. - let mut b = b.clone(); - b.trim(); - - let (a_degree, b_degree) = (self.degree(), b.degree()); - if self.is_zero() { - (Self::zero(1), Self::empty()) - } else if b.is_zero() { - panic!("Division by zero polynomial"); - } else if a_degree < b_degree { - (Self::zero(1), self.clone()) - } else { - // Now we know that self.degree() >= divisor.degree(); - let mut quotient = Self::zero(a_degree - b_degree + 1); - let mut remainder = self.clone(); - // Can unwrap here because we know self is not zero. - let divisor_leading_inv = b.lead().inverse(); - while !remainder.is_zero() && remainder.degree() >= b_degree { - let cur_q_coeff = remainder.lead() * divisor_leading_inv; - let cur_q_degree = remainder.degree() - b_degree; - quotient[cur_q_degree] = cur_q_coeff; - - for (i, &div_coeff) in b.iter().enumerate() { - remainder[cur_q_degree + i] = - remainder[cur_q_degree + i] - (cur_q_coeff * div_coeff); - } - remainder.trim(); - } - (quotient, remainder) - } - } - - /// Computes the inverse of `self` modulo `x^n`. - fn inv_mod_xn(&self, n: usize) -> Self { - assert!(self[0].is_nonzero(), "Inverse doesn't exist."); - let mut h = self.clone(); - if h.len() < n { - h.pad(n); - } - let mut a = Self::empty(); - a.0.push(h[0].inverse()); - for i in 0..log2_ceil(n) { - let l = 1 << i; - let h0 = h[..l].to_vec().into(); - let mut h1: Polynomial = h[l..].to_vec().into(); - let mut c = a.mul(&h0); - if l == c.len() { - c = Self::zero(1); - } else { - c.drain(0..l); - } - h1.trim(); - let mut tmp = a.mul(&h1); - tmp = tmp.add(&c); - tmp.iter_mut().for_each(|x| *x = -(*x)); - tmp.trim(); - let mut b = a.mul(&tmp); - b.trim(); - if b.len() > l { - b.drain(l..); - } - a.0.extend_from_slice(&b[..]); - } - a.drain(n..); - a - } - - /// Polynomial division. - /// Returns `(q,r)` the quotient and remainder of the polynomial division of `a` by `b`. - /// Algorithm from http://people.csail.mit.edu/madhu/ST12/scribe/lect06.pdf - pub fn polynomial_division(&self, b: &Self) -> (Self, Self) { - let (a_degree, b_degree) = (self.degree(), b.degree()); - if self.is_zero() { - (Self::zero(1), Self::empty()) - } else if b.is_zero() { - panic!("Division by zero polynomial"); - } else if a_degree < b_degree { - (Self::zero(1), self.clone()) - } else if b_degree == 0 { - (self.scalar_mul(b[0].inverse()), Self::empty()) - } else { - let rev_b = b.rev(); - let rev_b_inv = rev_b.inv_mod_xn(a_degree - b_degree + 1); - let rev_q: Polynomial = rev_b_inv - .mul(&self.rev()[..=a_degree - b_degree].to_vec().into())[..=a_degree - b_degree] - .to_vec() - .into(); - let mut q = rev_q.rev(); - let mut qb = q.mul(b); - qb.pad(self.len()); - let mut r = self.add(&qb.neg()); - q.trim(); - r.trim(); - (q, r) - } - } - - // Divides a polynomial `a` by `Z_H = X^n - 1`. Assumes `Z_H | a`, otherwise result is meaningless. - pub fn divide_by_z_h(&self, n: usize) -> Self { - if self.is_zero() { - return self.clone(); - } - let mut a_trim = self.clone(); - a_trim.trim(); - let g = F::MULTIPLICATIVE_GROUP_GENERATOR; - let mut g_pow = F::ONE; - // Multiply the i-th coefficient of `a` by `g^i`. Then `new_a(w^j) = old_a(g.w^j)`. - a_trim.iter_mut().for_each(|x| { - *x = (*x) * g_pow; - g_pow = g * g_pow; - }); - let d = a_trim.degree(); - let root = F::primitive_root_of_unity(log2_ceil(d + 1)); - let precomputation = fft_precompute(d + 1); - // Equals to the evaluation of `a` on `{g.w^i}`. - let mut a_eval = a_trim.eval_domain(&precomputation); - // Compute the denominators `1/(g^n.w^(n*i) - 1)` using batch inversion. - let denominator_g = g.exp_usize(n); - let root_n = root.exp_usize(n); - let mut root_pow = F::ONE; - let denominators = (0..a_eval.len()) - .map(|i| { - if i != 0 { - root_pow = root_pow * root_n; - } - denominator_g * root_pow - F::ONE - }) - .collect::>(); - let denominators_inv = F::batch_multiplicative_inverse(&denominators); - // Divide every element of `a_eval` by the corresponding denominator. - // Then, `a_eval` is the evaluation of `a/Z_H` on `{g.w^i}`. - a_eval - .iter_mut() - .zip(denominators_inv.iter()) - .for_each(|(x, &d)| { - *x = (*x) * d; - }); - // `p` is the interpolating polynomial of `a_eval` on `{w^i}`. - let mut p = Self::from_evaluations(&a_eval, &precomputation); - // We need to scale it by `g^(-i)` to get the interpolating polynomial of `a_eval` on `{g.w^i}`, - // a.k.a `a/Z_H`. - let g_inv = g.inverse(); - let mut g_inv_pow = F::ONE; - p.iter_mut().for_each(|x| { - *x = (*x) * g_inv_pow; - g_inv_pow = g_inv_pow * g_inv; - }); - p - } -} - -#[cfg(test)] -mod test { - use super::*; - use crate::field::crandall_field::CrandallField; - use rand::{thread_rng, Rng}; - use std::time::Instant; - - #[test] - fn test_polynomial_multiplication() { - type F = CrandallField; - let mut rng = thread_rng(); - let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000)); - let a = Polynomial(F::rand_vec(a_deg)); - let b = Polynomial(F::rand_vec(b_deg)); - let m1 = a.mul(&b); - let m2 = a.mul(&b); - for _ in 0..1000 { - let x = F::rand(); - assert_eq!(m1.eval(x), a.eval(x) * b.eval(x)); - assert_eq!(m2.eval(x), a.eval(x) * b.eval(x)); - } - } - - #[test] - fn test_inv_mod_xn() { - type F = CrandallField; - let mut rng = thread_rng(); - let a_deg = rng.gen_range(1, 1_000); - let n = rng.gen_range(1, 1_000); - let a = Polynomial(F::rand_vec(a_deg)); - let b = a.inv_mod_xn(n); - let mut m = a.mul(&b); - m.drain(n..); - m.trim(); - assert_eq!( - m, - Polynomial(vec![F::ONE]), - "a: {:#?}, b:{:#?}, n:{:#?}, m:{:#?}", - a, - b, - n, - m - ); - } - - #[test] - fn test_polynomial_long_division() { - type F = CrandallField; - let mut rng = thread_rng(); - let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000)); - let a = Polynomial(F::rand_vec(a_deg)); - let b = Polynomial(F::rand_vec(b_deg)); - let (q, r) = a.polynomial_long_division(&b); - for _ in 0..1000 { - let x = F::rand(); - assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x)); - } - } - - #[test] - fn test_polynomial_division() { - type F = CrandallField; - let mut rng = thread_rng(); - let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000)); - let a = Polynomial(F::rand_vec(a_deg)); - let b = Polynomial(F::rand_vec(b_deg)); - let (q, r) = a.polynomial_division(&b); - for _ in 0..1000 { - let x = F::rand(); - assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x)); - } - } - - #[test] - fn test_polynomial_division_by_constant() { - type F = CrandallField; - let mut rng = thread_rng(); - let a_deg = rng.gen_range(1, 10_000); - let a = Polynomial(F::rand_vec(a_deg)); - let b = Polynomial::from(vec![F::rand()]); - let (q, r) = a.polynomial_division(&b); - for _ in 0..1000 { - let x = F::rand(); - assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x)); - } - } - - #[test] - fn test_division_by_z_h() { - type F = CrandallField; - let mut rng = thread_rng(); - let a_deg = rng.gen_range(1, 10_000); - let n = rng.gen_range(1, a_deg); - let mut a = Polynomial(F::rand_vec(a_deg)); - a.trim(); - let z_h = { - let mut z_h_vec = vec![F::ZERO; n + 1]; - z_h_vec[n] = F::ONE; - z_h_vec[0] = F::NEG_ONE; - Polynomial(z_h_vec) - }; - let m = a.mul(&z_h); - let now = Instant::now(); - let mut a_test = m.divide_by_z_h(n); - a_test.trim(); - println!("Division time: {:?}", now.elapsed()); - assert_eq!(a, a_test); - } - - #[test] - fn divide_zero_poly_by_z_h() { - let zero_poly = Polynomial::::empty(); - zero_poly.divide_by_z_h(16); - } - - // Test to see which polynomial division method is faster for divisions of the type - // `(X^n - 1)/(X - a) - #[test] - fn test_division_linear() { - type F = CrandallField; - let mut rng = thread_rng(); - let l = 14; - let n = 1 << l; - let g = F::primitive_root_of_unity(l); - let xn_minus_one = { - let mut xn_min_one_vec = vec![F::ZERO; n + 1]; - xn_min_one_vec[n] = F::ONE; - xn_min_one_vec[0] = F::NEG_ONE; - Polynomial(xn_min_one_vec) - }; - - let a = g.exp_usize(rng.gen_range(0, n)); - let denom = Polynomial(vec![-a, F::ONE]); - let now = Instant::now(); - xn_minus_one.polynomial_division(&denom); - println!("Division time: {:?}", now.elapsed()); - let now = Instant::now(); - xn_minus_one.polynomial_long_division(&denom); - println!("Division time: {:?}", now.elapsed()); - } - - #[test] - fn eq() { - type F = CrandallField; - assert_eq!(Polynomial::(vec![]), Polynomial(vec![])); - assert_eq!(Polynomial::(vec![F::ZERO]), Polynomial(vec![F::ZERO])); - assert_eq!(Polynomial::(vec![]), Polynomial(vec![F::ZERO])); - assert_eq!(Polynomial::(vec![F::ZERO]), Polynomial(vec![])); - assert_eq!( - Polynomial::(vec![F::ZERO]), - Polynomial(vec![F::ZERO, F::ZERO]) - ); - assert_eq!( - Polynomial::(vec![F::ONE]), - Polynomial(vec![F::ONE, F::ZERO]) - ); - assert_ne!(Polynomial::(vec![]), Polynomial(vec![F::ONE])); - assert_ne!( - Polynomial::(vec![F::ZERO]), - Polynomial(vec![F::ZERO, F::ONE]) - ); - assert_ne!( - Polynomial::(vec![F::ZERO]), - Polynomial(vec![F::ONE, F::ZERO]) - ); - } -} diff --git a/src/polynomial/polynomial.rs b/src/polynomial/polynomial.rs index 62b9e8b8..f0e9ee08 100644 --- a/src/polynomial/polynomial.rs +++ b/src/polynomial/polynomial.rs @@ -3,8 +3,7 @@ use std::ops::{Add, Mul, Sub}; use crate::field::fft::{fft, ifft}; use crate::field::field::Field; -use crate::polynomial::old_polynomial::Polynomial; -use crate::util::log2_strict; +use crate::util::{log2_ceil, log2_strict}; /// A polynomial in point-value form. /// @@ -49,7 +48,7 @@ impl From> for PolynomialValues { } /// A polynomial in coefficient form. -#[derive(Clone, Debug, Eq, PartialEq)] +#[derive(Clone, Debug)] pub struct PolynomialCoeffs { pub(crate) coeffs: Vec, } @@ -67,6 +66,10 @@ impl PolynomialCoeffs { PolynomialCoeffs { coeffs } } + pub(crate) fn empty() -> Self { + Self::new(Vec::new()) + } + pub(crate) fn zero(len: usize) -> Self { Self::new(vec![F::ZERO; len]) } @@ -112,6 +115,7 @@ impl PolynomialCoeffs { } pub(crate) fn padded(&self, new_len: usize) -> Self { + assert!(new_len >= self.len()); let mut coeffs = self.coeffs.clone(); coeffs.resize(new_len, F::ZERO); Self { coeffs } @@ -119,10 +123,16 @@ impl PolynomialCoeffs { /// Removes leading zero coefficients. pub fn trim(&mut self) { - self.coeffs.drain(self.degree_plus_one()..); + self.coeffs.truncate(self.degree_plus_one()); } - /// Degree of the polynomial + 1. + /// Removes leading zero coefficients. + pub fn trimmed(&self) -> Self { + let coeffs = self.coeffs[..self.degree_plus_one()].to_vec(); + Self { coeffs } + } + + /// Degree of the polynomial + 1, or 0 for a polynomial with no non-zero coefficients. pub(crate) fn degree_plus_one(&self) -> usize { (0usize..self.len()) .rev() @@ -130,11 +140,50 @@ impl PolynomialCoeffs { .map_or(0, |i| i + 1) } - pub(crate) fn div_rem(&self, rhs: &Self) -> (Self, Self) { - let lhs = Polynomial::from(self.clone()); - let rhs = Polynomial::from(rhs.clone()); - let (q, r) = lhs.polynomial_division(&rhs); - (q.into(), r.into()) + /// Leading coefficient. + pub fn lead(&self) -> F { + self.coeffs + .iter() + .rev() + .find(|x| x.is_nonzero()) + .map_or(F::ZERO, |x| *x) + } + + /// Computes the inverse of `self` modulo `x^n`. + pub(crate) fn inv_mod_xn(&self, n: usize) -> Self { + assert!(self.coeffs[0].is_nonzero(), "Inverse doesn't exist."); + let mut h = self.padded(n); + let mut a = Self::empty(); + a.coeffs.push(h.coeffs[0].inverse()); + for i in 0..log2_ceil(n) { + let l = 1 << i; + let h0 = h.coeffs[..l].to_vec().into(); + let mut h1: Self = h.coeffs[l..].to_vec().into(); + let mut c = &a * &h0; + if l == c.len() { + c = Self::zero(1); + } else { + c.coeffs.drain(0..l); + } + h1.trim(); + let mut tmp = a.mul(&h1); + tmp = tmp.add(&c); + tmp.coeffs.iter_mut().for_each(|x| *x = -(*x)); + tmp.trim(); + let mut b = &a * &tmp; + b.trim(); + if b.len() > l { + b.coeffs.drain(l..); + } + a.coeffs.extend_from_slice(&b.coeffs); + } + a.coeffs.drain(n..); + a + } + + /// Reverse the order of the coefficients, not taking into account the leading zero coefficients. + pub(crate) fn rev(&self) -> Self { + Self::new(self.trimmed().coeffs.into_iter().rev().collect()) } pub fn fft(self) -> PolynomialValues { @@ -152,6 +201,22 @@ impl PolynomialCoeffs { } } +impl PartialEq for PolynomialCoeffs { + fn eq(&self, other: &Self) -> bool { + let max_terms = self.coeffs.len().max(other.coeffs.len()); + for i in 0..max_terms { + let self_i = self.coeffs.get(i).cloned().unwrap_or(F::ZERO); + let other_i = other.coeffs.get(i).cloned().unwrap_or(F::ZERO); + if self_i != other_i { + return false; + } + } + true + } +} + +impl Eq for PolynomialCoeffs {} + impl From> for PolynomialCoeffs { fn from(coeffs: Vec) -> Self { Self::new(coeffs) @@ -215,10 +280,40 @@ impl Mul for &PolynomialCoeffs { #[cfg(test)] mod tests { + use std::time::Instant; + + use rand::{thread_rng, Rng}; + use crate::field::crandall_field::CrandallField; use super::*; + #[test] + fn test_trimmed() { + type F = CrandallField; + + assert_eq!( + PolynomialCoeffs:: { coeffs: vec![] }.trimmed(), + PolynomialCoeffs:: { coeffs: vec![] } + ); + assert_eq!( + PolynomialCoeffs:: { + coeffs: vec![F::ZERO] + } + .trimmed(), + PolynomialCoeffs:: { coeffs: vec![] } + ); + assert_eq!( + PolynomialCoeffs:: { + coeffs: vec![F::ONE, F::TWO, F::ZERO, F::ZERO] + } + .trimmed(), + PolynomialCoeffs:: { + coeffs: vec![F::ONE, F::TWO] + } + ); + } + #[test] fn test_coset_fft() { type F = CrandallField; @@ -237,4 +332,179 @@ mod tests { assert_eq!(coset_evals, naive_coset_evals); } + + #[test] + fn test_polynomial_multiplication() { + type F = CrandallField; + let mut rng = thread_rng(); + let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000)); + let a = PolynomialCoeffs::new(F::rand_vec(a_deg)); + let b = PolynomialCoeffs::new(F::rand_vec(b_deg)); + let m1 = &a * &b; + let m2 = &a * &b; + for _ in 0..1000 { + let x = F::rand(); + assert_eq!(m1.eval(x), a.eval(x) * b.eval(x)); + assert_eq!(m2.eval(x), a.eval(x) * b.eval(x)); + } + } + + #[test] + fn test_inv_mod_xn() { + type F = CrandallField; + let mut rng = thread_rng(); + let a_deg = rng.gen_range(1, 1_000); + let n = rng.gen_range(1, 1_000); + let a = PolynomialCoeffs::new(F::rand_vec(a_deg)); + let b = a.inv_mod_xn(n); + let mut m = &a * &b; + m.coeffs.drain(n..); + m.trim(); + assert_eq!( + m, + PolynomialCoeffs::new(vec![F::ONE]), + "a: {:#?}, b:{:#?}, n:{:#?}, m:{:#?}", + a, + b, + n, + m + ); + } + + #[test] + fn test_polynomial_long_division() { + type F = CrandallField; + let mut rng = thread_rng(); + let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000)); + let a = PolynomialCoeffs::new(F::rand_vec(a_deg)); + let b = PolynomialCoeffs::new(F::rand_vec(b_deg)); + let (q, r) = a.div_rem_long_division(&b); + for _ in 0..1000 { + let x = F::rand(); + assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x)); + } + } + + #[test] + fn test_polynomial_division() { + type F = CrandallField; + let mut rng = thread_rng(); + let (a_deg, b_deg) = (rng.gen_range(1, 10_000), rng.gen_range(1, 10_000)); + let a = PolynomialCoeffs::new(F::rand_vec(a_deg)); + let b = PolynomialCoeffs::new(F::rand_vec(b_deg)); + let (q, r) = a.div_rem(&b); + for _ in 0..1000 { + let x = F::rand(); + assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x)); + } + } + + #[test] + fn test_polynomial_division_by_constant() { + type F = CrandallField; + let mut rng = thread_rng(); + let a_deg = rng.gen_range(1, 10_000); + let a = PolynomialCoeffs::new(F::rand_vec(a_deg)); + let b = PolynomialCoeffs::from(vec![F::rand()]); + let (q, r) = a.div_rem(&b); + for _ in 0..1000 { + let x = F::rand(); + assert_eq!(a.eval(x), b.eval(x) * q.eval(x) + r.eval(x)); + } + } + + #[test] + fn test_division_by_z_h() { + type F = CrandallField; + let mut rng = thread_rng(); + let a_deg = rng.gen_range(1, 10_000); + let n = rng.gen_range(1, a_deg); + let mut a = PolynomialCoeffs::new(F::rand_vec(a_deg)); + a.trim(); + let z_h = { + let mut z_h_vec = vec![F::ZERO; n + 1]; + z_h_vec[n] = F::ONE; + z_h_vec[0] = F::NEG_ONE; + PolynomialCoeffs::new(z_h_vec) + }; + let m = &a * &z_h; + let now = Instant::now(); + let mut a_test = m.divide_by_z_h(n); + a_test.trim(); + println!("Division time: {:?}", now.elapsed()); + assert_eq!(a, a_test); + } + + #[test] + fn divide_zero_poly_by_z_h() { + let zero_poly = PolynomialCoeffs::::empty(); + zero_poly.divide_by_z_h(16); + } + + // Test to see which polynomial division method is faster for divisions of the type + // `(X^n - 1)/(X - a) + #[test] + fn test_division_linear() { + type F = CrandallField; + let mut rng = thread_rng(); + let l = 14; + let n = 1 << l; + let g = F::primitive_root_of_unity(l); + let xn_minus_one = { + let mut xn_min_one_vec = vec![F::ZERO; n + 1]; + xn_min_one_vec[n] = F::ONE; + xn_min_one_vec[0] = F::NEG_ONE; + PolynomialCoeffs::new(xn_min_one_vec) + }; + + let a = g.exp_usize(rng.gen_range(0, n)); + let denom = PolynomialCoeffs::new(vec![-a, F::ONE]); + let now = Instant::now(); + xn_minus_one.div_rem(&denom); + println!("Division time: {:?}", now.elapsed()); + let now = Instant::now(); + xn_minus_one.div_rem_long_division(&denom); + println!("Division time: {:?}", now.elapsed()); + } + + #[test] + fn eq() { + type F = CrandallField; + assert_eq!( + PolynomialCoeffs::::new(vec![]), + PolynomialCoeffs::new(vec![]) + ); + assert_eq!( + PolynomialCoeffs::::new(vec![F::ZERO]), + PolynomialCoeffs::new(vec![F::ZERO]) + ); + assert_eq!( + PolynomialCoeffs::::new(vec![]), + PolynomialCoeffs::new(vec![F::ZERO]) + ); + assert_eq!( + PolynomialCoeffs::::new(vec![F::ZERO]), + PolynomialCoeffs::new(vec![]) + ); + assert_eq!( + PolynomialCoeffs::::new(vec![F::ZERO]), + PolynomialCoeffs::new(vec![F::ZERO, F::ZERO]) + ); + assert_eq!( + PolynomialCoeffs::::new(vec![F::ONE]), + PolynomialCoeffs::new(vec![F::ONE, F::ZERO]) + ); + assert_ne!( + PolynomialCoeffs::::new(vec![]), + PolynomialCoeffs::new(vec![F::ONE]) + ); + assert_ne!( + PolynomialCoeffs::::new(vec![F::ZERO]), + PolynomialCoeffs::new(vec![F::ZERO, F::ONE]) + ); + assert_ne!( + PolynomialCoeffs::::new(vec![F::ZERO]), + PolynomialCoeffs::new(vec![F::ONE, F::ZERO]) + ); + } } diff --git a/src/prover.rs b/src/prover.rs index 434d0084..6eb2e7c5 100644 --- a/src/prover.rs +++ b/src/prover.rs @@ -10,7 +10,6 @@ use crate::generator::generate_partial_witness; use crate::plonk_challenger::Challenger; use crate::plonk_common::{eval_l_1, evaluate_gate_constraints, reduce_with_powers_multi}; use crate::polynomial::commitment::ListPolynomialCommitment; -use crate::polynomial::division::divide_by_z_h; use crate::polynomial::polynomial::{PolynomialCoeffs, PolynomialValues}; use crate::proof::Proof; use crate::util::transpose; @@ -106,7 +105,7 @@ pub(crate) fn prove( let mut all_quotient_poly_chunks = Vec::with_capacity(num_checks * quotient_degree); for vanishing_poly in vanishing_polys.into_iter() { let vanishing_poly_coeff = ifft(vanishing_poly); - let quotient_poly_coeff = divide_by_z_h(vanishing_poly_coeff, degree); + let quotient_poly_coeff = vanishing_poly_coeff.divide_by_z_h(degree); // Split t into degree-n chunks. let quotient_poly_coeff_chunks = quotient_poly_coeff.chunks(degree); all_quotient_poly_chunks.extend(quotient_poly_coeff_chunks); From 7f445686eedb85f6d47f929e9563ee4bb0e401f3 Mon Sep 17 00:00:00 2001 From: Daniel Lubarov Date: Thu, 13 May 2021 15:44:36 -0700 Subject: [PATCH 2/2] Tweaks --- src/polynomial/division.rs | 42 +++++++++++++++++++++++++++++++----- src/polynomial/polynomial.rs | 32 --------------------------- 2 files changed, 37 insertions(+), 37 deletions(-) diff --git a/src/polynomial/division.rs b/src/polynomial/division.rs index c72b798c..1dcc8425 100644 --- a/src/polynomial/division.rs +++ b/src/polynomial/division.rs @@ -1,16 +1,16 @@ use crate::field::fft::{fft, ifft}; use crate::field::field::Field; use crate::polynomial::polynomial::PolynomialCoeffs; -use crate::util::log2_strict; +use crate::util::{log2_strict, log2_ceil}; impl PolynomialCoeffs { /// Polynomial division. /// Returns `(q, r)`, the quotient and remainder of the polynomial division of `a` by `b`. pub(crate) fn div_rem(&self, b: &Self) -> (Self, Self) { let (a_degree_plug_1, b_degree_plus_1) = (self.degree_plus_one(), b.degree_plus_one()); - if self.is_zero() { + if a_degree_plug_1 == 0 { (Self::zero(1), Self::empty()) - } else if b.is_zero() { + } else if b_degree_plus_1 == 0 { panic!("Division by zero polynomial"); } else if a_degree_plug_1 < b_degree_plus_1 { (Self::zero(1), self.clone()) @@ -41,9 +41,9 @@ impl PolynomialCoeffs { let b = b.trimmed(); let (a_degree_plus_1, b_degree_plus_1) = (self.degree_plus_one(), b.degree_plus_one()); - if self.is_zero() { + if a_degree_plus_1 == 0 { (Self::zero(1), Self::empty()) - } else if b.is_zero() { + } else if b_degree_plus_1 == 0 { panic!("Division by zero polynomial"); } else if a_degree_plus_1 < b_degree_plus_1 { (Self::zero(1), self.clone()) @@ -124,6 +124,38 @@ impl PolynomialCoeffs { }); p } + + /// Computes the inverse of `self` modulo `x^n`. + pub(crate) fn inv_mod_xn(&self, n: usize) -> Self { + assert!(self.coeffs[0].is_nonzero(), "Inverse doesn't exist."); + let mut h = self.padded(n); + let mut a = Self::empty(); + a.coeffs.push(h.coeffs[0].inverse()); + for i in 0..log2_ceil(n) { + let l = 1 << i; + let h0 = h.coeffs[..l].to_vec().into(); + let mut h1: Self = h.coeffs[l..].to_vec().into(); + let mut c = &a * &h0; + if l == c.len() { + c = Self::zero(1); + } else { + c.coeffs.drain(0..l); + } + h1.trim(); + let mut tmp = &a * &h1; + tmp = &tmp + &c; + tmp.coeffs.iter_mut().for_each(|x| *x = -(*x)); + tmp.trim(); + let mut b = &a * &tmp; + b.trim(); + if b.len() > l { + b.coeffs.drain(l..); + } + a.coeffs.extend_from_slice(&b.coeffs); + } + a.coeffs.drain(n..); + a + } } #[cfg(test)] diff --git a/src/polynomial/polynomial.rs b/src/polynomial/polynomial.rs index f0e9ee08..630a4816 100644 --- a/src/polynomial/polynomial.rs +++ b/src/polynomial/polynomial.rs @@ -149,38 +149,6 @@ impl PolynomialCoeffs { .map_or(F::ZERO, |x| *x) } - /// Computes the inverse of `self` modulo `x^n`. - pub(crate) fn inv_mod_xn(&self, n: usize) -> Self { - assert!(self.coeffs[0].is_nonzero(), "Inverse doesn't exist."); - let mut h = self.padded(n); - let mut a = Self::empty(); - a.coeffs.push(h.coeffs[0].inverse()); - for i in 0..log2_ceil(n) { - let l = 1 << i; - let h0 = h.coeffs[..l].to_vec().into(); - let mut h1: Self = h.coeffs[l..].to_vec().into(); - let mut c = &a * &h0; - if l == c.len() { - c = Self::zero(1); - } else { - c.coeffs.drain(0..l); - } - h1.trim(); - let mut tmp = a.mul(&h1); - tmp = tmp.add(&c); - tmp.coeffs.iter_mut().for_each(|x| *x = -(*x)); - tmp.trim(); - let mut b = &a * &tmp; - b.trim(); - if b.len() > l { - b.coeffs.drain(l..); - } - a.coeffs.extend_from_slice(&b.coeffs); - } - a.coeffs.drain(n..); - a - } - /// Reverse the order of the coefficients, not taking into account the leading zero coefficients. pub(crate) fn rev(&self) -> Self { Self::new(self.trimmed().coeffs.into_iter().rev().collect())