Add Z terms in vanishing poly

Cargo.toml

@@ -8,6 +8,8 @@ edition = "2018"
 env_logger = "0.8.3"
 log = "0.4.14"
 num = "0.3"
+rand = "0.7.3"
+rand_chacha = "0.2.2"
 rayon = "1.5.0"
 unroll = "0.1.5"
 

src/circuit_builder.rs

@@ -14,6 +14,7 @@ use crate::polynomial::polynomial::PolynomialValues;
 use crate::target::Target;
 use crate::util::{log2_strict, transpose, transpose_poly_values};
 use crate::wire::Wire;
+use crate::partition::get_subgroup_shift;
 
 pub struct CircuitBuilder<F: Field> {
     pub(crate) config: CircuitConfig,
@@ -148,7 +149,7 @@ impl<F: Field> CircuitBuilder<F> {
     }
 
     fn sigma_vecs(&self) -> Vec<PolynomialValues<F>> {
-        vec![PolynomialValues::zero(self.gate_instances.len())] // TODO
+        vec![PolynomialValues::zero(self.gate_instances.len()); self.config.num_routed_wires] // TODO
     }
 
     /// Builds a "full circuit", with both prover and verifier data.
@@ -166,10 +167,11 @@ impl<F: Field> CircuitBuilder<F> {
 
         let sigma_vecs = self.sigma_vecs();
         let sigma_ldes = PolynomialValues::lde_multiple(sigma_vecs, self.config.rate_bits);
-        let sigmas_root = merkle_root_bit_rev_order(transpose_poly_values(sigma_ldes));
+        let sigma_ldes_t = transpose_poly_values(sigma_ldes);
+        let sigmas_root = merkle_root_bit_rev_order(sigma_ldes_t.clone());
 
         let generators = self.get_generators();
-        let prover_only = ProverOnlyCircuitData { generators, constant_ldes_t };
+        let prover_only = ProverOnlyCircuitData { generators, constant_ldes_t, sigma_ldes_t };
         let verifier_only = VerifierOnlyCircuitData {};
 
         // The HashSet of gates will have a non-deterministic order. When converting to a Vec, we
@@ -182,6 +184,10 @@ impl<F: Field> CircuitBuilder<F> {
             .max()
             .expect("No gates?");
 
+        let k_is = (0..self.config.num_routed_wires)
+            .map(get_subgroup_shift)
+            .collect();
+
         let common = CommonCircuitData {
             config: self.config,
             degree_bits: log2_strict(degree),
@@ -189,6 +195,7 @@ impl<F: Field> CircuitBuilder<F> {
             num_gate_constraints,
             constants_root,
             sigmas_root,
+            k_is,
         };
 
         info!("Building circuit took {}s", start.elapsed().as_secs_f32());

src/circuit_data.rs

@@ -81,6 +81,8 @@ impl<F: Field> VerifierCircuitData<F> {
 pub(crate) struct ProverOnlyCircuitData<F: Field> {
     pub generators: Vec<Box<dyn WitnessGenerator<F>>>,
     pub constant_ldes_t: Vec<Vec<F>>,
+    /// Transpose of LDEs of sigma polynomials (in the context of Plonk's permutation argument).
+    pub sigma_ldes_t: Vec<Vec<F>>,
 }
 
 /// Circuit data required by the verifier, but not the prover.
@@ -102,6 +104,9 @@ pub(crate) struct CommonCircuitData<F: Field> {
 
     /// A commitment to each permutation polynomial.
     pub(crate) sigmas_root: Hash<F>,
+
+    /// {k_i}. See `get_subgroup_shift`.
+    pub(crate) k_is: Vec<F>,
 }
 
 impl<F: Field> CommonCircuitData<F> {

src/field/crandall_field.rs

@@ -6,11 +6,6 @@ use num::Integer;
 
 use crate::field::field::Field;
 
-/// P = 2**64 - EPSILON
-///   = 2**64 - 9 * 2**28 + 1
-///   = 2**28 * (2**36 - 9) + 1
-const P: u64 = 18446744071293632513;
-
 /// EPSILON = 9 * 2**28 - 1
 const EPSILON: u64 = 2415919103;
 
@@ -19,6 +14,13 @@ const TWO_ADICITY: usize = 28;
 const POWER_OF_TWO_GENERATOR: CrandallField = CrandallField(10281950781551402419);
 
 /// A field designed for use with the Crandall reduction algorithm.
+///
+/// Its order is
+/// ```
+/// P = 2**64 - EPSILON
+///   = 2**64 - 9 * 2**28 + 1
+///   = 2**28 * (2**36 - 9) + 1
+/// ```
 // TODO: [Partial]Eq should compare canonical representations.
 #[derive(Copy, Clone)]
 pub struct CrandallField(pub u64);
@@ -47,8 +49,9 @@ impl Field for CrandallField {
     const ZERO: Self = Self(0);
     const ONE: Self = Self(1);
     const TWO: Self = Self(2);
-    const NEG_ONE: Self = Self(P - 1);
+    const NEG_ONE: Self = Self(Self::ORDER - 1);
 
+    const ORDER: u64 = 18446744071293632513;
     const MULTIPLICATIVE_SUBGROUP_GENERATOR: Self = Self(5); // TODO: Double check.
 
     #[inline(always)]
@@ -70,7 +73,7 @@ impl Field for CrandallField {
         // Applications to Cryptography".
 
         let mut u = self.0;
-        let mut v = P;
+        let mut v = Self::ORDER;
         let mut b = 1;
         let mut c = 0;
 
@@ -78,7 +81,7 @@ impl Field for CrandallField {
             while u.is_even() {
                 u >>= 1;
                 if b.is_odd() {
-                    b += P;
+                    b += Self::ORDER;
                 }
                 b >>= 1;
             }
@@ -86,7 +89,7 @@ impl Field for CrandallField {
             while v.is_even() {
                 v >>= 1;
                 if c.is_odd() {
-                    c += P;
+                    c += Self::ORDER;
                 }
                 c >>= 1;
             }
@@ -94,13 +97,13 @@ impl Field for CrandallField {
             if u < v {
                 v -= u;
                 if c < b {
-                    c += P;
+                    c += Self::ORDER;
                 }
                 c -= b;
             } else {
                 u -= v;
                 if b < c {
-                    b += P;
+                    b += Self::ORDER;
                 }
                 b -= c;
             }
@@ -145,8 +148,8 @@ impl Neg for CrandallField {
 
     #[inline]
     fn neg(self) -> Self {
-        let (diff, under) = P.overflowing_sub(self.0);
-        Self(diff.overflowing_add((under as u64) * P).0)
+        let (diff, under) = Self::ORDER.overflowing_sub(self.0);
+        Self(diff.overflowing_add((under as u64) * Self::ORDER).0)
     }
 }
 
@@ -156,7 +159,7 @@ impl Add for CrandallField {
     #[inline]
     fn add(self, rhs: Self) -> Self {
         let (sum, over) = self.0.overflowing_add(rhs.0);
-        Self(sum.overflowing_sub((over as u64) * P).0)
+        Self(sum.overflowing_sub((over as u64) * Self::ORDER).0)
     }
 }
 
@@ -172,7 +175,7 @@ impl Sub for CrandallField {
     #[inline]
     fn sub(self, rhs: Self) -> Self {
         let (diff, under) = self.0.overflowing_sub(rhs.0);
-        Self(diff.overflowing_add((under as u64) * P).0)
+        Self(diff.overflowing_add((under as u64) * Self::ORDER).0)
     }
 }
 

src/field/field.rs

@@ -1,5 +1,6 @@
 use std::fmt::{Debug, Display};
 use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
+use rand::Rng;
 
 /// A finite field with prime order less than 2^64.
 pub trait Field: 'static
@@ -23,11 +24,24 @@ pub trait Field: 'static
     const TWO: Self;
     const NEG_ONE: Self;
 
+    const ORDER: u64;
     const MULTIPLICATIVE_SUBGROUP_GENERATOR: Self;
 
-    fn sq(&self) -> Self;
+    fn is_zero(&self) -> bool {
+        *self == Self::ZERO
+    }
 
-    fn cube(&self) -> Self;
+    fn is_one(&self) -> bool {
+        *self == Self::ONE
+    }
+
+    fn sq(&self) -> Self {
+        *self * *self
+    }
+
+    fn cube(&self) -> Self {
+        *self * *self * *self
+    }
 
     /// Compute the multiplicative inverse of this field element.
     fn try_inverse(&self) -> Option<Self>;
@@ -97,4 +111,8 @@ pub trait Field: 'static
     fn exp_usize(&self, power: usize) -> Self {
         self.exp(Self::from_canonical_usize(power))
     }
+
+    fn rand_from_rng<R: Rng>(rng: &mut R) -> Self {
+        Self::from_canonical_u64(rng.gen_range(0, Self::ORDER))
+    }
 }

src/main.rs

@@ -26,6 +26,7 @@ mod gates;
 mod generator;
 mod gmimc;
 mod hash;
+mod partition;
 mod plonk_common;
 mod polynomial;
 mod proof;

src/partition.rs

@@ -0,0 +1,23 @@
+use rand::SeedableRng;
+use rand_chacha::ChaCha8Rng;
+
+use crate::field::field::Field;
+
+/// Returns `k_i`, the multiplier used in `S_ID_i` in the context of Plonk's permutation argument.
+// TODO: This is copied from plonky1, but may need revisiting. Is the random approach still OK now
+// that our field is smaller?
+pub(crate) fn get_subgroup_shift<F: Field>(i: usize) -> F {
+    // The optimized variant of Plonk's permutation argument calls for NUM_ROUTED_WIRES shifts,
+    // k_1, ..., k_n, which result in distinct cosets. The paper suggests a method which is
+    // fairly straightforward when only three shifts are needed, but seems a bit complex and
+    // expensive if more are needed.
+
+    // We will "cheat" and just use random field elements. Since our subgroup has |F*|/degree
+    // possible cosets, the probability of a collision is negligible for large fields.
+
+    // Unlike what's shown in the Plonk paper, we do not set k_1=1 to "randomize" the
+    // sigmas polynomials evaluations and making them fit in both fields with high probability.
+    // TODO: Go back to k_1=1 if we change the way we deal with values not fitting in both fields.
+    let mut rng = ChaCha8Rng::seed_from_u64(i as u64);
+    F::rand_from_rng(&mut rng)
+}

src/plonk_common.rs

@@ -2,6 +2,26 @@ use crate::circuit_builder::CircuitBuilder;
 use crate::field::field::Field;
 use crate::target::Target;
 
+/// Evaluate the polynomial which vanishes on any multiplicative subgroup of a given order `n`.
+pub(crate) fn eval_zero_poly<F: Field>(n: usize, x: F) -> F {
+    // Z(x) = x^n - 1
+    x.exp_usize(n) - F::ONE
+}
+
+/// Evaluate the Lagrange basis `L_1` with `L_1(1) = 1`, and `L_1(x) = 0` for other members of an
+/// order `n` multiplicative subgroup.
+pub(crate) fn eval_l_1<F: Field>(n: usize, x: F) -> F {
+    if x.is_one() {
+        // The code below would divide by zero, since we have (x - 1) in both the numerator and
+        // denominator.
+        return F::ONE;
+    }
+
+    // L_1(x) = (x^n - 1) / (n * (x - 1))
+    //        = Z(x) / (n * (x - 1))
+    eval_zero_poly(n, x) / (F::from_canonical_usize(n) * (x - F::ONE))
+}
+
 pub(crate) fn reduce_with_powers<F: Field>(terms: Vec<F>, alpha: F) -> F {
     let mut sum = F::ZERO;
     for &term in terms.iter().rev() {

src/prover.rs

@@ -9,7 +9,7 @@ use crate::field::fft::{fft, ifft};
 use crate::field::field::Field;
 use crate::generator::generate_partial_witness;
 use crate::hash::merkle_root_bit_rev_order;
-use crate::plonk_common::reduce_with_powers;
+use crate::plonk_common::{reduce_with_powers, eval_l_1};
 use crate::polynomial::division::divide_by_z_h;
 use crate::polynomial::polynomial::{PolynomialCoeffs, PolynomialValues};
 use crate::proof::Proof;
@@ -56,11 +56,13 @@ pub(crate) fn prove<F: Field>(
     let plonk_z_ldes_t = transpose_poly_values(plonk_z_ldes);
     let plonk_z_root = merkle_root_bit_rev_order(plonk_z_ldes_t.clone());
 
+    let beta = F::ZERO; // TODO
+    let gamma = F::ZERO; // TODO
     let alpha = F::ZERO; // TODO
 
     let start_vanishing_poly = Instant::now();
     let vanishing_poly = compute_vanishing_poly(
-        common_data, prover_data, wire_ldes_t, plonk_z_ldes_t, alpha);
+        common_data, prover_data, wire_ldes_t, plonk_z_ldes_t, beta, gamma, alpha);
     info!("Computing vanishing poly took {}s", start_vanishing_poly.elapsed().as_secs_f32());
 
     let quotient_poly_start = Instant::now();
@@ -102,6 +104,8 @@ fn compute_vanishing_poly<F: Field>(
     prover_data: &ProverOnlyCircuitData<F>,
     wire_ldes_t: Vec<Vec<F>>,
     plonk_z_lde_t: Vec<Vec<F>>,
+    beta: F,
+    gamma: F,
     alpha: F,
 ) -> PolynomialValues<F> {
     let lde_size = common_data.lde_size();
@@ -119,6 +123,7 @@ fn compute_vanishing_poly<F: Field>(
         let next_constants = &prover_data.constant_ldes_t[i_next];
         let local_plonk_zs = &plonk_z_lde_t[i];
         let next_plonk_zs = &plonk_z_lde_t[i_next];
+        let s_sigmas = &prover_data.sigma_ldes_t[i];
 
         debug_assert_eq!(local_wires.len(), common_data.config.num_wires);
         debug_assert_eq!(local_plonk_zs.len(), common_data.config.num_checks);
@@ -130,7 +135,7 @@ fn compute_vanishing_poly<F: Field>(
             next_wires,
         };
         result.push(compute_vanishing_poly_entry(
-            common_data, vars, local_plonk_zs, next_plonk_zs, alpha));
+            common_data, point, vars, local_plonk_zs, next_plonk_zs, s_sigmas, beta, gamma, alpha));
 
         point *= lde_gen;
     }
@@ -138,17 +143,52 @@ fn compute_vanishing_poly<F: Field>(
     PolynomialValues::new(result)
 }
 
+/// Evaluate the vanishing polynomial at `x`. In this context, the vanishing polynomial is a random
+/// linear combination of gate constraints, plus some other terms relating to the permutation
+/// argument. All such terms should vanish on `H`.
 fn compute_vanishing_poly_entry<F: Field>(
     common_data: &CommonCircuitData<F>,
+    x: F,
     vars: EvaluationVars<F>,
     local_plonk_zs: &[F],
     next_plonk_zs: &[F],
+    s_sigmas: &[F],
+    beta: F,
+    gamma: F,
     alpha: F,
 ) -> F {
-    let mut constraints = Vec::with_capacity(common_data.total_constraints());
-    // TODO: Add Z constraints.
-    constraints.extend(common_data.evaluate(vars));
-    reduce_with_powers(constraints, alpha)
+    let constraint_terms = common_data.evaluate(vars);
+
+    // The L_1(x) (Z(x) - 1) vanishing terms.
+    let mut vanishing_z_1_terms = Vec::new();
+    // The Z(x) f'(x) - g'(x) Z(g x) terms.
+    let mut vanishing_v_shift_terms = Vec::new();
+
+    for i in 0..common_data.config.num_checks {
+        let z_x = local_plonk_zs[i];
+        let z_gz = next_plonk_zs[i];
+        vanishing_z_1_terms.push(eval_l_1(common_data.degree(), x) * (z_x - F::ONE));
+
+        let mut f_prime = F::ONE;
+        let mut g_prime = F::ONE;
+        for j in 0..common_data.config.num_routed_wires {
+            let wire_value = vars.local_wires[j];
+            let k_i = common_data.k_is[j];
+            let s_id = k_i * x;
+            let s_sigma = s_sigmas[j];
+            f_prime *= wire_value + beta * s_id + gamma;
+            g_prime *= wire_value + beta * s_sigma + gamma;
+        }
+        vanishing_v_shift_terms.push(f_prime * z_x - g_prime * z_gz);
+    }
+
+    let vanishing_terms = [
+        vanishing_z_1_terms,
+        vanishing_v_shift_terms,
+        constraint_terms,
+    ].concat();
+
+    reduce_with_powers(vanishing_terms, alpha)
 }
 
 fn compute_wire_lde<F: Field>(