Batch all (#89)

- use sumcheck to batch open PCS
- split Prod and witness into two batches
- benchmark code

subroutines/Cargo.toml

@@ -0,0 +1,57 @@
+[package]
+name = "subroutines"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+
+ark-ff = { version = "^0.3.0", default-features = false }
+ark-std = { version = "^0.3.0", default-features = false }
+ark-ec = { version = "^0.3.0", default-features = false }
+ark-poly = { version = "^0.3.0", default-features = false }
+ark-serialize =  { version = "^0.3.0", default-features = false }
+ark-bls12-381 = { version = "0.3.0", default-features = false, features = [ "curve" ] }
+
+rand_chacha = { version = "0.3.0", default-features = false }
+displaydoc = { version = "0.2.3", default-features = false }
+rayon = { version = "1.5.2", default-features = false, optional = true }
+derivative = { version = "2", features = ["use_core"] }
+itertools = { version = "0.10.4", optional = true }
+
+transcript = { path = "../transcript" }
+arithmetic = { path = "../arithmetic" }
+util = { path = "../util" }
+
+# # Benchmarks
+# [[bench]]
+# name = "poly-iop-benches"
+# path = "benches/iop_bench.rs"
+# harness = false
+
+# Benchmarks
+[[bench]]
+name = "pcs-benches"
+path = "benches/pcs_bench.rs"
+harness = false
+
+[features]
+# default = [ "parallel", "print-trace" ]
+default = [ "parallel" ]
+# extensive sanity checks that are useful for debugging
+extensive_sanity_checks = [ ]
+parallel = [ 
+    "rayon",
+    "itertools",
+    "ark-std/parallel", 
+    "ark-ff/parallel",  
+    "ark-poly/parallel", 
+    "ark-ec/parallel",
+    "util/parallel",
+    "arithmetic/parallel",
+    ]
+print-trace = [ 
+    "arithmetic/print-trace",
+    "ark-std/print-trace",
+    ] 

subroutines/benches/iop_bench.rs

@@ -0,0 +1,270 @@
+use arithmetic::{identity_permutation_mle, VPAuxInfo, VirtualPolynomial};
+use ark_bls12_381::{Bls12_381, Fr};
+use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+use ark_std::test_rng;
+use std::{marker::PhantomData, rc::Rc, time::Instant};
+use subroutines::{
+    pcs::{prelude::MultilinearKzgPCS, PolynomialCommitmentScheme},
+    poly_iop::prelude::{
+        PermutationCheck, PolyIOP, PolyIOPErrors, ProductCheck, SumCheck, ZeroCheck,
+    },
+};
+
+type KZG = MultilinearKzgPCS<Bls12_381>;
+
+fn main() -> Result<(), PolyIOPErrors> {
+    bench_permutation_check()?;
+    println!("\n\n");
+    bench_sum_check()?;
+    println!("\n\n");
+    bench_prod_check()?;
+    println!("\n\n");
+    bench_zero_check()
+}
+
+fn bench_sum_check() -> Result<(), PolyIOPErrors> {
+    let mut rng = test_rng();
+    for degree in 2..4 {
+        for nv in 4..25 {
+            let repetition = if nv < 10 {
+                100
+            } else if nv < 20 {
+                50
+            } else {
+                10
+            };
+
+            let (poly, asserted_sum) =
+                VirtualPolynomial::rand(nv, (degree, degree + 1), 2, &mut rng)?;
+            let poly_info = poly.aux_info.clone();
+            let proof = {
+                let start = Instant::now();
+                for _ in 0..repetition {
+                    let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+                    let _proof = <PolyIOP<Fr> as SumCheck<Fr>>::prove(&poly, &mut transcript)?;
+                }
+
+                println!(
+                    "sum check proving time for {} variables and {} degree: {} ns",
+                    nv,
+                    degree,
+                    start.elapsed().as_nanos() / repetition as u128
+                );
+                let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+                <PolyIOP<Fr> as SumCheck<Fr>>::prove(&poly, &mut transcript)?
+            };
+
+            {
+                let start = Instant::now();
+
+                for _ in 0..repetition {
+                    let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+                    let _subclaim = <PolyIOP<Fr> as SumCheck<Fr>>::verify(
+                        asserted_sum,
+                        &proof,
+                        &poly_info,
+                        &mut transcript,
+                    )?;
+                }
+                println!(
+                    "sum check verification time for {} variables and {} degree: {} ns",
+                    nv,
+                    degree,
+                    start.elapsed().as_nanos() / repetition as u128
+                );
+            }
+
+            println!("====================================");
+        }
+    }
+    Ok(())
+}
+
+fn bench_zero_check() -> Result<(), PolyIOPErrors> {
+    let mut rng = test_rng();
+    for degree in 2..4 {
+        for nv in 4..20 {
+            let repetition = if nv < 10 {
+                100
+            } else if nv < 20 {
+                50
+            } else {
+                10
+            };
+
+            let poly = VirtualPolynomial::rand_zero(nv, (degree, degree + 1), 2, &mut rng)?;
+            let poly_info = poly.aux_info.clone();
+            let proof = {
+                let start = Instant::now();
+                let mut transcript = <PolyIOP<Fr> as ZeroCheck<Fr>>::init_transcript();
+                transcript.append_message(b"testing", b"initializing transcript for testing")?;
+                let proof = <PolyIOP<Fr> as ZeroCheck<Fr>>::prove(&poly, &mut transcript)?;
+
+                println!(
+                    "zero check proving time for {} variables and {} degree: {} ns",
+                    nv,
+                    degree,
+                    start.elapsed().as_nanos() / repetition as u128
+                );
+                proof
+            };
+
+            {
+                let start = Instant::now();
+                let mut transcript = <PolyIOP<Fr> as ZeroCheck<Fr>>::init_transcript();
+                transcript.append_message(b"testing", b"initializing transcript for testing")?;
+                let _zero_subclaim =
+                    <PolyIOP<Fr> as ZeroCheck<Fr>>::verify(&proof, &poly_info, &mut transcript)?;
+                println!(
+                    "zero check verification time for {} variables and {} degree: {} ns",
+                    nv,
+                    degree,
+                    start.elapsed().as_nanos() / repetition as u128
+                );
+            }
+
+            println!("====================================");
+        }
+    }
+    Ok(())
+}
+
+fn bench_permutation_check() -> Result<(), PolyIOPErrors> {
+    let mut rng = test_rng();
+
+    for nv in 4..20 {
+        let srs = KZG::gen_srs_for_testing(&mut rng, nv + 1)?;
+        let (pcs_param, _) = KZG::trim(&srs, None, Some(nv + 1))?;
+
+        let repetition = if nv < 10 {
+            100
+        } else if nv < 20 {
+            50
+        } else {
+            10
+        };
+
+        let w = Rc::new(DenseMultilinearExtension::rand(nv, &mut rng));
+
+        // s_perm is the identity map
+        let s_perm = identity_permutation_mle(nv);
+
+        let proof = {
+            let start = Instant::now();
+            let mut transcript =
+                <PolyIOP<Fr> as PermutationCheck<Bls12_381, KZG>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+
+            let (proof, _q_x) = <PolyIOP<Fr> as PermutationCheck<Bls12_381, KZG>>::prove(
+                &pcs_param,
+                &w,
+                &w,
+                &s_perm,
+                &mut transcript,
+            )?;
+
+            println!(
+                "permutation check proving time for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+            proof
+        };
+
+        {
+            let poly_info = VPAuxInfo {
+                max_degree: 2,
+                num_variables: nv,
+                phantom: PhantomData::default(),
+            };
+
+            let start = Instant::now();
+            let mut transcript =
+                <PolyIOP<Fr> as PermutationCheck<Bls12_381, KZG>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+            let _perm_check_sum_claim = <PolyIOP<Fr> as PermutationCheck<Bls12_381, KZG>>::verify(
+                &proof,
+                &poly_info,
+                &mut transcript,
+            )?;
+            println!(
+                "permutation check verification time for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+        }
+
+        println!("====================================");
+    }
+
+    Ok(())
+}
+
+fn bench_prod_check() -> Result<(), PolyIOPErrors> {
+    let mut rng = test_rng();
+
+    for nv in 4..20 {
+        let srs = KZG::gen_srs_for_testing(&mut rng, nv + 1)?;
+        let (pcs_param, _) = KZG::trim(&srs, None, Some(nv + 1))?;
+
+        let repetition = if nv < 10 {
+            100
+        } else if nv < 20 {
+            50
+        } else {
+            10
+        };
+
+        let f: DenseMultilinearExtension<Fr> = DenseMultilinearExtension::rand(nv, &mut rng);
+        let mut g = f.clone();
+        g.evaluations.reverse();
+        let f = Rc::new(f);
+        let g = Rc::new(g);
+
+        let proof = {
+            let start = Instant::now();
+            let mut transcript = <PolyIOP<Fr> as ProductCheck<Bls12_381, KZG>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+
+            let (proof, _prod_x) = <PolyIOP<Fr> as ProductCheck<Bls12_381, KZG>>::prove(
+                &pcs_param,
+                &f,
+                &g,
+                &mut transcript,
+            )?;
+
+            println!(
+                "product check proving time for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+            proof
+        };
+
+        {
+            let poly_info = VPAuxInfo {
+                max_degree: 2,
+                num_variables: nv,
+                phantom: PhantomData::default(),
+            };
+
+            let start = Instant::now();
+            let mut transcript = <PolyIOP<Fr> as ProductCheck<Bls12_381, KZG>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+            let _perm_check_sum_claim = <PolyIOP<Fr> as ProductCheck<Bls12_381, KZG>>::verify(
+                &proof,
+                &poly_info,
+                &mut transcript,
+            )?;
+            println!(
+                "product check verification time for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+        }
+
+        println!("====================================");
+    }
+
+    Ok(())
+}

subroutines/benches/pcs_bench.rs

@@ -0,0 +1,85 @@
+use ark_bls12_381::{Bls12_381, Fr};
+use ark_ff::UniformRand;
+use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+use ark_std::{rc::Rc, test_rng};
+use std::time::Instant;
+use subroutines::pcs::{
+    prelude::{MultilinearKzgPCS, PCSError, PolynomialCommitmentScheme},
+    StructuredReferenceString,
+};
+
+fn main() -> Result<(), PCSError> {
+    bench_pcs()
+}
+
+fn bench_pcs() -> Result<(), PCSError> {
+    let mut rng = test_rng();
+
+    // normal polynomials
+    let uni_params = MultilinearKzgPCS::<Bls12_381>::gen_srs_for_testing(&mut rng, 24)?;
+
+    for nv in 4..25 {
+        let repetition = if nv < 10 {
+            10
+        } else if nv < 20 {
+            5
+        } else {
+            2
+        };
+
+        let poly = Rc::new(DenseMultilinearExtension::rand(nv, &mut rng));
+        let (ck, vk) = uni_params.trim(nv)?;
+
+        let point: Vec<_> = (0..nv).map(|_| Fr::rand(&mut rng)).collect();
+
+        // commit
+        let com = {
+            let start = Instant::now();
+            for _ in 0..repetition {
+                let _commit = MultilinearKzgPCS::commit(&ck, &poly)?;
+            }
+
+            println!(
+                "KZG commit for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+
+            MultilinearKzgPCS::commit(&ck, &poly)?
+        };
+
+        // open
+        let (proof, value) = {
+            let start = Instant::now();
+            for _ in 0..repetition {
+                let _open = MultilinearKzgPCS::open(&ck, &poly, &point)?;
+            }
+
+            println!(
+                "KZG open for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+            MultilinearKzgPCS::open(&ck, &poly, &point)?
+        };
+
+        // verify
+        {
+            let start = Instant::now();
+            for _ in 0..repetition {
+                assert!(MultilinearKzgPCS::verify(
+                    &vk, &com, &point, &value, &proof
+                )?);
+            }
+            println!(
+                "KZG verify for {} variables: {} ns",
+                nv,
+                start.elapsed().as_nanos() / repetition as u128
+            );
+        }
+
+        println!("====================================");
+    }
+
+    Ok(())
+}

subroutines/src/lib.rs

@@ -0,0 +1,5 @@
+pub mod pcs;
+pub mod poly_iop;
+
+pub use pcs::prelude::*;
+pub use poly_iop::prelude::*;

subroutines/src/pcs/errors.rs

@@ -0,0 +1,50 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Error module.
+
+use arithmetic::ArithErrors;
+use ark_serialize::SerializationError;
+use ark_std::string::String;
+use displaydoc::Display;
+use transcript::TranscriptError;
+
+/// A `enum` specifying the possible failure modes of the PCS.
+#[derive(Display, Debug)]
+pub enum PCSError {
+    /// Invalid Prover: {0}
+    InvalidProver(String),
+    /// Invalid Verifier: {0}
+    InvalidVerifier(String),
+    /// Invalid Proof: {0}
+    InvalidProof(String),
+    /// Invalid parameters: {0}
+    InvalidParameters(String),
+    /// An error during (de)serialization: {0}
+    SerializationError(SerializationError),
+    /// Transcript error {0}
+    TranscriptError(TranscriptError),
+    /// ArithErrors error {0}
+    ArithErrors(ArithErrors),
+}
+
+impl From<SerializationError> for PCSError {
+    fn from(e: ark_serialize::SerializationError) -> Self {
+        Self::SerializationError(e)
+    }
+}
+
+impl From<TranscriptError> for PCSError {
+    fn from(e: TranscriptError) -> Self {
+        Self::TranscriptError(e)
+    }
+}
+
+impl From<ArithErrors> for PCSError {
+    fn from(e: ArithErrors) -> Self {
+        Self::ArithErrors(e)
+    }
+}

subroutines/src/pcs/mod.rs

@@ -0,0 +1,172 @@
+mod errors;
+mod multilinear_kzg;
+mod structs;
+mod univariate_kzg;
+
+pub mod prelude;
+
+use ark_ec::PairingEngine;
+use ark_ff::Field;
+use ark_serialize::{CanonicalDeserialize, CanonicalSerialize};
+use ark_std::rand::{CryptoRng, RngCore};
+use errors::PCSError;
+use std::{borrow::Borrow, fmt::Debug, hash::Hash};
+use transcript::IOPTranscript;
+
+/// This trait defines APIs for polynomial commitment schemes.
+/// Note that for our usage of PCS, we do not require the hiding property.
+pub trait PolynomialCommitmentScheme<E: PairingEngine> {
+    /// Prover parameters
+    type ProverParam: Clone;
+    /// Verifier parameters
+    type VerifierParam: Clone + CanonicalSerialize + CanonicalDeserialize;
+    /// Structured reference string
+    type SRS: Clone + Debug;
+    /// Polynomial and its associated types
+    type Polynomial: Clone
+        + Debug
+        + Hash
+        + PartialEq
+        + Eq
+        + CanonicalSerialize
+        + CanonicalDeserialize;
+    /// Polynomial input domain
+    type Point: Clone + Ord + Debug + Sync + Hash + PartialEq + Eq;
+    /// Polynomial Evaluation
+    type Evaluation: Field;
+    /// Commitments
+    type Commitment: Clone + CanonicalSerialize + CanonicalDeserialize + Debug + PartialEq + Eq;
+    /// Proofs
+    type Proof: Clone + CanonicalSerialize + CanonicalDeserialize + Debug + PartialEq + Eq;
+    /// Batch proofs
+    type BatchProof;
+
+    /// Build SRS for testing.
+    ///
+    /// - For univariate polynomials, `supported_size` is the maximum degree.
+    /// - For multilinear polynomials, `supported_size` is the number of
+    ///   variables.
+    ///
+    /// WARNING: THIS FUNCTION IS FOR TESTING PURPOSE ONLY.
+    /// THE OUTPUT SRS SHOULD NOT BE USED IN PRODUCTION.
+    fn gen_srs_for_testing<R: RngCore + CryptoRng>(
+        rng: &mut R,
+        supported_size: usize,
+    ) -> Result<Self::SRS, PCSError>;
+
+    /// Trim the universal parameters to specialize the public parameters.
+    /// Input both `supported_degree` for univariate and
+    /// `supported_num_vars` for multilinear.
+    /// ## Note on function signature
+    /// Usually, data structure like SRS and ProverParam are huge and users
+    /// might wish to keep them in heap using different kinds of smart pointers
+    /// (instead of only in stack) therefore our `impl Borrow<_>` interface
+    /// allows for passing in any pointer type, e.g.: `trim(srs: &Self::SRS,
+    /// ..)` or `trim(srs: Box<Self::SRS>, ..)` or `trim(srs: Arc<Self::SRS>,
+    /// ..)` etc.
+    fn trim(
+        srs: impl Borrow<Self::SRS>,
+        supported_degree: Option<usize>,
+        supported_num_vars: Option<usize>,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), PCSError>;
+
+    /// Generate a commitment for a polynomial
+    /// ## Note on function signature
+    /// Usually, data structure like SRS and ProverParam are huge and users
+    /// might wish to keep them in heap using different kinds of smart pointers
+    /// (instead of only in stack) therefore our `impl Borrow<_>` interface
+    /// allows for passing in any pointer type, e.g.: `commit(prover_param:
+    /// &Self::ProverParam, ..)` or `commit(prover_param:
+    /// Box<Self::ProverParam>, ..)` or `commit(prover_param:
+    /// Arc<Self::ProverParam>, ..)` etc.
+    fn commit(
+        prover_param: impl Borrow<Self::ProverParam>,
+        poly: &Self::Polynomial,
+    ) -> Result<Self::Commitment, PCSError>;
+
+    /// On input a polynomial `p` and a point `point`, outputs a proof for the
+    /// same.
+    fn open(
+        prover_param: impl Borrow<Self::ProverParam>,
+        polynomial: &Self::Polynomial,
+        point: &Self::Point,
+    ) -> Result<(Self::Proof, Self::Evaluation), PCSError>;
+
+    /// Input a list of multilinear extensions, and a same number of points, and
+    /// a transcript, compute a multi-opening for all the polynomials.
+    fn multi_open(
+        _prover_param: impl Borrow<Self::ProverParam>,
+        _polynomials: &[Self::Polynomial],
+        _points: &[Self::Point],
+        _evals: &[Self::Evaluation],
+        _transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<Self::BatchProof, PCSError> {
+        // the reason we use unimplemented!() is to enable developers to implement the
+        // trait without always implementing the batching APIs.
+        unimplemented!()
+    }
+
+    /// Verifies that `value` is the evaluation at `x` of the polynomial
+    /// committed inside `comm`.
+    fn verify(
+        verifier_param: &Self::VerifierParam,
+        commitment: &Self::Commitment,
+        point: &Self::Point,
+        value: &E::Fr,
+        proof: &Self::Proof,
+    ) -> Result<bool, PCSError>;
+
+    /// Verifies that `value_i` is the evaluation at `x_i` of the polynomial
+    /// `poly_i` committed inside `comm`.
+    fn batch_verify(
+        _verifier_param: &Self::VerifierParam,
+        _commitments: &[Self::Commitment],
+        _points: &[Self::Point],
+        _batch_proof: &Self::BatchProof,
+        _transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<bool, PCSError> {
+        // the reason we use unimplemented!() is to enable developers to implement the
+        // trait without always implementing the batching APIs.
+        unimplemented!()
+    }
+}
+
+/// API definitions for structured reference string
+pub trait StructuredReferenceString<E: PairingEngine>: Sized {
+    /// Prover parameters
+    type ProverParam;
+    /// Verifier parameters
+    type VerifierParam;
+
+    /// Extract the prover parameters from the public parameters.
+    fn extract_prover_param(&self, supported_size: usize) -> Self::ProverParam;
+    /// Extract the verifier parameters from the public parameters.
+    fn extract_verifier_param(&self, supported_size: usize) -> Self::VerifierParam;
+
+    /// Trim the universal parameters to specialize the public parameters
+    /// for polynomials to the given `supported_size`, and
+    /// returns committer key and verifier key.
+    ///
+    /// - For univariate polynomials, `supported_size` is the maximum degree.
+    /// - For multilinear polynomials, `supported_size` is 2 to the number of
+    ///   variables.
+    ///
+    /// `supported_log_size` should be in range `1..=params.log_size`
+    fn trim(
+        &self,
+        supported_size: usize,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), PCSError>;
+
+    /// Build SRS for testing.
+    ///
+    /// - For univariate polynomials, `supported_size` is the maximum degree.
+    /// - For multilinear polynomials, `supported_size` is the number of
+    ///   variables.
+    ///
+    /// WARNING: THIS FUNCTION IS FOR TESTING PURPOSE ONLY.
+    /// THE OUTPUT SRS SHOULD NOT BE USED IN PRODUCTION.
+    fn gen_srs_for_testing<R: RngCore + CryptoRng>(
+        rng: &mut R,
+        supported_size: usize,
+    ) -> Result<Self, PCSError>;
+}

subroutines/src/pcs/multilinear_kzg/batching.rs

@@ -0,0 +1,329 @@
+//! Sumcheck based batch opening and verify commitment.
+// TODO: refactoring this code to somewhere else
+// currently IOP depends on PCS because perm check requires commitment.
+// The sumcheck based batch opening therefore cannot stay in the PCS repo --
+// which creates a cyclic dependency.
+
+use crate::{
+    pcs::{
+        multilinear_kzg::util::eq_eval,
+        prelude::{Commitment, PCSError},
+        PolynomialCommitmentScheme,
+    },
+    poly_iop::{prelude::SumCheck, PolyIOP},
+    IOPProof,
+};
+use arithmetic::{
+    build_eq_x_r_vec, fix_last_variables, DenseMultilinearExtension, VPAuxInfo, VirtualPolynomial,
+};
+use ark_ec::{AffineCurve, PairingEngine, ProjectiveCurve};
+use ark_std::{end_timer, log2, start_timer, One, Zero};
+use std::{marker::PhantomData, rc::Rc};
+use transcript::IOPTranscript;
+
+#[derive(Clone, Debug, Default, PartialEq, Eq)]
+pub struct BatchProof<E, PCS>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<E>,
+{
+    /// A sum check proof proving tilde g's sum
+    pub(crate) sum_check_proof: IOPProof<E::Fr>,
+    /// f_i(point_i)
+    pub f_i_eval_at_point_i: Vec<E::Fr>,
+    /// proof for g'(a_2)
+    pub(crate) g_prime_proof: PCS::Proof,
+}
+
+/// Steps:
+/// 1. get challenge point t from transcript
+/// 2. build eq(t,i) for i in [0..k]
+/// 3. build \tilde g(i, b) = eq(t, i) * f_i(b)
+/// 4. compute \tilde eq
+/// 5. run sumcheck on \tilde eq * \tilde g(i, b)
+/// 6. build g'(a2) where (a1, a2) is the sumcheck's point
+pub(crate) fn multi_open_internal<E, PCS>(
+    prover_param: &PCS::ProverParam,
+    polynomials: &[PCS::Polynomial],
+    points: &[PCS::Point],
+    evals: &[PCS::Evaluation],
+    transcript: &mut IOPTranscript<E::Fr>,
+) -> Result<BatchProof<E, PCS>, PCSError>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<
+        E,
+        Polynomial = Rc<DenseMultilinearExtension<E::Fr>>,
+        Point = Vec<E::Fr>,
+        Evaluation = E::Fr,
+    >,
+{
+    let open_timer = start_timer!(|| format!("multi open {} points", points.len()));
+
+    // TODO: sanity checks
+    let num_var = polynomials[0].num_vars;
+    let k = polynomials.len();
+    let ell = log2(k) as usize;
+    let merged_num_var = num_var + ell;
+
+    // challenge point t
+    let t = transcript.get_and_append_challenge_vectors("t".as_ref(), ell)?;
+
+    // eq(t, i) for i in [0..k]
+    let eq_t_i_list = build_eq_x_r_vec(t.as_ref())?;
+
+    // \tilde g(i, b) = eq(t, i) * f_i(b)
+    let timer = start_timer!(|| format!("compute tilde g for {} points", points.len()));
+    let mut tilde_g_eval = vec![E::Fr::zero(); 1 << (ell + num_var)];
+    let block_size = 1 << num_var;
+    for (index, f_i) in polynomials.iter().enumerate() {
+        for (j, &f_i_eval) in f_i.iter().enumerate() {
+            tilde_g_eval[index * block_size + j] = f_i_eval * eq_t_i_list[index];
+        }
+    }
+    let tilde_g = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        merged_num_var,
+        tilde_g_eval,
+    ));
+    end_timer!(timer);
+
+    let timer = start_timer!(|| format!("compute tilde eq for {} points", points.len()));
+    let mut tilde_eq_eval = vec![E::Fr::zero(); 1 << (ell + num_var)];
+    for (index, point) in points.iter().enumerate() {
+        let eq_b_zi = build_eq_x_r_vec(point)?;
+        let start = index * block_size;
+        tilde_eq_eval[start..start + block_size].copy_from_slice(eq_b_zi.as_slice());
+    }
+    let tilde_eq = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        merged_num_var,
+        tilde_eq_eval,
+    ));
+    end_timer!(timer);
+
+    // built the virtual polynomial for SumCheck
+    let timer = start_timer!(|| format!("sum check prove of {} variables", num_var + ell));
+
+    let step = start_timer!(|| "add mle");
+    let mut sum_check_vp = VirtualPolynomial::new(num_var + ell);
+    sum_check_vp.add_mle_list([tilde_g.clone(), tilde_eq], E::Fr::one())?;
+    end_timer!(step);
+
+    let proof = match <PolyIOP<E::Fr> as SumCheck<E::Fr>>::prove(&sum_check_vp, transcript) {
+        Ok(p) => p,
+        Err(_e) => {
+            // cannot wrap IOPError with PCSError due to cyclic dependency
+            return Err(PCSError::InvalidProver(
+                "Sumcheck in batch proving Failed".to_string(),
+            ));
+        },
+    };
+
+    end_timer!(timer);
+
+    // (a1, a2) := sumcheck's point
+    let step = start_timer!(|| "open at a2");
+    let a1 = &proof.point[num_var..];
+    let a2 = &proof.point[..num_var];
+    end_timer!(step);
+
+    // build g'(a2)
+    let step = start_timer!(|| "evaluate at a2");
+    let g_prime = Rc::new(fix_last_variables(&tilde_g, a1));
+    end_timer!(step);
+
+    let step = start_timer!(|| "pcs open");
+    let (g_prime_proof, _g_prime_eval) = PCS::open(prover_param, &g_prime, a2.to_vec().as_ref())?;
+    // assert_eq!(g_prime_eval, tilde_g_eval);
+    end_timer!(step);
+
+    let step = start_timer!(|| "evaluate fi(pi)");
+    end_timer!(step);
+    end_timer!(open_timer);
+
+    Ok(BatchProof {
+        sum_check_proof: proof,
+        f_i_eval_at_point_i: evals.to_vec(),
+        g_prime_proof,
+    })
+}
+
+/// Steps:
+/// 1. get challenge point t from transcript
+/// 2. build g' commitment
+/// 3. ensure \sum_i eq(t, <i>) * f_i_evals matches the sum via SumCheck
+/// verification 4. verify commitment
+pub(crate) fn batch_verify_internal<E, PCS>(
+    verifier_param: &PCS::VerifierParam,
+    f_i_commitments: &[Commitment<E>],
+    points: &[PCS::Point],
+    proof: &BatchProof<E, PCS>,
+    transcript: &mut IOPTranscript<E::Fr>,
+) -> Result<bool, PCSError>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<
+        E,
+        Polynomial = Rc<DenseMultilinearExtension<E::Fr>>,
+        Point = Vec<E::Fr>,
+        Evaluation = E::Fr,
+        Commitment = Commitment<E>,
+    >,
+{
+    let open_timer = start_timer!(|| "batch verification");
+
+    // TODO: sanity checks
+
+    let k = f_i_commitments.len();
+    let ell = log2(k) as usize;
+    let num_var = proof.sum_check_proof.point.len() - ell;
+
+    // challenge point t
+    let t = transcript.get_and_append_challenge_vectors("t".as_ref(), ell)?;
+
+    // sum check point (a1, a2)
+    let a1 = &proof.sum_check_proof.point[num_var..];
+    let a2 = &proof.sum_check_proof.point[..num_var];
+
+    // build g' commitment
+    let eq_a1_list = build_eq_x_r_vec(a1)?;
+    let eq_t_list = build_eq_x_r_vec(t.as_ref())?;
+
+    let mut g_prime_commit = E::G1Affine::zero().into_projective();
+    for i in 0..k {
+        let tmp = eq_a1_list[i] * eq_t_list[i];
+        g_prime_commit += &f_i_commitments[i].0.mul(tmp);
+    }
+
+    // ensure \sum_i eq(t, <i>) * f_i_evals matches the sum via SumCheck
+    // verification
+    let mut sum = E::Fr::zero();
+    for (i, &e) in eq_t_list.iter().enumerate().take(k) {
+        sum += e * proof.f_i_eval_at_point_i[i];
+    }
+    let aux_info = VPAuxInfo {
+        max_degree: 2,
+        num_variables: num_var + ell,
+        phantom: PhantomData,
+    };
+    let subclaim = match <PolyIOP<E::Fr> as SumCheck<E::Fr>>::verify(
+        sum,
+        &proof.sum_check_proof,
+        &aux_info,
+        transcript,
+    ) {
+        Ok(p) => p,
+        Err(_e) => {
+            // cannot wrap IOPError with PCSError due to cyclic dependency
+            return Err(PCSError::InvalidProver(
+                "Sumcheck in batch verification failed".to_string(),
+            ));
+        },
+    };
+    let mut eq_tilde_eval = E::Fr::zero();
+    for (point, &coef) in points.iter().zip(eq_a1_list.iter()) {
+        eq_tilde_eval += coef * eq_eval(a2, point)?;
+    }
+    let tilde_g_eval = subclaim.expected_evaluation / eq_tilde_eval;
+
+    // verify commitment
+    let res = PCS::verify(
+        verifier_param,
+        &Commitment(g_prime_commit.into_affine()),
+        a2.to_vec().as_ref(),
+        &tilde_g_eval,
+        &proof.g_prime_proof,
+    )?;
+
+    end_timer!(open_timer);
+    Ok(res)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::pcs::{
+        prelude::{MultilinearKzgPCS, MultilinearUniversalParams},
+        StructuredReferenceString,
+    };
+    use arithmetic::get_batched_nv;
+    use ark_bls12_381::Bls12_381 as E;
+    use ark_ec::PairingEngine;
+    use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+    use ark_std::{
+        rand::{CryptoRng, RngCore},
+        test_rng,
+        vec::Vec,
+        UniformRand,
+    };
+
+    type Fr = <E as PairingEngine>::Fr;
+
+    fn test_multi_open_helper<R: RngCore + CryptoRng>(
+        ml_params: &MultilinearUniversalParams<E>,
+        polys: &[Rc<DenseMultilinearExtension<Fr>>],
+        rng: &mut R,
+    ) -> Result<(), PCSError> {
+        let merged_nv = get_batched_nv(polys[0].num_vars(), polys.len());
+        let (ml_ck, ml_vk) = ml_params.trim(merged_nv)?;
+
+        let mut points = Vec::new();
+        for poly in polys.iter() {
+            let point = (0..poly.num_vars())
+                .map(|_| Fr::rand(rng))
+                .collect::<Vec<Fr>>();
+            points.push(point);
+        }
+
+        let evals = polys
+            .iter()
+            .zip(points.iter())
+            .map(|(f, p)| f.evaluate(p).unwrap())
+            .collect::<Vec<_>>();
+
+        let commitments = polys
+            .iter()
+            .map(|poly| MultilinearKzgPCS::commit(&ml_ck.clone(), poly).unwrap())
+            .collect::<Vec<_>>();
+
+        let mut transcript = IOPTranscript::new("test transcript".as_ref());
+        transcript.append_field_element("init".as_ref(), &Fr::zero())?;
+
+        let batch_proof = multi_open_internal::<E, MultilinearKzgPCS<E>>(
+            &ml_ck,
+            polys,
+            &points,
+            &evals,
+            &mut transcript,
+        )?;
+
+        // good path
+        let mut transcript = IOPTranscript::new("test transcript".as_ref());
+        transcript.append_field_element("init".as_ref(), &Fr::zero())?;
+        assert!(batch_verify_internal::<E, MultilinearKzgPCS<E>>(
+            &ml_vk,
+            &commitments,
+            &points,
+            &batch_proof,
+            &mut transcript
+        )?);
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_multi_open_internal() -> Result<(), PCSError> {
+        let mut rng = test_rng();
+
+        let ml_params = MultilinearUniversalParams::<E>::gen_srs_for_testing(&mut rng, 20)?;
+        for num_poly in 5..6 {
+            for nv in 15..16 {
+                let polys1: Vec<_> = (0..num_poly)
+                    .map(|_| Rc::new(DenseMultilinearExtension::rand(nv, &mut rng)))
+                    .collect();
+                test_multi_open_helper(&ml_params, &polys1, &mut rng)?;
+            }
+        }
+
+        Ok(())
+    }
+}

subroutines/src/pcs/multilinear_kzg/mod.rs

@@ -0,0 +1,405 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Main module for multilinear KZG commitment scheme
+
+pub(crate) mod batching;
+pub(crate) mod srs;
+pub(crate) mod util;
+
+use crate::{
+    pcs::{prelude::Commitment, PCSError, PolynomialCommitmentScheme, StructuredReferenceString},
+    BatchProof,
+};
+use arithmetic::evaluate_opt;
+use ark_ec::{
+    msm::{FixedBaseMSM, VariableBaseMSM},
+    AffineCurve, PairingEngine, ProjectiveCurve,
+};
+use ark_ff::PrimeField;
+use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+use ark_serialize::{CanonicalDeserialize, CanonicalSerialize, Read, SerializationError, Write};
+use ark_std::{
+    borrow::Borrow,
+    end_timer, format,
+    marker::PhantomData,
+    rand::{CryptoRng, RngCore},
+    rc::Rc,
+    start_timer,
+    string::ToString,
+    vec,
+    vec::Vec,
+    One, Zero,
+};
+// use batching::{batch_verify_internal, multi_open_internal};
+use srs::{MultilinearProverParam, MultilinearUniversalParams, MultilinearVerifierParam};
+use transcript::IOPTranscript;
+
+use self::batching::{batch_verify_internal, multi_open_internal};
+
+/// KZG Polynomial Commitment Scheme on multilinear polynomials.
+pub struct MultilinearKzgPCS<E: PairingEngine> {
+    #[doc(hidden)]
+    phantom: PhantomData<E>,
+}
+
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug, PartialEq, Eq)]
+/// proof of opening
+pub struct MultilinearKzgProof<E: PairingEngine> {
+    /// Evaluation of quotients
+    pub proofs: Vec<E::G1Affine>,
+}
+
+impl<E: PairingEngine> PolynomialCommitmentScheme<E> for MultilinearKzgPCS<E> {
+    // Parameters
+    type ProverParam = MultilinearProverParam<E>;
+    type VerifierParam = MultilinearVerifierParam<E>;
+    type SRS = MultilinearUniversalParams<E>;
+    // Polynomial and its associated types
+    type Polynomial = Rc<DenseMultilinearExtension<E::Fr>>;
+    type Point = Vec<E::Fr>;
+    type Evaluation = E::Fr;
+    // Commitments and proofs
+    type Commitment = Commitment<E>;
+    type Proof = MultilinearKzgProof<E>;
+    type BatchProof = BatchProof<E, Self>;
+
+    /// Build SRS for testing.
+    ///
+    /// - For univariate polynomials, `log_size` is the log of maximum degree.
+    /// - For multilinear polynomials, `log_size` is the number of variables.
+    ///
+    /// WARNING: THIS FUNCTION IS FOR TESTING PURPOSE ONLY.
+    /// THE OUTPUT SRS SHOULD NOT BE USED IN PRODUCTION.
+    fn gen_srs_for_testing<R: RngCore + CryptoRng>(
+        rng: &mut R,
+        log_size: usize,
+    ) -> Result<Self::SRS, PCSError> {
+        MultilinearUniversalParams::<E>::gen_srs_for_testing(rng, log_size)
+    }
+
+    /// Trim the universal parameters to specialize the public parameters.
+    /// Input both `supported_log_degree` for univariate and
+    /// `supported_num_vars` for multilinear.
+    fn trim(
+        srs: impl Borrow<Self::SRS>,
+        supported_degree: Option<usize>,
+        supported_num_vars: Option<usize>,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), PCSError> {
+        assert!(supported_degree.is_none());
+
+        let supported_num_vars = match supported_num_vars {
+            Some(p) => p,
+            None => {
+                return Err(PCSError::InvalidParameters(
+                    "multilinear should receive a num_var param".to_string(),
+                ))
+            },
+        };
+        let (ml_ck, ml_vk) = srs.borrow().trim(supported_num_vars)?;
+
+        Ok((ml_ck, ml_vk))
+    }
+
+    /// Generate a commitment for a polynomial.
+    ///
+    /// This function takes `2^num_vars` number of scalar multiplications over
+    /// G1.
+    fn commit(
+        prover_param: impl Borrow<Self::ProverParam>,
+        poly: &Self::Polynomial,
+    ) -> Result<Self::Commitment, PCSError> {
+        let prover_param = prover_param.borrow();
+        let commit_timer = start_timer!(|| "commit");
+        if prover_param.num_vars < poly.num_vars {
+            return Err(PCSError::InvalidParameters(format!(
+                "MlE length ({}) exceeds param limit ({})",
+                poly.num_vars, prover_param.num_vars
+            )));
+        }
+        let ignored = prover_param.num_vars - poly.num_vars;
+        let scalars: Vec<_> = poly
+            .to_evaluations()
+            .into_iter()
+            .map(|x| x.into_repr())
+            .collect();
+        let commitment = VariableBaseMSM::multi_scalar_mul(
+            &prover_param.powers_of_g[ignored].evals,
+            scalars.as_slice(),
+        )
+        .into_affine();
+
+        end_timer!(commit_timer);
+        Ok(Commitment(commitment))
+    }
+
+    /// On input a polynomial `p` and a point `point`, outputs a proof for the
+    /// same. This function does not need to take the evaluation value as an
+    /// input.
+    ///
+    /// This function takes 2^{num_var +1} number of scalar multiplications over
+    /// G1:
+    /// - it prodceeds with `num_var` number of rounds,
+    /// - at round i, we compute an MSM for `2^{num_var - i + 1}` number of G2
+    ///   elements.
+    fn open(
+        prover_param: impl Borrow<Self::ProverParam>,
+        polynomial: &Self::Polynomial,
+        point: &Self::Point,
+    ) -> Result<(Self::Proof, Self::Evaluation), PCSError> {
+        open_internal(prover_param.borrow(), polynomial, point)
+    }
+
+    /// Input a list of multilinear extensions, and a same number of points, and
+    /// a transcript, compute a multi-opening for all the polynomials.
+    fn multi_open(
+        prover_param: impl Borrow<Self::ProverParam>,
+        polynomials: &[Self::Polynomial],
+        points: &[Self::Point],
+        evals: &[Self::Evaluation],
+        transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<BatchProof<E, Self>, PCSError> {
+        multi_open_internal(
+            prover_param.borrow(),
+            polynomials,
+            points,
+            evals,
+            transcript,
+        )
+    }
+
+    /// Verifies that `value` is the evaluation at `x` of the polynomial
+    /// committed inside `comm`.
+    ///
+    /// This function takes
+    /// - num_var number of pairing product.
+    /// - num_var number of MSM
+    fn verify(
+        verifier_param: &Self::VerifierParam,
+        commitment: &Self::Commitment,
+        point: &Self::Point,
+        value: &E::Fr,
+        proof: &Self::Proof,
+    ) -> Result<bool, PCSError> {
+        verify_internal(verifier_param, commitment, point, value, proof)
+    }
+
+    /// Verifies that `value_i` is the evaluation at `x_i` of the polynomial
+    /// `poly_i` committed inside `comm`.
+    fn batch_verify(
+        verifier_param: &Self::VerifierParam,
+        commitments: &[Self::Commitment],
+        points: &[Self::Point],
+        batch_proof: &Self::BatchProof,
+        transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<bool, PCSError> {
+        batch_verify_internal(verifier_param, commitments, points, batch_proof, transcript)
+    }
+}
+
+/// On input a polynomial `p` and a point `point`, outputs a proof for the
+/// same. This function does not need to take the evaluation value as an
+/// input.
+///
+/// This function takes 2^{num_var +1} number of scalar multiplications over
+/// G1:
+/// - it proceeds with `num_var` number of rounds,
+/// - at round i, we compute an MSM for `2^{num_var - i + 1}` number of G2
+///   elements.
+fn open_internal<E: PairingEngine>(
+    prover_param: &MultilinearProverParam<E>,
+    polynomial: &DenseMultilinearExtension<E::Fr>,
+    point: &[E::Fr],
+) -> Result<(MultilinearKzgProof<E>, E::Fr), PCSError> {
+    let open_timer = start_timer!(|| format!("open mle with {} variable", polynomial.num_vars));
+
+    if polynomial.num_vars() > prover_param.num_vars {
+        return Err(PCSError::InvalidParameters(format!(
+            "Polynomial num_vars {} exceed the limit {}",
+            polynomial.num_vars, prover_param.num_vars
+        )));
+    }
+
+    if polynomial.num_vars() != point.len() {
+        return Err(PCSError::InvalidParameters(format!(
+            "Polynomial num_vars {} does not match point len {}",
+            polynomial.num_vars,
+            point.len()
+        )));
+    }
+
+    let nv = polynomial.num_vars();
+    let ignored = prover_param.num_vars - nv;
+    let mut r: Vec<Vec<E::Fr>> = (0..nv + 1).map(|_| Vec::new()).collect();
+    let mut q: Vec<Vec<E::Fr>> = (0..nv + 1).map(|_| Vec::new()).collect();
+
+    r[nv] = polynomial.to_evaluations();
+
+    let mut proofs = Vec::new();
+
+    for (i, (&point_at_k, gi)) in point
+        .iter()
+        .zip(prover_param.powers_of_g[ignored..].iter())
+        .take(nv)
+        .enumerate()
+    {
+        let ith_round = start_timer!(|| format!("{}-th round", i));
+
+        let k = nv - i;
+        let cur_dim = 1 << (k - 1);
+        let mut cur_q = vec![E::Fr::zero(); cur_dim];
+        let mut cur_r = vec![E::Fr::zero(); cur_dim];
+        let one_minus_point_at_k = E::Fr::one() - point_at_k;
+
+        let ith_round_eval = start_timer!(|| format!("{}-th round eval", i));
+        for b in 0..(1 << (k - 1)) {
+            // q_b = pre_r [2^b + 1] - pre_r [2^b]
+            cur_q[b] = r[k][(b << 1) + 1] - r[k][b << 1];
+
+            // r_b = pre_r [2^b]*(1-p) + pre_r [2^b + 1] * p
+            cur_r[b] = r[k][b << 1] * one_minus_point_at_k + (r[k][(b << 1) + 1] * point_at_k);
+        }
+        end_timer!(ith_round_eval);
+        let scalars: Vec<_> = (0..(1 << k)).map(|x| cur_q[x >> 1].into_repr()).collect();
+
+        q[k] = cur_q;
+        r[k - 1] = cur_r;
+
+        // this is a MSM over G1 and is likely to be the bottleneck
+        proofs.push(VariableBaseMSM::multi_scalar_mul(&gi.evals, &scalars).into_affine());
+        end_timer!(ith_round);
+    }
+    let eval = evaluate_opt(polynomial, point);
+    end_timer!(open_timer);
+    Ok((MultilinearKzgProof { proofs }, eval))
+}
+
+/// Verifies that `value` is the evaluation at `x` of the polynomial
+/// committed inside `comm`.
+///
+/// This function takes
+/// - num_var number of pairing product.
+/// - num_var number of MSM
+fn verify_internal<E: PairingEngine>(
+    verifier_param: &MultilinearVerifierParam<E>,
+    commitment: &Commitment<E>,
+    point: &[E::Fr],
+    value: &E::Fr,
+    proof: &MultilinearKzgProof<E>,
+) -> Result<bool, PCSError> {
+    let verify_timer = start_timer!(|| "verify");
+    let num_var = point.len();
+
+    if num_var > verifier_param.num_vars {
+        return Err(PCSError::InvalidParameters(format!(
+            "point length ({}) exceeds param limit ({})",
+            num_var, verifier_param.num_vars
+        )));
+    }
+
+    let ignored = verifier_param.num_vars - num_var;
+    let prepare_inputs_timer = start_timer!(|| "prepare pairing inputs");
+
+    let scalar_size = E::Fr::size_in_bits();
+    let window_size = FixedBaseMSM::get_mul_window_size(num_var);
+
+    let h_table = FixedBaseMSM::get_window_table(
+        scalar_size,
+        window_size,
+        verifier_param.h.into_projective(),
+    );
+    let h_mul: Vec<E::G2Projective> =
+        FixedBaseMSM::multi_scalar_mul(scalar_size, window_size, &h_table, point);
+
+    let h_vec: Vec<_> = (0..num_var)
+        .map(|i| verifier_param.h_mask[ignored + i].into_projective() - h_mul[i])
+        .collect();
+    let h_vec: Vec<E::G2Affine> = E::G2Projective::batch_normalization_into_affine(&h_vec);
+    end_timer!(prepare_inputs_timer);
+
+    let pairing_product_timer = start_timer!(|| "pairing product");
+
+    let mut pairings: Vec<_> = proof
+        .proofs
+        .iter()
+        .map(|&x| E::G1Prepared::from(x))
+        .zip(h_vec.into_iter().take(num_var).map(E::G2Prepared::from))
+        .collect();
+
+    pairings.push((
+        E::G1Prepared::from(
+            (verifier_param.g.mul(*value) - commitment.0.into_projective()).into_affine(),
+        ),
+        E::G2Prepared::from(verifier_param.h),
+    ));
+
+    let res = E::product_of_pairings(pairings.iter()) == E::Fqk::one();
+
+    end_timer!(pairing_product_timer);
+    end_timer!(verify_timer);
+    Ok(res)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ark_bls12_381::Bls12_381;
+    use ark_ec::PairingEngine;
+    use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+    use ark_std::{rand::RngCore, test_rng, vec::Vec, UniformRand};
+
+    type E = Bls12_381;
+    type Fr = <E as PairingEngine>::Fr;
+
+    fn test_single_helper<R: RngCore + CryptoRng>(
+        params: &MultilinearUniversalParams<E>,
+        poly: &Rc<DenseMultilinearExtension<Fr>>,
+        rng: &mut R,
+    ) -> Result<(), PCSError> {
+        let nv = poly.num_vars();
+        assert_ne!(nv, 0);
+        let (ck, vk) = MultilinearKzgPCS::trim(params, None, Some(nv + 1))?;
+        let point: Vec<_> = (0..nv).map(|_| Fr::rand(rng)).collect();
+        let com = MultilinearKzgPCS::commit(&ck, poly)?;
+        let (proof, value) = MultilinearKzgPCS::open(&ck, poly, &point)?;
+
+        assert!(MultilinearKzgPCS::verify(
+            &vk, &com, &point, &value, &proof
+        )?);
+
+        let value = Fr::rand(rng);
+        assert!(!MultilinearKzgPCS::verify(
+            &vk, &com, &point, &value, &proof
+        )?);
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_single_commit() -> Result<(), PCSError> {
+        let mut rng = test_rng();
+
+        let params = MultilinearKzgPCS::<E>::gen_srs_for_testing(&mut rng, 10)?;
+
+        // normal polynomials
+        let poly1 = Rc::new(DenseMultilinearExtension::rand(8, &mut rng));
+        test_single_helper(&params, &poly1, &mut rng)?;
+
+        // single-variate polynomials
+        let poly2 = Rc::new(DenseMultilinearExtension::rand(1, &mut rng));
+        test_single_helper(&params, &poly2, &mut rng)?;
+
+        Ok(())
+    }
+
+    #[test]
+    fn setup_commit_verify_constant_polynomial() {
+        let mut rng = test_rng();
+
+        // normal polynomials
+        assert!(MultilinearKzgPCS::<E>::gen_srs_for_testing(&mut rng, 0).is_err());
+    }
+}

subroutines/src/pcs/multilinear_kzg/srs.rs

@@ -0,0 +1,253 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Implementing Structured Reference Strings for multilinear polynomial KZG
+use crate::pcs::{
+    multilinear_kzg::util::eq_extension, prelude::PCSError, StructuredReferenceString,
+};
+use ark_ec::{msm::FixedBaseMSM, AffineCurve, PairingEngine, ProjectiveCurve};
+use ark_ff::{Field, PrimeField};
+use ark_poly::DenseMultilinearExtension;
+use ark_serialize::{CanonicalDeserialize, CanonicalSerialize, Read, SerializationError, Write};
+use ark_std::{
+    collections::LinkedList,
+    end_timer, format,
+    rand::{CryptoRng, RngCore},
+    start_timer,
+    string::ToString,
+    vec::Vec,
+    UniformRand,
+};
+use core::iter::FromIterator;
+
+/// Evaluations over {0,1}^n for G1 or G2
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug)]
+pub struct Evaluations<C: AffineCurve> {
+    /// The evaluations.
+    pub evals: Vec<C>,
+}
+
+/// Universal Parameter
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug)]
+pub struct MultilinearUniversalParams<E: PairingEngine> {
+    /// prover parameters
+    pub prover_param: MultilinearProverParam<E>,
+    /// h^randomness: h^t1, h^t2, ..., **h^{t_nv}**
+    pub h_mask: Vec<E::G2Affine>,
+}
+
+/// Prover Parameters
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug)]
+pub struct MultilinearProverParam<E: PairingEngine> {
+    /// number of variables
+    pub num_vars: usize,
+    /// `pp_{num_vars}`, `pp_{num_vars - 1}`, `pp_{num_vars - 2}`, ..., defined
+    /// by XZZPD19
+    pub powers_of_g: Vec<Evaluations<E::G1Affine>>,
+    /// generator for G1
+    pub g: E::G1Affine,
+    /// generator for G2
+    pub h: E::G2Affine,
+}
+
+/// Verifier Parameters
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug)]
+pub struct MultilinearVerifierParam<E: PairingEngine> {
+    /// number of variables
+    pub num_vars: usize,
+    /// generator of G1
+    pub g: E::G1Affine,
+    /// generator of G2
+    pub h: E::G2Affine,
+    /// h^randomness: h^t1, h^t2, ..., **h^{t_nv}**
+    pub h_mask: Vec<E::G2Affine>,
+}
+
+impl<E: PairingEngine> StructuredReferenceString<E> for MultilinearUniversalParams<E> {
+    type ProverParam = MultilinearProverParam<E>;
+    type VerifierParam = MultilinearVerifierParam<E>;
+
+    /// Extract the prover parameters from the public parameters.
+    fn extract_prover_param(&self, supported_num_vars: usize) -> Self::ProverParam {
+        let to_reduce = self.prover_param.num_vars - supported_num_vars;
+
+        Self::ProverParam {
+            powers_of_g: self.prover_param.powers_of_g[to_reduce..].to_vec(),
+            g: self.prover_param.g,
+            h: self.prover_param.h,
+            num_vars: supported_num_vars,
+        }
+    }
+
+    /// Extract the verifier parameters from the public parameters.
+    fn extract_verifier_param(&self, supported_num_vars: usize) -> Self::VerifierParam {
+        let to_reduce = self.prover_param.num_vars - supported_num_vars;
+        Self::VerifierParam {
+            num_vars: supported_num_vars,
+            g: self.prover_param.g,
+            h: self.prover_param.h,
+            h_mask: self.h_mask[to_reduce..].to_vec(),
+        }
+    }
+
+    /// Trim the universal parameters to specialize the public parameters
+    /// for multilinear polynomials to the given `supported_num_vars`, and
+    /// returns committer key and verifier key. `supported_num_vars` should
+    /// be in range `1..=params.num_vars`
+    fn trim(
+        &self,
+        supported_num_vars: usize,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), PCSError> {
+        if supported_num_vars > self.prover_param.num_vars {
+            return Err(PCSError::InvalidParameters(format!(
+                "SRS does not support target number of vars {}",
+                supported_num_vars
+            )));
+        }
+
+        let to_reduce = self.prover_param.num_vars - supported_num_vars;
+        let ck = Self::ProverParam {
+            powers_of_g: self.prover_param.powers_of_g[to_reduce..].to_vec(),
+            g: self.prover_param.g,
+            h: self.prover_param.h,
+            num_vars: supported_num_vars,
+        };
+        let vk = Self::VerifierParam {
+            num_vars: supported_num_vars,
+            g: self.prover_param.g,
+            h: self.prover_param.h,
+            h_mask: self.h_mask[to_reduce..].to_vec(),
+        };
+        Ok((ck, vk))
+    }
+
+    /// Build SRS for testing.
+    /// WARNING: THIS FUNCTION IS FOR TESTING PURPOSE ONLY.
+    /// THE OUTPUT SRS SHOULD NOT BE USED IN PRODUCTION.
+    fn gen_srs_for_testing<R: RngCore + CryptoRng>(
+        rng: &mut R,
+        num_vars: usize,
+    ) -> Result<Self, PCSError> {
+        if num_vars == 0 {
+            return Err(PCSError::InvalidParameters(
+                "constant polynomial not supported".to_string(),
+            ));
+        }
+
+        let total_timer = start_timer!(|| "SRS generation");
+
+        let pp_generation_timer = start_timer!(|| "Prover Param generation");
+
+        let g = E::G1Projective::rand(rng);
+        let h = E::G2Projective::rand(rng);
+
+        let mut powers_of_g = Vec::new();
+
+        let t: Vec<_> = (0..num_vars).map(|_| E::Fr::rand(rng)).collect();
+        let scalar_bits = E::Fr::size_in_bits();
+
+        let mut eq: LinkedList<DenseMultilinearExtension<E::Fr>> =
+            LinkedList::from_iter(eq_extension(&t).into_iter());
+        let mut eq_arr = LinkedList::new();
+        let mut base = eq.pop_back().unwrap().evaluations;
+
+        for i in (0..num_vars).rev() {
+            eq_arr.push_front(remove_dummy_variable(&base, i)?);
+            if i != 0 {
+                let mul = eq.pop_back().unwrap().evaluations;
+                base = base
+                    .into_iter()
+                    .zip(mul.into_iter())
+                    .map(|(a, b)| a * b)
+                    .collect();
+            }
+        }
+
+        let mut pp_powers = Vec::new();
+        let mut total_scalars = 0;
+        for i in 0..num_vars {
+            let eq = eq_arr.pop_front().unwrap();
+            let pp_k_powers = (0..(1 << (num_vars - i))).map(|x| eq[x]);
+            pp_powers.extend(pp_k_powers);
+            total_scalars += 1 << (num_vars - i);
+        }
+        let window_size = FixedBaseMSM::get_mul_window_size(total_scalars);
+        let g_table = FixedBaseMSM::get_window_table(scalar_bits, window_size, g);
+
+        let pp_g = E::G1Projective::batch_normalization_into_affine(
+            &FixedBaseMSM::multi_scalar_mul(scalar_bits, window_size, &g_table, &pp_powers),
+        );
+
+        let mut start = 0;
+        for i in 0..num_vars {
+            let size = 1 << (num_vars - i);
+            let pp_k_g = Evaluations {
+                evals: pp_g[start..(start + size)].to_vec(),
+            };
+            powers_of_g.push(pp_k_g);
+            start += size;
+        }
+
+        let pp = Self::ProverParam {
+            num_vars,
+            g: g.into_affine(),
+            h: h.into_affine(),
+            powers_of_g,
+        };
+
+        end_timer!(pp_generation_timer);
+
+        let vp_generation_timer = start_timer!(|| "VP generation");
+        let h_mask = {
+            let window_size = FixedBaseMSM::get_mul_window_size(num_vars);
+            let h_table = FixedBaseMSM::get_window_table(scalar_bits, window_size, h);
+            E::G2Projective::batch_normalization_into_affine(&FixedBaseMSM::multi_scalar_mul(
+                scalar_bits,
+                window_size,
+                &h_table,
+                &t,
+            ))
+        };
+        end_timer!(vp_generation_timer);
+        end_timer!(total_timer);
+        Ok(Self {
+            prover_param: pp,
+            h_mask,
+        })
+    }
+}
+
+/// fix first `pad` variables of `poly` represented in evaluation form to zero
+fn remove_dummy_variable<F: Field>(poly: &[F], pad: usize) -> Result<Vec<F>, PCSError> {
+    if pad == 0 {
+        return Ok(poly.to_vec());
+    }
+    if !poly.len().is_power_of_two() {
+        return Err(PCSError::InvalidParameters(
+            "Size of polynomial should be power of two.".to_string(),
+        ));
+    }
+    let nv = ark_std::log2(poly.len()) as usize - pad;
+    Ok((0..(1 << nv)).map(|x| poly[x << pad]).collect())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ark_bls12_381::Bls12_381;
+    use ark_std::test_rng;
+    type E = Bls12_381;
+
+    #[test]
+    fn test_srs_gen() -> Result<(), PCSError> {
+        let mut rng = test_rng();
+        for nv in 4..10 {
+            let _ = MultilinearUniversalParams::<E>::gen_srs_for_testing(&mut rng, nv)?;
+        }
+
+        Ok(())
+    }
+}

subroutines/src/pcs/multilinear_kzg/util.rs

@@ -0,0 +1,51 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Useful utilities for KZG PCS
+use ark_ff::PrimeField;
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer, vec::Vec};
+
+use crate::PCSError;
+
+/// Generate eq(t,x), a product of multilinear polynomials with fixed t.
+/// eq(a,b) is takes extensions of a,b in {0,1}^num_vars such that if a and b in
+/// {0,1}^num_vars are equal then this polynomial evaluates to 1.
+pub(crate) fn eq_extension<F: PrimeField>(t: &[F]) -> Vec<DenseMultilinearExtension<F>> {
+    let start = start_timer!(|| "eq extension");
+
+    let dim = t.len();
+    let mut result = Vec::new();
+    for (i, &ti) in t.iter().enumerate().take(dim) {
+        let mut poly = Vec::with_capacity(1 << dim);
+        for x in 0..(1 << dim) {
+            let xi = if x >> i & 1 == 1 { F::one() } else { F::zero() };
+            let ti_xi = ti * xi;
+            poly.push(ti_xi + ti_xi - xi - ti + F::one());
+        }
+        result.push(DenseMultilinearExtension::from_evaluations_vec(dim, poly));
+    }
+
+    end_timer!(start);
+    result
+}
+
+/// Evaluate eq polynomial. use the public one later
+pub(crate) fn eq_eval<F: PrimeField>(x: &[F], y: &[F]) -> Result<F, PCSError> {
+    if x.len() != y.len() {
+        return Err(PCSError::InvalidParameters(
+            "x and y have different length".to_string(),
+        ));
+    }
+    let start = start_timer!(|| "eq_eval");
+    let mut res = F::one();
+    for (&xi, &yi) in x.iter().zip(y.iter()) {
+        let xi_yi = xi * yi;
+        res *= xi_yi + xi_yi - xi - yi + F::one();
+    }
+    end_timer!(start);
+    Ok(res)
+}

subroutines/src/pcs/prelude.rs

@@ -0,0 +1,21 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Prelude
+pub use crate::pcs::{
+    errors::PCSError,
+    multilinear_kzg::{
+        batching::BatchProof,
+        srs::{MultilinearProverParam, MultilinearUniversalParams, MultilinearVerifierParam},
+        MultilinearKzgPCS, MultilinearKzgProof,
+    },
+    structs::Commitment,
+    univariate_kzg::{
+        srs::{UnivariateProverParam, UnivariateUniversalParams, UnivariateVerifierParam},
+        UnivariateKzgBatchProof, UnivariateKzgPCS, UnivariateKzgProof,
+    },
+    PolynomialCommitmentScheme, StructuredReferenceString,
+};

subroutines/src/pcs/readme.md

@@ -0,0 +1,7 @@
+KZG based multilinear polynomial commitment
+-----
+
+# Compiling features:
+- `parallel`: use multi-threading when possible.
+- `print-trace`: print out user friendly information about the running time for each micro component.
+- `extensive_sanity_checks`: runs additional sanity checks that is not essential and will slow down the scheme.

subroutines/src/pcs/structs.rs

@@ -0,0 +1,25 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+use ark_ec::PairingEngine;
+use ark_serialize::{CanonicalDeserialize, CanonicalSerialize, Read, SerializationError, Write};
+use derivative::Derivative;
+
+#[derive(Derivative, CanonicalSerialize, CanonicalDeserialize)]
+#[derivative(
+    Default(bound = ""),
+    Hash(bound = ""),
+    Clone(bound = ""),
+    Copy(bound = ""),
+    Debug(bound = ""),
+    PartialEq(bound = ""),
+    Eq(bound = "")
+)]
+/// A commitment is an Affine point.
+pub struct Commitment<E: PairingEngine>(
+    /// the actual commitment is an affine point.
+    pub E::G1Affine,
+);

subroutines/src/pcs/univariate_kzg/mod.rs

@@ -0,0 +1,269 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Main module for univariate KZG commitment scheme
+
+use crate::pcs::{
+    prelude::Commitment, PCSError, PolynomialCommitmentScheme, StructuredReferenceString,
+};
+use ark_ec::{msm::VariableBaseMSM, AffineCurve, PairingEngine, ProjectiveCurve};
+use ark_ff::PrimeField;
+use ark_poly::{univariate::DensePolynomial, Polynomial, UVPolynomial};
+use ark_serialize::{CanonicalDeserialize, CanonicalSerialize, Read, SerializationError, Write};
+use ark_std::{
+    borrow::Borrow,
+    end_timer, format,
+    marker::PhantomData,
+    rand::{CryptoRng, RngCore},
+    start_timer,
+    string::ToString,
+    vec,
+    vec::Vec,
+    One,
+};
+use srs::{UnivariateProverParam, UnivariateUniversalParams, UnivariateVerifierParam};
+
+pub(crate) mod srs;
+
+/// KZG Polynomial Commitment Scheme on univariate polynomial.
+pub struct UnivariateKzgPCS<E: PairingEngine> {
+    #[doc(hidden)]
+    phantom: PhantomData<E>,
+}
+
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug, PartialEq, Eq)]
+/// proof of opening
+pub struct UnivariateKzgProof<E: PairingEngine> {
+    /// Evaluation of quotients
+    pub proof: E::G1Affine,
+}
+/// batch proof
+pub type UnivariateKzgBatchProof<E> = Vec<UnivariateKzgProof<E>>;
+
+impl<E: PairingEngine> PolynomialCommitmentScheme<E> for UnivariateKzgPCS<E> {
+    // Parameters
+    type ProverParam = UnivariateProverParam<E::G1Affine>;
+    type VerifierParam = UnivariateVerifierParam<E>;
+    type SRS = UnivariateUniversalParams<E>;
+    // Polynomial and its associated types
+    type Polynomial = DensePolynomial<E::Fr>;
+    type Point = E::Fr;
+    type Evaluation = E::Fr;
+    // Polynomial and its associated types
+    type Commitment = Commitment<E>;
+    type Proof = UnivariateKzgProof<E>;
+
+    // We do not implement batch univariate KZG at the current version.
+    type BatchProof = ();
+
+    /// Build SRS for testing.
+    ///
+    /// - For univariate polynomials, `supported_size` is the maximum degree.
+    ///
+    /// WARNING: THIS FUNCTION IS FOR TESTING PURPOSE ONLY.
+    /// THE OUTPUT SRS SHOULD NOT BE USED IN PRODUCTION.
+    fn gen_srs_for_testing<R: RngCore + CryptoRng>(
+        rng: &mut R,
+        supported_size: usize,
+    ) -> Result<Self::SRS, PCSError> {
+        Self::SRS::gen_srs_for_testing(rng, supported_size)
+    }
+
+    /// Trim the universal parameters to specialize the public parameters.
+    /// Input `max_degree` for univariate.
+    /// `supported_num_vars` must be None or an error is returned.
+    fn trim(
+        srs: impl Borrow<Self::SRS>,
+        supported_degree: Option<usize>,
+        supported_num_vars: Option<usize>,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), PCSError> {
+        assert!(supported_num_vars.is_none());
+        if supported_num_vars.is_some() {
+            return Err(PCSError::InvalidParameters(
+                "univariate should not receive a num_var param".to_string(),
+            ));
+        }
+        srs.borrow().trim(supported_degree.unwrap())
+    }
+
+    /// Generate a commitment for a polynomial
+    /// Note that the scheme is not hidding
+    fn commit(
+        prover_param: impl Borrow<Self::ProverParam>,
+        poly: &Self::Polynomial,
+    ) -> Result<Self::Commitment, PCSError> {
+        let prover_param = prover_param.borrow();
+        let commit_time =
+            start_timer!(|| format!("Committing to polynomial of degree {} ", poly.degree()));
+
+        if poly.degree() >= prover_param.powers_of_g.len() {
+            return Err(PCSError::InvalidParameters(format!(
+                "uni poly degree {} is larger than allowed {}",
+                poly.degree(),
+                prover_param.powers_of_g.len()
+            )));
+        }
+
+        let (num_leading_zeros, plain_coeffs) = skip_leading_zeros_and_convert_to_bigints(poly);
+
+        let msm_time = start_timer!(|| "MSM to compute commitment to plaintext poly");
+        let commitment = VariableBaseMSM::multi_scalar_mul(
+            &prover_param.powers_of_g[num_leading_zeros..],
+            &plain_coeffs,
+        )
+        .into_affine();
+        end_timer!(msm_time);
+
+        end_timer!(commit_time);
+        Ok(Commitment(commitment))
+    }
+
+    /// On input a polynomial `p` and a point `point`, outputs a proof for the
+    /// same.
+    fn open(
+        prover_param: impl Borrow<Self::ProverParam>,
+        polynomial: &Self::Polynomial,
+        point: &Self::Point,
+    ) -> Result<(Self::Proof, Self::Evaluation), PCSError> {
+        let open_time =
+            start_timer!(|| format!("Opening polynomial of degree {}", polynomial.degree()));
+        let divisor = Self::Polynomial::from_coefficients_vec(vec![-*point, E::Fr::one()]);
+
+        let witness_time = start_timer!(|| "Computing witness polynomial");
+        let witness_polynomial = polynomial / &divisor;
+        end_timer!(witness_time);
+
+        let (num_leading_zeros, witness_coeffs) =
+            skip_leading_zeros_and_convert_to_bigints(&witness_polynomial);
+
+        let proof = VariableBaseMSM::multi_scalar_mul(
+            &prover_param.borrow().powers_of_g[num_leading_zeros..],
+            &witness_coeffs,
+        )
+        .into_affine();
+
+        let eval = polynomial.evaluate(point);
+
+        end_timer!(open_time);
+        Ok((Self::Proof { proof }, eval))
+    }
+
+    /// Verifies that `value` is the evaluation at `x` of the polynomial
+    /// committed inside `comm`.
+    fn verify(
+        verifier_param: &Self::VerifierParam,
+        commitment: &Self::Commitment,
+        point: &Self::Point,
+        value: &E::Fr,
+        proof: &Self::Proof,
+    ) -> Result<bool, PCSError> {
+        let check_time = start_timer!(|| "Checking evaluation");
+        let pairing_inputs: Vec<(E::G1Prepared, E::G2Prepared)> = vec![
+            (
+                (verifier_param.g.mul(value.into_repr())
+                    - proof.proof.mul(point.into_repr())
+                    - commitment.0.into_projective())
+                .into_affine()
+                .into(),
+                verifier_param.h.into(),
+            ),
+            (proof.proof.into(), verifier_param.beta_h.into()),
+        ];
+
+        let res = E::product_of_pairings(pairing_inputs.iter()).is_one();
+
+        end_timer!(check_time, || format!("Result: {}", res));
+        Ok(res)
+    }
+}
+
+fn skip_leading_zeros_and_convert_to_bigints<F: PrimeField, P: UVPolynomial<F>>(
+    p: &P,
+) -> (usize, Vec<F::BigInt>) {
+    let mut num_leading_zeros = 0;
+    while num_leading_zeros < p.coeffs().len() && p.coeffs()[num_leading_zeros].is_zero() {
+        num_leading_zeros += 1;
+    }
+    let coeffs = convert_to_bigints(&p.coeffs()[num_leading_zeros..]);
+    (num_leading_zeros, coeffs)
+}
+
+fn convert_to_bigints<F: PrimeField>(p: &[F]) -> Vec<F::BigInt> {
+    let to_bigint_time = start_timer!(|| "Converting polynomial coeffs to bigints");
+    let coeffs = p.iter().map(|s| s.into_repr()).collect::<Vec<_>>();
+    end_timer!(to_bigint_time);
+    coeffs
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::StructuredReferenceString;
+    use ark_bls12_381::Bls12_381;
+    use ark_ec::PairingEngine;
+    use ark_poly::univariate::DensePolynomial;
+    use ark_std::{test_rng, UniformRand};
+
+    fn end_to_end_test_template<E>() -> Result<(), PCSError>
+    where
+        E: PairingEngine,
+    {
+        let rng = &mut test_rng();
+        for _ in 0..100 {
+            let mut degree = 0;
+            while degree <= 1 {
+                degree = usize::rand(rng) % 20;
+            }
+            let pp = UnivariateKzgPCS::<E>::gen_srs_for_testing(rng, degree)?;
+            let (ck, vk) = pp.trim(degree)?;
+            let p = <DensePolynomial<E::Fr> as UVPolynomial<E::Fr>>::rand(degree, rng);
+            let comm = UnivariateKzgPCS::<E>::commit(&ck, &p)?;
+            let point = E::Fr::rand(rng);
+            let (proof, value) = UnivariateKzgPCS::<E>::open(&ck, &p, &point)?;
+            assert!(
+                UnivariateKzgPCS::<E>::verify(&vk, &comm, &point, &value, &proof)?,
+                "proof was incorrect for max_degree = {}, polynomial_degree = {}",
+                degree,
+                p.degree(),
+            );
+        }
+        Ok(())
+    }
+
+    fn linear_polynomial_test_template<E>() -> Result<(), PCSError>
+    where
+        E: PairingEngine,
+    {
+        let rng = &mut test_rng();
+        for _ in 0..100 {
+            let degree = 50;
+
+            let pp = UnivariateKzgPCS::<E>::gen_srs_for_testing(rng, degree)?;
+            let (ck, vk) = pp.trim(degree)?;
+            let p = <DensePolynomial<E::Fr> as UVPolynomial<E::Fr>>::rand(degree, rng);
+            let comm = UnivariateKzgPCS::<E>::commit(&ck, &p)?;
+            let point = E::Fr::rand(rng);
+            let (proof, value) = UnivariateKzgPCS::<E>::open(&ck, &p, &point)?;
+            assert!(
+                UnivariateKzgPCS::<E>::verify(&vk, &comm, &point, &value, &proof)?,
+                "proof was incorrect for max_degree = {}, polynomial_degree = {}",
+                degree,
+                p.degree(),
+            );
+        }
+        Ok(())
+    }
+
+    #[test]
+    fn end_to_end_test() {
+        end_to_end_test_template::<Bls12_381>().expect("test failed for bls12-381");
+    }
+
+    #[test]
+    fn linear_polynomial_test() {
+        linear_polynomial_test_template::<Bls12_381>().expect("test failed for bls12-381");
+    }
+}

subroutines/src/pcs/univariate_kzg/srs.rs

@@ -0,0 +1,156 @@
+// Copyright (c) 2022 Espresso Systems (espressosys.com)
+// This file is part of the Jellyfish library.
+
+// You should have received a copy of the MIT License
+// along with the Jellyfish library. If not, see <https://mit-license.org/>.
+
+//! Implementing Structured Reference Strings for univariate polynomial KZG
+
+use crate::pcs::{PCSError, StructuredReferenceString};
+use ark_ec::{msm::FixedBaseMSM, AffineCurve, PairingEngine, ProjectiveCurve};
+use ark_ff::PrimeField;
+use ark_serialize::{CanonicalDeserialize, CanonicalSerialize, Read, SerializationError, Write};
+use ark_std::{
+    end_timer,
+    rand::{CryptoRng, RngCore},
+    start_timer, vec,
+    vec::Vec,
+    One, UniformRand,
+};
+use derivative::Derivative;
+
+/// `UniversalParams` are the universal parameters for the KZG10 scheme.
+// Adapted from
+// https://github.com/arkworks-rs/poly-commit/blob/master/src/kzg10/data_structures.rs#L20
+#[derive(Debug, Clone, Eq, PartialEq, CanonicalSerialize, CanonicalDeserialize, Default)]
+pub struct UnivariateUniversalParams<E: PairingEngine> {
+    /// Group elements of the form `{ \beta^i G }`, where `i` ranges from 0 to
+    /// `degree`.
+    pub powers_of_g: Vec<E::G1Affine>,
+    /// The generator of G2.
+    pub h: E::G2Affine,
+    /// \beta times the above generator of G2.
+    pub beta_h: E::G2Affine,
+}
+
+impl<E: PairingEngine> UnivariateUniversalParams<E> {
+    /// Returns the maximum supported degree
+    pub fn max_degree(&self) -> usize {
+        self.powers_of_g.len()
+    }
+}
+
+/// `UnivariateProverParam` is used to generate a proof
+#[derive(CanonicalSerialize, CanonicalDeserialize, Clone, Debug, Eq, PartialEq, Default)]
+pub struct UnivariateProverParam<C: AffineCurve> {
+    /// Parameters
+    pub powers_of_g: Vec<C>,
+}
+
+/// `UnivariateVerifierParam` is used to check evaluation proofs for a given
+/// commitment.
+#[derive(Derivative, CanonicalSerialize, CanonicalDeserialize)]
+#[derivative(
+    Default(bound = ""),
+    Clone(bound = ""),
+    Copy(bound = ""),
+    Debug(bound = ""),
+    PartialEq(bound = ""),
+    Eq(bound = "")
+)]
+pub struct UnivariateVerifierParam<E: PairingEngine> {
+    /// The generator of G1.
+    pub g: E::G1Affine,
+    /// The generator of G2.
+    pub h: E::G2Affine,
+    /// \beta times the above generator of G2.
+    pub beta_h: E::G2Affine,
+}
+
+impl<E: PairingEngine> StructuredReferenceString<E> for UnivariateUniversalParams<E> {
+    type ProverParam = UnivariateProverParam<E::G1Affine>;
+    type VerifierParam = UnivariateVerifierParam<E>;
+
+    /// Extract the prover parameters from the public parameters.
+    fn extract_prover_param(&self, supported_size: usize) -> Self::ProverParam {
+        let powers_of_g = self.powers_of_g[..=supported_size].to_vec();
+
+        Self::ProverParam { powers_of_g }
+    }
+
+    /// Extract the verifier parameters from the public parameters.
+    fn extract_verifier_param(&self, _supported_size: usize) -> Self::VerifierParam {
+        Self::VerifierParam {
+            g: self.powers_of_g[0],
+            h: self.h,
+            beta_h: self.beta_h,
+        }
+    }
+
+    /// Trim the universal parameters to specialize the public parameters
+    /// for univariate polynomials to the given `supported_size`, and
+    /// returns committer key and verifier key. `supported_size` should
+    /// be in range `1..params.len()`
+    fn trim(
+        &self,
+        supported_size: usize,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), PCSError> {
+        let powers_of_g = self.powers_of_g[..=supported_size].to_vec();
+
+        let pk = Self::ProverParam { powers_of_g };
+        let vk = Self::VerifierParam {
+            g: self.powers_of_g[0],
+            h: self.h,
+            beta_h: self.beta_h,
+        };
+        Ok((pk, vk))
+    }
+
+    /// Build SRS for testing.
+    /// WARNING: THIS FUNCTION IS FOR TESTING PURPOSE ONLY.
+    /// THE OUTPUT SRS SHOULD NOT BE USED IN PRODUCTION.
+    fn gen_srs_for_testing<R: RngCore + CryptoRng>(
+        rng: &mut R,
+        max_degree: usize,
+    ) -> Result<Self, PCSError> {
+        let setup_time = start_timer!(|| format!("KZG10::Setup with degree {}", max_degree));
+        let beta = E::Fr::rand(rng);
+        let g = E::G1Projective::rand(rng);
+        let h = E::G2Projective::rand(rng);
+
+        let mut powers_of_beta = vec![E::Fr::one()];
+
+        let mut cur = beta;
+        for _ in 0..max_degree {
+            powers_of_beta.push(cur);
+            cur *= &beta;
+        }
+
+        let window_size = FixedBaseMSM::get_mul_window_size(max_degree + 1);
+
+        let scalar_bits = E::Fr::size_in_bits();
+        let g_time = start_timer!(|| "Generating powers of G");
+        // TODO: parallelization
+        let g_table = FixedBaseMSM::get_window_table(scalar_bits, window_size, g);
+        let powers_of_g = FixedBaseMSM::multi_scalar_mul::<E::G1Projective>(
+            scalar_bits,
+            window_size,
+            &g_table,
+            &powers_of_beta,
+        );
+        end_timer!(g_time);
+
+        let powers_of_g = E::G1Projective::batch_normalization_into_affine(&powers_of_g);
+
+        let h = h.into_affine();
+        let beta_h = h.mul(beta).into_affine();
+
+        let pp = Self {
+            powers_of_g,
+            h,
+            beta_h,
+        };
+        end_timer!(setup_time);
+        Ok(pp)
+    }
+}

subroutines/src/poly_iop/errors.rs

@@ -0,0 +1,56 @@
+//! Error module.
+
+use crate::pcs::prelude::PCSError;
+use arithmetic::ArithErrors;
+use ark_std::string::String;
+use displaydoc::Display;
+use transcript::TranscriptError;
+
+/// A `enum` specifying the possible failure modes of the PolyIOP.
+#[derive(Display, Debug)]
+pub enum PolyIOPErrors {
+    /// Invalid Prover: {0}
+    InvalidProver(String),
+    /// Invalid Verifier: {0}
+    InvalidVerifier(String),
+    /// Invalid Proof: {0}
+    InvalidProof(String),
+    /// Invalid parameters: {0}
+    InvalidParameters(String),
+    /// Invalid challenge: {0}
+    InvalidChallenge(String),
+    /// Should not arrive to this point
+    ShouldNotArrive,
+    /// An error during (de)serialization: {0}
+    SerializationErrors(ark_serialize::SerializationError),
+    /// Transcript Error: {0}
+    TranscriptErrors(TranscriptError),
+    /// Arithmetic Error: {0}
+    ArithmeticErrors(ArithErrors),
+    /// PCS error {0}
+    PCSErrors(PCSError),
+}
+
+impl From<ark_serialize::SerializationError> for PolyIOPErrors {
+    fn from(e: ark_serialize::SerializationError) -> Self {
+        Self::SerializationErrors(e)
+    }
+}
+
+impl From<TranscriptError> for PolyIOPErrors {
+    fn from(e: TranscriptError) -> Self {
+        Self::TranscriptErrors(e)
+    }
+}
+
+impl From<ArithErrors> for PolyIOPErrors {
+    fn from(e: ArithErrors) -> Self {
+        Self::ArithmeticErrors(e)
+    }
+}
+
+impl From<PCSError> for PolyIOPErrors {
+    fn from(e: PCSError) -> Self {
+        Self::PCSErrors(e)
+    }
+}

subroutines/src/poly_iop/mod.rs

@@ -0,0 +1,30 @@
+use ark_ff::PrimeField;
+use std::marker::PhantomData;
+
+mod errors;
+mod perm_check;
+pub mod prelude;
+mod prod_check;
+mod structs;
+mod sum_check;
+mod utils;
+mod zero_check;
+
+#[derive(Clone, Debug, Default, Copy, PartialEq, Eq)]
+/// Struct for PolyIOP protocol.
+/// It has an associated type `F` that defines the prime field the multi-variate
+/// polynomial operates on.
+///
+/// An PolyIOP may be instantiated with one of the following:
+/// - SumCheck protocol.
+/// - ZeroCheck protocol.
+/// - PermutationCheck protocol.
+///
+/// Those individual protocol may have similar or identical APIs.
+/// The systematic way to invoke specific protocol is, for example
+///     `<PolyIOP<F> as SumCheck<F>>::prove()`
+pub struct PolyIOP<F: PrimeField> {
+    /// Associated field
+    #[doc(hidden)]
+    phantom: PhantomData<F>,
+}

subroutines/src/poly_iop/perm_check/mod.rs

@@ -0,0 +1,276 @@
+//! Main module for the Permutation Check protocol
+
+use self::util::computer_num_and_denom;
+use crate::{
+    pcs::PolynomialCommitmentScheme,
+    poly_iop::{errors::PolyIOPErrors, prelude::ProductCheck, PolyIOP},
+};
+use ark_ec::PairingEngine;
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer};
+use std::rc::Rc;
+use transcript::IOPTranscript;
+
+/// A permutation subclaim consists of
+/// - the SubClaim from the ProductCheck
+/// - Challenges beta and gamma
+#[derive(Clone, Debug, Default, PartialEq)]
+pub struct PermutationCheckSubClaim<E, PCS, PC>
+where
+    E: PairingEngine,
+    PC: ProductCheck<E, PCS>,
+    PCS: PolynomialCommitmentScheme<E>,
+{
+    /// the SubClaim from the ProductCheck
+    pub product_check_sub_claim: PC::ProductCheckSubClaim,
+    /// Challenges beta and gamma
+    pub challenges: (E::Fr, E::Fr),
+}
+
+pub mod util;
+
+/// A PermutationCheck w.r.t. `(f, g, perm)`
+/// proves that g is a permutation of f under
+/// permutation `perm`
+/// It is derived from ProductCheck.
+///
+/// A Permutation Check IOP takes the following steps:
+///
+/// Inputs:
+/// - f(x)
+/// - g(x)
+/// - permutation s_perm(x)
+pub trait PermutationCheck<E, PCS>: ProductCheck<E, PCS>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<E>,
+{
+    type PermutationCheckSubClaim;
+    type PermutationProof;
+
+    /// Initialize the system with a transcript
+    ///
+    /// This function is optional -- in the case where a PermutationCheck is
+    /// an building block for a more complex protocol, the transcript
+    /// may be initialized by this complex protocol, and passed to the
+    /// PermutationCheck prover/verifier.
+    fn init_transcript() -> Self::Transcript;
+
+    /// Inputs:
+    /// - f(x)
+    /// - g(x)
+    /// - permutation s_perm(x)
+    /// Outputs:
+    /// - a permutation check proof proving that g is a permutation of f under
+    ///   s_perm
+    /// - the product polynomial build during product check
+    ///
+    /// Cost: O(N)
+    fn prove(
+        pcs_param: &PCS::ProverParam,
+        fx: &Self::MultilinearExtension,
+        gx: &Self::MultilinearExtension,
+        s_perm: &Self::MultilinearExtension,
+        transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<(Self::PermutationProof, Self::MultilinearExtension), PolyIOPErrors>;
+
+    /// Verify that an MLE g(x) is a permutation of
+    /// MLE f(x) over a permutation given by s_perm.
+    fn verify(
+        proof: &Self::PermutationProof,
+        aux_info: &Self::VPAuxInfo,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::PermutationCheckSubClaim, PolyIOPErrors>;
+}
+
+impl<E, PCS> PermutationCheck<E, PCS> for PolyIOP<E::Fr>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<E, Polynomial = Rc<DenseMultilinearExtension<E::Fr>>>,
+{
+    type PermutationCheckSubClaim = PermutationCheckSubClaim<E, PCS, Self>;
+    type PermutationProof = Self::ProductCheckProof;
+
+    fn init_transcript() -> Self::Transcript {
+        IOPTranscript::<E::Fr>::new(b"Initializing PermutationCheck transcript")
+    }
+
+    fn prove(
+        pcs_param: &PCS::ProverParam,
+        fx: &Self::MultilinearExtension,
+        gx: &Self::MultilinearExtension,
+        s_perm: &Self::MultilinearExtension,
+        transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<(Self::PermutationProof, Self::MultilinearExtension), PolyIOPErrors> {
+        let start = start_timer!(|| "Permutation check prove");
+        if fx.num_vars != gx.num_vars {
+            return Err(PolyIOPErrors::InvalidParameters(
+                "fx and gx have different number of variables".to_string(),
+            ));
+        }
+
+        if fx.num_vars != s_perm.num_vars {
+            return Err(PolyIOPErrors::InvalidParameters(
+                "fx and s_perm have different number of variables".to_string(),
+            ));
+        }
+
+        // generate challenge `beta` and `gamma` from current transcript
+        let beta = transcript.get_and_append_challenge(b"beta")?;
+        let gamma = transcript.get_and_append_challenge(b"gamma")?;
+        let (numerator, denominator) = computer_num_and_denom(&beta, &gamma, fx, gx, s_perm)?;
+
+        // invoke product check on numerator and denominator
+        let (proof, prod_poly) =
+            <Self as ProductCheck<E, PCS>>::prove(pcs_param, &numerator, &denominator, transcript)?;
+
+        end_timer!(start);
+        Ok((proof, prod_poly))
+    }
+
+    /// Verify that an MLE g(x) is a permutation of an
+    /// MLE f(x) over a permutation given by s_perm.
+    fn verify(
+        proof: &Self::PermutationProof,
+        aux_info: &Self::VPAuxInfo,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::PermutationCheckSubClaim, PolyIOPErrors> {
+        let start = start_timer!(|| "Permutation check verify");
+
+        let beta = transcript.get_and_append_challenge(b"beta")?;
+        let gamma = transcript.get_and_append_challenge(b"gamma")?;
+
+        // invoke the zero check on the iop_proof
+        let product_check_sub_claim =
+            <Self as ProductCheck<E, PCS>>::verify(proof, aux_info, transcript)?;
+
+        end_timer!(start);
+        Ok(PermutationCheckSubClaim {
+            product_check_sub_claim,
+            challenges: (beta, gamma),
+        })
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::PermutationCheck;
+    use crate::{
+        pcs::{prelude::MultilinearKzgPCS, PolynomialCommitmentScheme},
+        poly_iop::{errors::PolyIOPErrors, PolyIOP},
+    };
+    use arithmetic::{evaluate_opt, identity_permutation_mle, random_permutation_mle, VPAuxInfo};
+    use ark_bls12_381::Bls12_381;
+    use ark_ec::PairingEngine;
+    use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+    use ark_std::test_rng;
+    use std::{marker::PhantomData, rc::Rc};
+
+    type KZG = MultilinearKzgPCS<Bls12_381>;
+
+    fn test_permutation_check_helper<E, PCS>(
+        pcs_param: &PCS::ProverParam,
+        fx: &Rc<DenseMultilinearExtension<E::Fr>>,
+        gx: &Rc<DenseMultilinearExtension<E::Fr>>,
+        s_perm: &Rc<DenseMultilinearExtension<E::Fr>>,
+    ) -> Result<(), PolyIOPErrors>
+    where
+        E: PairingEngine,
+        PCS: PolynomialCommitmentScheme<E, Polynomial = Rc<DenseMultilinearExtension<E::Fr>>>,
+    {
+        let nv = fx.num_vars;
+        let poly_info = VPAuxInfo {
+            max_degree: 2,
+            num_variables: nv,
+            phantom: PhantomData::default(),
+        };
+
+        // prover
+        let mut transcript = <PolyIOP<E::Fr> as PermutationCheck<E, PCS>>::init_transcript();
+        transcript.append_message(b"testing", b"initializing transcript for testing")?;
+        let (proof, prod_x) = <PolyIOP<E::Fr> as PermutationCheck<E, PCS>>::prove(
+            pcs_param,
+            fx,
+            gx,
+            s_perm,
+            &mut transcript,
+        )?;
+
+        // verifier
+        let mut transcript = <PolyIOP<E::Fr> as PermutationCheck<E, PCS>>::init_transcript();
+        transcript.append_message(b"testing", b"initializing transcript for testing")?;
+        let perm_check_sub_claim = <PolyIOP<E::Fr> as PermutationCheck<E, PCS>>::verify(
+            &proof,
+            &poly_info,
+            &mut transcript,
+        )?;
+
+        // check product subclaim
+        if evaluate_opt(
+            &prod_x,
+            &perm_check_sub_claim.product_check_sub_claim.final_query.0,
+        ) != perm_check_sub_claim.product_check_sub_claim.final_query.1
+        {
+            return Err(PolyIOPErrors::InvalidVerifier("wrong subclaim".to_string()));
+        };
+
+        Ok(())
+    }
+
+    fn test_permutation_check(nv: usize) -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+
+        let srs = MultilinearKzgPCS::<Bls12_381>::gen_srs_for_testing(&mut rng, nv + 1)?;
+        let (pcs_param, _) = MultilinearKzgPCS::<Bls12_381>::trim(&srs, None, Some(nv + 1))?;
+
+        {
+            // good path: w is a permutation of w itself under the identify map
+            let w = Rc::new(DenseMultilinearExtension::rand(nv, &mut rng));
+            // s_perm is the identity map
+            let s_perm = identity_permutation_mle(nv);
+            test_permutation_check_helper::<Bls12_381, KZG>(&pcs_param, &w, &w, &s_perm)?;
+        }
+
+        {
+            // bad path 1: w is a not permutation of w itself under a random map
+            let w = Rc::new(DenseMultilinearExtension::rand(nv, &mut rng));
+            // s_perm is a random map
+            let s_perm = random_permutation_mle(nv, &mut rng);
+
+            assert!(
+                test_permutation_check_helper::<Bls12_381, KZG>(&pcs_param, &w, &w, &s_perm)
+                    .is_err()
+            );
+        }
+
+        {
+            // bad path 2: f is a not permutation of g under a identity map
+            let f = Rc::new(DenseMultilinearExtension::rand(nv, &mut rng));
+            let g = Rc::new(DenseMultilinearExtension::rand(nv, &mut rng));
+            // s_perm is the identity map
+            let s_perm = identity_permutation_mle(nv);
+
+            assert!(
+                test_permutation_check_helper::<Bls12_381, KZG>(&pcs_param, &f, &g, &s_perm)
+                    .is_err()
+            );
+        }
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_trivial_polynomial() -> Result<(), PolyIOPErrors> {
+        test_permutation_check(1)
+    }
+    #[test]
+    fn test_normal_polynomial() -> Result<(), PolyIOPErrors> {
+        test_permutation_check(5)
+    }
+
+    #[test]
+    fn zero_polynomial_should_error() -> Result<(), PolyIOPErrors> {
+        assert!(test_permutation_check(0).is_err());
+        Ok(())
+    }
+}

subroutines/src/poly_iop/perm_check/util.rs

@@ -0,0 +1,61 @@
+//! This module implements useful functions for the permutation check protocol.
+
+use crate::poly_iop::errors::PolyIOPErrors;
+use arithmetic::identity_permutation_mle;
+use ark_ff::PrimeField;
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer};
+use std::rc::Rc;
+
+/// Returns the evaluations of two MLEs:
+/// - numerator
+/// - denominator
+///
+///  where
+///  - beta and gamma are challenges
+///  - f(x), g(x), s_id(x), s_perm(x) are mle-s
+///
+/// - numerator is the MLE for `f(x) + \beta s_id(x) + \gamma`
+/// - denominator is the MLE for `g(x) + \beta s_perm(x) + \gamma`
+#[allow(clippy::type_complexity)]
+pub(super) fn computer_num_and_denom<F: PrimeField>(
+    beta: &F,
+    gamma: &F,
+    fx: &DenseMultilinearExtension<F>,
+    gx: &DenseMultilinearExtension<F>,
+    s_perm: &DenseMultilinearExtension<F>,
+) -> Result<
+    (
+        Rc<DenseMultilinearExtension<F>>,
+        Rc<DenseMultilinearExtension<F>>,
+    ),
+    PolyIOPErrors,
+> {
+    let start = start_timer!(|| "compute numerator and denominator");
+
+    let num_vars = fx.num_vars;
+    let mut numerator_evals = vec![];
+    let mut denominator_evals = vec![];
+    let s_id = identity_permutation_mle::<F>(num_vars);
+
+    for (&fi, (&gi, (&s_id_i, &s_perm_i))) in
+        fx.iter().zip(gx.iter().zip(s_id.iter().zip(s_perm.iter())))
+    {
+        let numerator = fi + *beta * s_id_i + gamma;
+        let denominator = gi + *beta * s_perm_i + gamma;
+
+        numerator_evals.push(numerator);
+        denominator_evals.push(denominator);
+    }
+    let numerator = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        num_vars,
+        numerator_evals,
+    ));
+    let denominator = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        num_vars,
+        denominator_evals,
+    ));
+
+    end_timer!(start);
+    Ok((numerator, denominator))
+}

subroutines/src/poly_iop/prelude.rs

@@ -0,0 +1,4 @@
+pub use crate::poly_iop::{
+    errors::PolyIOPErrors, perm_check::PermutationCheck, prod_check::ProductCheck,
+    structs::IOPProof, sum_check::SumCheck, utils::*, zero_check::ZeroCheck, PolyIOP,
+};

subroutines/src/poly_iop/prod_check/mod.rs

@@ -0,0 +1,304 @@
+//! Main module for the Product Check protocol
+
+use crate::{
+    pcs::PolynomialCommitmentScheme,
+    poly_iop::{
+        errors::PolyIOPErrors,
+        prod_check::util::{compute_product_poly, prove_zero_check},
+        zero_check::ZeroCheck,
+        PolyIOP,
+    },
+};
+use arithmetic::VPAuxInfo;
+use ark_ec::PairingEngine;
+use ark_ff::{One, PrimeField, Zero};
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer};
+use std::rc::Rc;
+use transcript::IOPTranscript;
+
+mod util;
+
+/// A product-check proves that two n-variate multilinear polynomials `f(x),
+/// g(x)` satisfy:
+/// \prod_{x \in {0,1}^n} f(x) = \prod_{x \in {0,1}^n} g(x)
+///
+/// A ProductCheck is derived from ZeroCheck.
+///
+/// Prover steps:
+/// 1. build `prod(x0, ..., x_n)` from f and g,
+///    such that `prod(0, x1, ..., xn)` equals `f/g` over domain {0,1}^n
+/// 2. push commitments of `prod(x)` to the transcript,
+///    and `generate_challenge` from current transcript (generate alpha)
+/// 3. generate the zerocheck proof for the virtual polynomial
+///    prod(1, x) - prod(x, 0) * prod(x, 1) + alpha * (f(x) - prod(0, x) * g(x))
+///
+/// Verifier steps:
+/// 1. Extract commitments of `prod(x)` from the proof, push
+/// them to the transcript
+/// 2. `generate_challenge` from current transcript (generate alpha)
+/// 3. `verify` to verify the zerocheck proof and generate the subclaim for
+/// polynomial evaluations
+pub trait ProductCheck<E, PCS>: ZeroCheck<E::Fr>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<E>,
+{
+    type ProductCheckSubClaim;
+    type ProductCheckProof;
+
+    /// Initialize the system with a transcript
+    ///
+    /// This function is optional -- in the case where a ProductCheck is
+    /// an building block for a more complex protocol, the transcript
+    /// may be initialized by this complex protocol, and passed to the
+    /// ProductCheck prover/verifier.
+    fn init_transcript() -> Self::Transcript;
+
+    /// Generate a proof for product check, showing that witness multilinear
+    /// polynomials f(x), g(x) satisfy `\prod_{x \in {0,1}^n} f(x) =
+    /// \prod_{x \in {0,1}^n} g(x)`
+    ///
+    /// Inputs:
+    /// - fx: the numerator multilinear polynomial
+    /// - gx: the denominator multilinear polynomial
+    /// - transcript: the IOP transcript
+    /// - pk: PCS committing key
+    ///
+    /// Outputs
+    /// - the product check proof
+    /// - the product polynomial (used for testing)
+    ///
+    /// Cost: O(N)
+    fn prove(
+        pcs_param: &PCS::ProverParam,
+        fx: &Self::MultilinearExtension,
+        gx: &Self::MultilinearExtension,
+        transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<(Self::ProductCheckProof, Self::MultilinearExtension), PolyIOPErrors>;
+
+    /// Verify that for witness multilinear polynomials f(x), g(x)
+    /// it holds that `\prod_{x \in {0,1}^n} f(x) = \prod_{x \in {0,1}^n} g(x)`
+    fn verify(
+        proof: &Self::ProductCheckProof,
+        aux_info: &VPAuxInfo<E::Fr>,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::ProductCheckSubClaim, PolyIOPErrors>;
+}
+
+/// A product check subclaim consists of
+/// - A zero check IOP subclaim for
+/// `Q(x) = prod(1, x) - prod(x, 0) * prod(x, 1) + alpha * (f(x) - prod(0,
+/// x) * g(x)) = 0`
+/// - The random challenge `alpha`
+/// - A final query for `prod(1, ..., 1, 0) = 1`.
+// Note that this final query is in fact a constant that
+// is independent from the proof. So we should avoid
+// (de)serialize it.
+#[derive(Clone, Debug, Default, PartialEq)]
+pub struct ProductCheckSubClaim<F: PrimeField, ZC: ZeroCheck<F>> {
+    // the SubClaim from the ZeroCheck
+    pub zero_check_sub_claim: ZC::ZeroCheckSubClaim,
+    // final query which consists of
+    // - the vector `(1, ..., 1, 0)` (needs to be reversed because Arkwork's MLE uses big-endian
+    //   format for points)
+    // The expected final query evaluation is 1
+    pub final_query: (Vec<F>, F),
+    pub alpha: F,
+}
+
+/// A product check proof consists of
+/// - a zerocheck proof
+/// - a product polynomial commitment
+#[derive(Clone, Debug, Default, PartialEq)]
+pub struct ProductCheckProof<
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<E>,
+    ZC: ZeroCheck<E::Fr>,
+> {
+    pub zero_check_proof: ZC::ZeroCheckProof,
+    pub prod_x_comm: PCS::Commitment,
+}
+
+impl<E, PCS> ProductCheck<E, PCS> for PolyIOP<E::Fr>
+where
+    E: PairingEngine,
+    PCS: PolynomialCommitmentScheme<E, Polynomial = Rc<DenseMultilinearExtension<E::Fr>>>,
+{
+    type ProductCheckSubClaim = ProductCheckSubClaim<E::Fr, Self>;
+    type ProductCheckProof = ProductCheckProof<E, PCS, Self>;
+
+    fn init_transcript() -> Self::Transcript {
+        IOPTranscript::<E::Fr>::new(b"Initializing ProductCheck transcript")
+    }
+
+    fn prove(
+        pcs_param: &PCS::ProverParam,
+        fx: &Self::MultilinearExtension,
+        gx: &Self::MultilinearExtension,
+        transcript: &mut IOPTranscript<E::Fr>,
+    ) -> Result<(Self::ProductCheckProof, Self::MultilinearExtension), PolyIOPErrors> {
+        let start = start_timer!(|| "prod_check prove");
+
+        if fx.num_vars != gx.num_vars {
+            return Err(PolyIOPErrors::InvalidParameters(
+                "fx and gx have different number of variables".to_string(),
+            ));
+        }
+
+        // compute the product polynomial
+        let prod_x = compute_product_poly(fx, gx)?;
+
+        // generate challenge
+        let prod_x_comm = PCS::commit(pcs_param, &prod_x)?;
+        transcript.append_serializable_element(b"prod(x)", &prod_x_comm)?;
+        let alpha = transcript.get_and_append_challenge(b"alpha")?;
+
+        // build the zero-check proof
+        let (zero_check_proof, _) = prove_zero_check(fx, gx, &prod_x, &alpha, transcript)?;
+
+        end_timer!(start);
+
+        Ok((
+            ProductCheckProof {
+                zero_check_proof,
+                prod_x_comm,
+            },
+            prod_x,
+        ))
+    }
+
+    fn verify(
+        proof: &Self::ProductCheckProof,
+        aux_info: &VPAuxInfo<E::Fr>,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::ProductCheckSubClaim, PolyIOPErrors> {
+        let start = start_timer!(|| "prod_check verify");
+
+        // update transcript and generate challenge
+        transcript.append_serializable_element(b"prod(x)", &proof.prod_x_comm)?;
+        let alpha = transcript.get_and_append_challenge(b"alpha")?;
+
+        // invoke the zero check on the iop_proof
+        // the virtual poly info for Q(x)
+        let zero_check_sub_claim =
+            <Self as ZeroCheck<E::Fr>>::verify(&proof.zero_check_proof, aux_info, transcript)?;
+
+        // the final query is on prod_x, hence has length `num_vars` + 1
+        let mut final_query = vec![E::Fr::one(); aux_info.num_variables + 1];
+        // the point has to be reversed because Arkworks uses big-endian.
+        final_query[0] = E::Fr::zero();
+        let final_eval = E::Fr::one();
+
+        end_timer!(start);
+
+        Ok(ProductCheckSubClaim {
+            zero_check_sub_claim,
+            final_query: (final_query, final_eval),
+            alpha,
+        })
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::ProductCheck;
+    use crate::{
+        pcs::{prelude::MultilinearKzgPCS, PolynomialCommitmentScheme},
+        poly_iop::{errors::PolyIOPErrors, PolyIOP},
+    };
+    use arithmetic::VPAuxInfo;
+    use ark_bls12_381::{Bls12_381, Fr};
+    use ark_ec::PairingEngine;
+    use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+    use ark_std::test_rng;
+    use std::{marker::PhantomData, rc::Rc};
+
+    // f and g are guaranteed to have the same product
+    fn test_product_check_helper<E, PCS>(
+        f: &DenseMultilinearExtension<E::Fr>,
+        g: &DenseMultilinearExtension<E::Fr>,
+        pcs_param: &PCS::ProverParam,
+    ) -> Result<(), PolyIOPErrors>
+    where
+        E: PairingEngine,
+        PCS: PolynomialCommitmentScheme<E, Polynomial = Rc<DenseMultilinearExtension<E::Fr>>>,
+    {
+        let mut transcript = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::init_transcript();
+        transcript.append_message(b"testing", b"initializing transcript for testing")?;
+
+        let (proof, prod_x) = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::prove(
+            pcs_param,
+            &Rc::new(f.clone()),
+            &Rc::new(g.clone()),
+            &mut transcript,
+        )?;
+
+        let mut transcript = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::init_transcript();
+        transcript.append_message(b"testing", b"initializing transcript for testing")?;
+
+        let aux_info = VPAuxInfo {
+            max_degree: 2,
+            num_variables: f.num_vars,
+            phantom: PhantomData::default(),
+        };
+        let prod_subclaim =
+            <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::verify(&proof, &aux_info, &mut transcript)?;
+        assert_eq!(
+            prod_x.evaluate(&prod_subclaim.final_query.0).unwrap(),
+            prod_subclaim.final_query.1,
+            "different product"
+        );
+
+        // bad path
+        let mut transcript = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::init_transcript();
+        transcript.append_message(b"testing", b"initializing transcript for testing")?;
+
+        let h = f + g;
+        let (bad_proof, prod_x_bad) = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::prove(
+            pcs_param,
+            &Rc::new(f.clone()),
+            &Rc::new(h),
+            &mut transcript,
+        )?;
+
+        let mut transcript = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::init_transcript();
+        transcript.append_message(b"testing", b"initializing transcript for testing")?;
+        let bad_subclaim = <PolyIOP<E::Fr> as ProductCheck<E, PCS>>::verify(
+            &bad_proof,
+            &aux_info,
+            &mut transcript,
+        )?;
+        assert_ne!(
+            prod_x_bad.evaluate(&bad_subclaim.final_query.0).unwrap(),
+            bad_subclaim.final_query.1,
+            "can't detect wrong proof"
+        );
+
+        Ok(())
+    }
+
+    fn test_product_check(nv: usize) -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+
+        let f: DenseMultilinearExtension<Fr> = DenseMultilinearExtension::rand(nv, &mut rng);
+        let mut g = f.clone();
+        g.evaluations.reverse();
+
+        let srs = MultilinearKzgPCS::<Bls12_381>::gen_srs_for_testing(&mut rng, nv + 1)?;
+        let (pcs_param, _) = MultilinearKzgPCS::<Bls12_381>::trim(&srs, None, Some(nv + 1))?;
+
+        test_product_check_helper::<Bls12_381, MultilinearKzgPCS<Bls12_381>>(&f, &g, &pcs_param)?;
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_trivial_polynomial() -> Result<(), PolyIOPErrors> {
+        test_product_check(1)
+    }
+    #[test]
+    fn test_normal_polynomial() -> Result<(), PolyIOPErrors> {
+        test_product_check(10)
+    }
+}

subroutines/src/poly_iop/prod_check/util.rs

@@ -0,0 +1,202 @@
+//! This module implements useful functions for the product check protocol.
+
+use crate::poly_iop::{errors::PolyIOPErrors, structs::IOPProof, zero_check::ZeroCheck, PolyIOP};
+use arithmetic::{get_index, VirtualPolynomial};
+use ark_ff::PrimeField;
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer};
+use rayon::prelude::{IntoParallelRefIterator, ParallelIterator};
+use std::rc::Rc;
+use transcript::IOPTranscript;
+
+/// Compute the product polynomial `prod(x)` where
+///
+///  - `prod(0,x) := prod(0, x1, …, xn)` is the MLE over the
+/// evaluations of `f(x)/g(x)` on the boolean hypercube {0,1}^n
+///
+/// - `prod(1,x)` is a MLE over the evaluations of `prod(x, 0) * prod(x, 1)`
+/// on the boolean hypercube {0,1}^n
+///
+/// The caller needs to check num_vars matches in f and g
+/// Cost: linear in N.
+pub(super) fn compute_product_poly<F: PrimeField>(
+    fx: &Rc<DenseMultilinearExtension<F>>,
+    gx: &Rc<DenseMultilinearExtension<F>>,
+) -> Result<Rc<DenseMultilinearExtension<F>>, PolyIOPErrors> {
+    let start = start_timer!(|| "compute evaluations of prod polynomial");
+    let num_vars = fx.num_vars;
+
+    // ===================================
+    // prod(0, x)
+    // ===================================
+    let prod_0x_eval = compute_prod_0(fx, gx)?;
+
+    // ===================================
+    // prod(1, x)
+    // ===================================
+    //
+    // `prod(1, x)` can be computed via recursing the following formula for 2^n-1
+    // times
+    //
+    // `prod(1, x_1, ..., x_n) :=
+    //      prod(x_1, x_2, ..., x_n, 0) * prod(x_1, x_2, ..., x_n, 1)`
+    //
+    // At any given step, the right hand side of the equation
+    // is available via either eval_0x or the current view of eval_1x
+    let mut prod_1x_eval = vec![];
+    for x in 0..(1 << num_vars) - 1 {
+        // sign will decide if the evaluation should be looked up from eval_0x or
+        // eval_1x; x_zero_index is the index for the evaluation (x_2, ..., x_n,
+        // 0); x_one_index is the index for the evaluation (x_2, ..., x_n, 1);
+        let (x_zero_index, x_one_index, sign) = get_index(x, num_vars);
+        if !sign {
+            prod_1x_eval.push(prod_0x_eval[x_zero_index] * prod_0x_eval[x_one_index]);
+        } else {
+            // sanity check: if we are trying to look up from the eval_1x table,
+            // then the target index must already exist
+            if x_zero_index >= prod_1x_eval.len() || x_one_index >= prod_1x_eval.len() {
+                return Err(PolyIOPErrors::ShouldNotArrive);
+            }
+            prod_1x_eval.push(prod_1x_eval[x_zero_index] * prod_1x_eval[x_one_index]);
+        }
+    }
+
+    // prod(1, 1, ..., 1) := 0
+    prod_1x_eval.push(F::zero());
+
+    // ===================================
+    // prod(x)
+    // ===================================
+    // prod(x)'s evaluation is indeed `e := [eval_0x[..], eval_1x[..]].concat()`
+    let eval = [prod_0x_eval.as_slice(), prod_1x_eval.as_slice()].concat();
+
+    let prod_x = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        num_vars + 1,
+        eval,
+    ));
+
+    end_timer!(start);
+    Ok(prod_x)
+}
+
+/// generate the zerocheck proof for the virtual polynomial
+///    prod(1, x) - prod(x, 0) * prod(x, 1) + alpha * (f(x) - prod(0, x) * g(x))
+///
+/// Returns proof and Q(x) for testing purpose.
+///
+/// Cost: O(N)
+pub(super) fn prove_zero_check<F: PrimeField>(
+    fx: &Rc<DenseMultilinearExtension<F>>,
+    gx: &Rc<DenseMultilinearExtension<F>>,
+    prod_x: &Rc<DenseMultilinearExtension<F>>,
+    alpha: &F,
+    transcript: &mut IOPTranscript<F>,
+) -> Result<(IOPProof<F>, VirtualPolynomial<F>), PolyIOPErrors> {
+    let start = start_timer!(|| "zerocheck in product check");
+
+    let prod_partial_evals = build_prod_partial_eval(prod_x)?;
+    let prod_0x = prod_partial_evals[0].clone();
+    let prod_1x = prod_partial_evals[1].clone();
+    let prod_x0 = prod_partial_evals[2].clone();
+    let prod_x1 = prod_partial_evals[3].clone();
+
+    // compute g(x) * prod(0, x) * alpha
+    let mut q_x = VirtualPolynomial::new_from_mle(gx, F::one());
+    q_x.mul_by_mle(prod_0x, *alpha)?;
+
+    //   g(x) * prod(0, x) * alpha
+    // - f(x) * alpha
+    q_x.add_mle_list([fx.clone()], -*alpha)?;
+
+    // Q(x) := prod(1,x) - prod(x, 0) * prod(x, 1)
+    //       + alpha * (
+    //             g(x) * prod(0, x)
+    //           - f(x))
+    q_x.add_mle_list([prod_x0, prod_x1], -F::one())?;
+    q_x.add_mle_list([prod_1x], F::one())?;
+
+    let iop_proof = <PolyIOP<F> as ZeroCheck<F>>::prove(&q_x, transcript)?;
+
+    end_timer!(start);
+    Ok((iop_proof, q_x))
+}
+
+/// Helper function of the IOP.
+///
+/// Input:
+/// - prod(x)
+///
+/// Output: the following 4 polynomials
+/// - prod(0, x)
+/// - prod(1, x)
+/// - prod(x, 0)
+/// - prod(x, 1)
+fn build_prod_partial_eval<F: PrimeField>(
+    prod_x: &Rc<DenseMultilinearExtension<F>>,
+) -> Result<[Rc<DenseMultilinearExtension<F>>; 4], PolyIOPErrors> {
+    let start = start_timer!(|| "build partial prod polynomial");
+
+    let prod_x_eval = &prod_x.evaluations;
+    let num_vars = prod_x.num_vars - 1;
+
+    // prod(0, x)
+    let prod_0_x = Rc::new(DenseMultilinearExtension::from_evaluations_slice(
+        num_vars,
+        &prod_x_eval[0..1 << num_vars],
+    ));
+    // prod(1, x)
+    let prod_1_x = Rc::new(DenseMultilinearExtension::from_evaluations_slice(
+        num_vars,
+        &prod_x_eval[1 << num_vars..1 << (num_vars + 1)],
+    ));
+
+    // ===================================
+    // prod(x, 0) and prod(x, 1)
+    // ===================================
+    //
+    // now we compute eval_x0 and eval_x1
+    // eval_0x will be the odd coefficients of eval
+    // and eval_1x will be the even coefficients of eval
+    let mut eval_x0 = vec![];
+    let mut eval_x1 = vec![];
+    for (x, &prod_x) in prod_x_eval.iter().enumerate() {
+        if x & 1 == 0 {
+            eval_x0.push(prod_x);
+        } else {
+            eval_x1.push(prod_x);
+        }
+    }
+    let prod_x_0 = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        num_vars, eval_x0,
+    ));
+    let prod_x_1 = Rc::new(DenseMultilinearExtension::from_evaluations_vec(
+        num_vars, eval_x1,
+    ));
+
+    end_timer!(start);
+
+    Ok([prod_0_x, prod_1_x, prod_x_0, prod_x_1])
+}
+
+/// Returns the evaluations of
+/// - `prod(0,x) := prod(0, x1, …, xn)` which is the MLE over the
+/// evaluations of f(x)/g(x) on the boolean hypercube {0,1}^n:
+///
+/// The caller needs to check num_vars matches in f/g
+/// Cost: linear in N.
+fn compute_prod_0<F: PrimeField>(
+    fx: &DenseMultilinearExtension<F>,
+    gx: &DenseMultilinearExtension<F>,
+) -> Result<Vec<F>, PolyIOPErrors> {
+    let start = start_timer!(|| "compute prod(0,x)");
+
+    let input = fx
+        .iter()
+        .zip(gx.iter())
+        .map(|(&fi, &gi)| (fi, gi))
+        .collect::<Vec<_>>();
+    let prod_0x_evals = input.par_iter().map(|(x, y)| *x / *y).collect::<Vec<_>>();
+
+    end_timer!(start);
+    Ok(prod_0x_evals)
+}

subroutines/src/poly_iop/readme.md

@@ -0,0 +1,9 @@
+Poly IOP
+-----
+
+Implements the following protocols
+
+- sum checks
+- zero checks
+- product checks
+- permutation checks

subroutines/src/poly_iop/structs.rs

@@ -0,0 +1,45 @@
+//! This module defines structs that are shared by all sub protocols.
+
+use arithmetic::VirtualPolynomial;
+use ark_ff::PrimeField;
+use ark_serialize::{CanonicalSerialize, SerializationError, Write};
+
+/// An IOP proof is a collections of
+/// - messages from prover to verifier at each round through the interactive
+///   protocol.
+/// - a point that is generated by the transcript for evaluation
+#[derive(Clone, Debug, Default, PartialEq, Eq)]
+pub struct IOPProof<F: PrimeField> {
+    pub point: Vec<F>,
+    pub proofs: Vec<IOPProverMessage<F>>,
+}
+
+/// A message from the prover to the verifier at a given round
+/// is a list of evaluations.
+#[derive(Clone, Debug, Default, PartialEq, Eq, CanonicalSerialize)]
+pub struct IOPProverMessage<F: PrimeField> {
+    pub(crate) evaluations: Vec<F>,
+}
+
+/// Prover State of a PolyIOP.
+pub struct IOPProverState<F: PrimeField> {
+    /// sampled randomness given by the verifier
+    pub challenges: Vec<F>,
+    /// the current round number
+    pub(crate) round: usize,
+    /// pointer to the virtual polynomial
+    pub(crate) poly: VirtualPolynomial<F>,
+}
+
+/// Prover State of a PolyIOP
+pub struct IOPVerifierState<F: PrimeField> {
+    pub(crate) round: usize,
+    pub(crate) num_vars: usize,
+    pub(crate) max_degree: usize,
+    pub(crate) finished: bool,
+    /// a list storing the univariate polynomial in evaluation form sent by the
+    /// prover at each round
+    pub(crate) polynomials_received: Vec<Vec<F>>,
+    /// a list storing the randomness sampled by the verifier at each round
+    pub(crate) challenges: Vec<F>,
+}

subroutines/src/poly_iop/sum_check/mod.rs

@@ -0,0 +1,387 @@
+//! This module implements the sum check protocol.
+
+use crate::poly_iop::{
+    errors::PolyIOPErrors,
+    structs::{IOPProof, IOPProverState, IOPVerifierState},
+    PolyIOP,
+};
+use arithmetic::{VPAuxInfo, VirtualPolynomial};
+use ark_ff::PrimeField;
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer};
+use std::{fmt::Debug, rc::Rc};
+use transcript::IOPTranscript;
+
+mod prover;
+mod verifier;
+
+/// Trait for doing sum check protocols.
+pub trait SumCheck<F: PrimeField> {
+    type VirtualPolynomial;
+    type VPAuxInfo;
+    type MultilinearExtension;
+
+    type SumCheckProof: Clone + Debug + Default + PartialEq;
+    type Transcript;
+    type SumCheckSubClaim: Clone + Debug + Default + PartialEq;
+
+    /// Extract sum from the proof
+    fn extract_sum(proof: &Self::SumCheckProof) -> F;
+
+    /// Initialize the system with a transcript
+    ///
+    /// This function is optional -- in the case where a SumCheck is
+    /// an building block for a more complex protocol, the transcript
+    /// may be initialized by this complex protocol, and passed to the
+    /// SumCheck prover/verifier.
+    fn init_transcript() -> Self::Transcript;
+
+    /// Generate proof of the sum of polynomial over {0,1}^`num_vars`
+    ///
+    /// The polynomial is represented in the form of a VirtualPolynomial.
+    fn prove(
+        poly: &Self::VirtualPolynomial,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::SumCheckProof, PolyIOPErrors>;
+
+    /// Verify the claimed sum using the proof
+    fn verify(
+        sum: F,
+        proof: &Self::SumCheckProof,
+        aux_info: &Self::VPAuxInfo,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::SumCheckSubClaim, PolyIOPErrors>;
+}
+
+/// Trait for sum check protocol prover side APIs.
+pub trait SumCheckProver<F: PrimeField>
+where
+    Self: Sized,
+{
+    type VirtualPolynomial;
+    type ProverMessage;
+
+    /// Initialize the prover state to argue for the sum of the input polynomial
+    /// over {0,1}^`num_vars`.
+    fn prover_init(polynomial: &Self::VirtualPolynomial) -> Result<Self, PolyIOPErrors>;
+
+    /// Receive message from verifier, generate prover message, and proceed to
+    /// next round.
+    ///
+    /// Main algorithm used is from section 3.2 of [XZZPS19](https://eprint.iacr.org/2019/317.pdf#subsection.3.2).
+    fn prove_round_and_update_state(
+        &mut self,
+        challenge: &Option<F>,
+    ) -> Result<Self::ProverMessage, PolyIOPErrors>;
+}
+
+/// Trait for sum check protocol verifier side APIs.
+pub trait SumCheckVerifier<F: PrimeField> {
+    type VPAuxInfo;
+    type ProverMessage;
+    type Challenge;
+    type Transcript;
+    type SumCheckSubClaim;
+
+    /// Initialize the verifier's state.
+    fn verifier_init(index_info: &Self::VPAuxInfo) -> Self;
+
+    /// Run verifier for the current round, given a prover message.
+    ///
+    /// Note that `verify_round_and_update_state` only samples and stores
+    /// challenges; and update the verifier's state accordingly. The actual
+    /// verifications are deferred (in batch) to `check_and_generate_subclaim`
+    /// at the last step.
+    fn verify_round_and_update_state(
+        &mut self,
+        prover_msg: &Self::ProverMessage,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::Challenge, PolyIOPErrors>;
+
+    /// This function verifies the deferred checks in the interactive version of
+    /// the protocol; and generate the subclaim. Returns an error if the
+    /// proof failed to verify.
+    ///
+    /// If the asserted sum is correct, then the multilinear polynomial
+    /// evaluated at `subclaim.point` will be `subclaim.expected_evaluation`.
+    /// Otherwise, it is highly unlikely that those two will be equal.
+    /// Larger field size guarantees smaller soundness error.
+    fn check_and_generate_subclaim(
+        &self,
+        asserted_sum: &F,
+    ) -> Result<Self::SumCheckSubClaim, PolyIOPErrors>;
+}
+
+/// A SumCheckSubClaim is a claim generated by the verifier at the end of
+/// verification when it is convinced.
+#[derive(Clone, Debug, Default, PartialEq, Eq)]
+pub struct SumCheckSubClaim<F: PrimeField> {
+    /// the multi-dimensional point that this multilinear extension is evaluated
+    /// to
+    pub point: Vec<F>,
+    /// the expected evaluation
+    pub expected_evaluation: F,
+}
+
+impl<F: PrimeField> SumCheck<F> for PolyIOP<F> {
+    type SumCheckProof = IOPProof<F>;
+    type VirtualPolynomial = VirtualPolynomial<F>;
+    type VPAuxInfo = VPAuxInfo<F>;
+    type MultilinearExtension = Rc<DenseMultilinearExtension<F>>;
+    type SumCheckSubClaim = SumCheckSubClaim<F>;
+    type Transcript = IOPTranscript<F>;
+
+    fn extract_sum(proof: &Self::SumCheckProof) -> F {
+        let start = start_timer!(|| "extract sum");
+        let res = proof.proofs[0].evaluations[0] + proof.proofs[0].evaluations[1];
+        end_timer!(start);
+        res
+    }
+
+    fn init_transcript() -> Self::Transcript {
+        let start = start_timer!(|| "init transcript");
+        let res = IOPTranscript::<F>::new(b"Initializing SumCheck transcript");
+        end_timer!(start);
+        res
+    }
+
+    fn prove(
+        poly: &Self::VirtualPolynomial,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::SumCheckProof, PolyIOPErrors> {
+        let start = start_timer!(|| "sum check prove");
+
+        transcript.append_serializable_element(b"aux info", &poly.aux_info)?;
+
+        let mut prover_state = IOPProverState::prover_init(poly)?;
+        let mut challenge = None;
+        let mut prover_msgs = Vec::with_capacity(poly.aux_info.num_variables);
+        for _ in 0..poly.aux_info.num_variables {
+            let prover_msg =
+                IOPProverState::prove_round_and_update_state(&mut prover_state, &challenge)?;
+            transcript.append_serializable_element(b"prover msg", &prover_msg)?;
+            prover_msgs.push(prover_msg);
+            challenge = Some(transcript.get_and_append_challenge(b"Internal round")?);
+        }
+        // pushing the last challenge point to the state
+        if let Some(p) = challenge {
+            prover_state.challenges.push(p)
+        };
+
+        end_timer!(start);
+        Ok(IOPProof {
+            point: prover_state.challenges,
+            proofs: prover_msgs,
+        })
+    }
+
+    fn verify(
+        claimed_sum: F,
+        proof: &Self::SumCheckProof,
+        aux_info: &Self::VPAuxInfo,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::SumCheckSubClaim, PolyIOPErrors> {
+        let start = start_timer!(|| "sum check verify");
+
+        transcript.append_serializable_element(b"aux info", aux_info)?;
+        let mut verifier_state = IOPVerifierState::verifier_init(aux_info);
+        for i in 0..aux_info.num_variables {
+            let prover_msg = proof.proofs.get(i).expect("proof is incomplete");
+            transcript.append_serializable_element(b"prover msg", prover_msg)?;
+            IOPVerifierState::verify_round_and_update_state(
+                &mut verifier_state,
+                prover_msg,
+                transcript,
+            )?;
+        }
+
+        let res = IOPVerifierState::check_and_generate_subclaim(&verifier_state, &claimed_sum);
+
+        end_timer!(start);
+        res
+    }
+}
+
+#[cfg(test)]
+mod test {
+
+    use super::*;
+    use ark_bls12_381::Fr;
+    use ark_ff::UniformRand;
+    use ark_poly::{DenseMultilinearExtension, MultilinearExtension};
+    use ark_std::test_rng;
+    use std::rc::Rc;
+
+    fn test_sumcheck(
+        nv: usize,
+        num_multiplicands_range: (usize, usize),
+        num_products: usize,
+    ) -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+        let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+
+        let (poly, asserted_sum) =
+            VirtualPolynomial::rand(nv, num_multiplicands_range, num_products, &mut rng)?;
+        let proof = <PolyIOP<Fr> as SumCheck<Fr>>::prove(&poly, &mut transcript)?;
+        let poly_info = poly.aux_info.clone();
+        let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+        let subclaim = <PolyIOP<Fr> as SumCheck<Fr>>::verify(
+            asserted_sum,
+            &proof,
+            &poly_info,
+            &mut transcript,
+        )?;
+        assert!(
+            poly.evaluate(&subclaim.point).unwrap() == subclaim.expected_evaluation,
+            "wrong subclaim"
+        );
+        Ok(())
+    }
+
+    fn test_sumcheck_internal(
+        nv: usize,
+        num_multiplicands_range: (usize, usize),
+        num_products: usize,
+    ) -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+        let (poly, asserted_sum) =
+            VirtualPolynomial::<Fr>::rand(nv, num_multiplicands_range, num_products, &mut rng)?;
+        let poly_info = poly.aux_info.clone();
+        let mut prover_state = IOPProverState::prover_init(&poly)?;
+        let mut verifier_state = IOPVerifierState::verifier_init(&poly_info);
+        let mut challenge = None;
+        let mut transcript = IOPTranscript::new(b"a test transcript");
+        transcript
+            .append_message(b"testing", b"initializing transcript for testing")
+            .unwrap();
+        for _ in 0..poly.aux_info.num_variables {
+            let prover_message =
+                IOPProverState::prove_round_and_update_state(&mut prover_state, &challenge)
+                    .unwrap();
+
+            challenge = Some(
+                IOPVerifierState::verify_round_and_update_state(
+                    &mut verifier_state,
+                    &prover_message,
+                    &mut transcript,
+                )
+                .unwrap(),
+            );
+        }
+        let subclaim =
+            IOPVerifierState::check_and_generate_subclaim(&verifier_state, &asserted_sum)
+                .expect("fail to generate subclaim");
+        assert!(
+            poly.evaluate(&subclaim.point).unwrap() == subclaim.expected_evaluation,
+            "wrong subclaim"
+        );
+        Ok(())
+    }
+
+    #[test]
+    fn test_trivial_polynomial() -> Result<(), PolyIOPErrors> {
+        let nv = 1;
+        let num_multiplicands_range = (4, 13);
+        let num_products = 5;
+
+        test_sumcheck(nv, num_multiplicands_range, num_products)?;
+        test_sumcheck_internal(nv, num_multiplicands_range, num_products)
+    }
+    #[test]
+    fn test_normal_polynomial() -> Result<(), PolyIOPErrors> {
+        let nv = 12;
+        let num_multiplicands_range = (4, 9);
+        let num_products = 5;
+
+        test_sumcheck(nv, num_multiplicands_range, num_products)?;
+        test_sumcheck_internal(nv, num_multiplicands_range, num_products)
+    }
+    #[test]
+    fn zero_polynomial_should_error() {
+        let nv = 0;
+        let num_multiplicands_range = (4, 13);
+        let num_products = 5;
+
+        assert!(test_sumcheck(nv, num_multiplicands_range, num_products).is_err());
+        assert!(test_sumcheck_internal(nv, num_multiplicands_range, num_products).is_err());
+    }
+
+    #[test]
+    fn test_extract_sum() -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+        let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+        let (poly, asserted_sum) = VirtualPolynomial::<Fr>::rand(8, (3, 4), 3, &mut rng)?;
+
+        let proof = <PolyIOP<Fr> as SumCheck<Fr>>::prove(&poly, &mut transcript)?;
+        assert_eq!(
+            <PolyIOP<Fr> as SumCheck<Fr>>::extract_sum(&proof),
+            asserted_sum
+        );
+        Ok(())
+    }
+
+    #[test]
+    /// Test that the memory usage of shared-reference is linear to number of
+    /// unique MLExtensions instead of total number of multiplicands.
+    fn test_shared_reference() -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+        let ml_extensions: Vec<_> = (0..5)
+            .map(|_| Rc::new(DenseMultilinearExtension::<Fr>::rand(8, &mut rng)))
+            .collect();
+        let mut poly = VirtualPolynomial::new(8);
+        poly.add_mle_list(
+            vec![
+                ml_extensions[2].clone(),
+                ml_extensions[3].clone(),
+                ml_extensions[0].clone(),
+            ],
+            Fr::rand(&mut rng),
+        )?;
+        poly.add_mle_list(
+            vec![
+                ml_extensions[1].clone(),
+                ml_extensions[4].clone(),
+                ml_extensions[4].clone(),
+            ],
+            Fr::rand(&mut rng),
+        )?;
+        poly.add_mle_list(
+            vec![
+                ml_extensions[3].clone(),
+                ml_extensions[2].clone(),
+                ml_extensions[1].clone(),
+            ],
+            Fr::rand(&mut rng),
+        )?;
+        poly.add_mle_list(
+            vec![ml_extensions[0].clone(), ml_extensions[0].clone()],
+            Fr::rand(&mut rng),
+        )?;
+        poly.add_mle_list(vec![ml_extensions[4].clone()], Fr::rand(&mut rng))?;
+
+        assert_eq!(poly.flattened_ml_extensions.len(), 5);
+
+        // test memory usage for prover
+        let prover = IOPProverState::<Fr>::prover_init(&poly).unwrap();
+        assert_eq!(prover.poly.flattened_ml_extensions.len(), 5);
+        drop(prover);
+
+        let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+        let poly_info = poly.aux_info.clone();
+        let proof = <PolyIOP<Fr> as SumCheck<Fr>>::prove(&poly, &mut transcript)?;
+        let asserted_sum = <PolyIOP<Fr> as SumCheck<Fr>>::extract_sum(&proof);
+
+        let mut transcript = <PolyIOP<Fr> as SumCheck<Fr>>::init_transcript();
+        let subclaim = <PolyIOP<Fr> as SumCheck<Fr>>::verify(
+            asserted_sum,
+            &proof,
+            &poly_info,
+            &mut transcript,
+        )?;
+        assert!(
+            poly.evaluate(&subclaim.point)? == subclaim.expected_evaluation,
+            "wrong subclaim"
+        );
+        Ok(())
+    }
+}

subroutines/src/poly_iop/sum_check/prover.rs

@@ -0,0 +1,200 @@
+//! Prover subroutines for a SumCheck protocol.
+
+use super::SumCheckProver;
+use crate::poly_iop::{
+    errors::PolyIOPErrors,
+    structs::{IOPProverMessage, IOPProverState},
+};
+use arithmetic::{fix_variables, VirtualPolynomial};
+use ark_ff::PrimeField;
+use ark_poly::DenseMultilinearExtension;
+use ark_std::{end_timer, start_timer, vec::Vec};
+use rayon::prelude::IntoParallelIterator;
+use std::rc::Rc;
+
+#[cfg(feature = "parallel")]
+use rayon::iter::{IntoParallelRefMutIterator, ParallelIterator};
+
+impl<F: PrimeField> SumCheckProver<F> for IOPProverState<F> {
+    type VirtualPolynomial = VirtualPolynomial<F>;
+    type ProverMessage = IOPProverMessage<F>;
+
+    /// Initialize the prover state to argue for the sum of the input polynomial
+    /// over {0,1}^`num_vars`.
+    fn prover_init(polynomial: &Self::VirtualPolynomial) -> Result<Self, PolyIOPErrors> {
+        let start = start_timer!(|| "sum check prover init");
+        if polynomial.aux_info.num_variables == 0 {
+            return Err(PolyIOPErrors::InvalidParameters(
+                "Attempt to prove a constant.".to_string(),
+            ));
+        }
+        end_timer!(start);
+
+        Ok(Self {
+            challenges: Vec::with_capacity(polynomial.aux_info.num_variables),
+            round: 0,
+            poly: polynomial.clone(),
+        })
+    }
+
+    /// Receive message from verifier, generate prover message, and proceed to
+    /// next round.
+    ///
+    /// Main algorithm used is from section 3.2 of [XZZPS19](https://eprint.iacr.org/2019/317.pdf#subsection.3.2).
+    fn prove_round_and_update_state(
+        &mut self,
+        challenge: &Option<F>,
+    ) -> Result<Self::ProverMessage, PolyIOPErrors> {
+        // let start =
+        //     start_timer!(|| format!("sum check prove {}-th round and update state",
+        // self.round));
+
+        if self.round >= self.poly.aux_info.num_variables {
+            return Err(PolyIOPErrors::InvalidProver(
+                "Prover is not active".to_string(),
+            ));
+        }
+
+        // let fix_argument = start_timer!(|| "fix argument");
+
+        // Step 1:
+        // fix argument and evaluate f(x) over x_m = r; where r is the challenge
+        // for the current round, and m is the round number, indexed from 1
+        //
+        // i.e.:
+        // at round m <= n, for each mle g(x_1, ... x_n) within the flattened_mle
+        // which has already been evaluated to
+        //
+        //    g(r_1, ..., r_{m-1}, x_m ... x_n)
+        //
+        // eval g over r_m, and mutate g to g(r_1, ... r_m,, x_{m+1}... x_n)
+        let mut flattened_ml_extensions: Vec<DenseMultilinearExtension<F>> = self
+            .poly
+            .flattened_ml_extensions
+            .iter()
+            .map(|x| x.as_ref().clone())
+            .collect();
+
+        if let Some(chal) = challenge {
+            if self.round == 0 {
+                return Err(PolyIOPErrors::InvalidProver(
+                    "first round should be prover first.".to_string(),
+                ));
+            }
+            self.challenges.push(*chal);
+
+            let r = self.challenges[self.round - 1];
+            #[cfg(feature = "parallel")]
+            flattened_ml_extensions
+                .par_iter_mut()
+                .for_each(|mle| *mle = fix_variables(mle, &[r]));
+            #[cfg(not(feature = "parallel"))]
+            flattened_ml_extensions
+                .iter_mut()
+                .for_each(|mle| *mle = fix_variables(mle, &[r]));
+        } else if self.round > 0 {
+            return Err(PolyIOPErrors::InvalidProver(
+                "verifier message is empty".to_string(),
+            ));
+        }
+        // end_timer!(fix_argument);
+
+        self.round += 1;
+
+        let products_list = self.poly.products.clone();
+        let mut products_sum = vec![F::zero(); self.poly.aux_info.max_degree + 1];
+
+        // let compute_sum = start_timer!(|| "compute sum");
+
+        // Step 2: generate sum for the partial evaluated polynomial:
+        // f(r_1, ... r_m,, x_{m+1}... x_n)
+
+        #[cfg(feature = "parallel")]
+        {
+            let flag = (self.poly.aux_info.max_degree == 2)
+                && (products_list.len() == 1)
+                && (products_list[0].0 == F::one());
+            if flag {
+                for (t, e) in products_sum.iter_mut().enumerate() {
+                    let evals = (0..1 << (self.poly.aux_info.num_variables - self.round))
+                        .into_par_iter()
+                        .map(|b| {
+                            // evaluate P_round(t)
+                            let table0 = &flattened_ml_extensions[products_list[0].1[0]];
+                            let table1 = &flattened_ml_extensions[products_list[0].1[1]];
+                            if t == 0 {
+                                table0[b << 1] * table1[b << 1]
+                            } else if t == 1 {
+                                table0[(b << 1) + 1] * table1[(b << 1) + 1]
+                            } else {
+                                (table0[(b << 1) + 1] + table0[(b << 1) + 1] - table0[b << 1])
+                                    * (table1[(b << 1) + 1] + table1[(b << 1) + 1] - table1[b << 1])
+                            }
+                        })
+                        .collect::<Vec<F>>();
+                    for val in evals.iter() {
+                        *e += val
+                    }
+                }
+            } else {
+                for (t, e) in products_sum.iter_mut().enumerate() {
+                    let t = F::from(t as u64);
+                    let products = (0..1 << (self.poly.aux_info.num_variables - self.round))
+                        .into_par_iter()
+                        .map(|b| {
+                            // evaluate P_round(t)
+                            let mut tmp = F::zero();
+                            products_list.iter().for_each(|(coefficient, products)| {
+                                let num_mles = products.len();
+                                let mut product = *coefficient;
+                                for &f in products.iter().take(num_mles) {
+                                    let table = &flattened_ml_extensions[f]; // f's range is checked in init
+                                    product *=
+                                        table[b << 1] + (table[(b << 1) + 1] - table[b << 1]) * t;
+                                }
+                                tmp += product;
+                            });
+
+                            tmp
+                        })
+                        .collect::<Vec<F>>();
+
+                    for i in products.iter() {
+                        *e += i
+                    }
+                }
+            }
+        }
+
+        #[cfg(not(feature = "parallel"))]
+        products_sum.iter_mut().enumerate().for_each(|(t, e)| {
+            let t = F::from(t as u64);
+            let one_minus_t = F::one() - t;
+
+            for b in 0..1 << (self.poly.aux_info.num_variables - self.round) {
+                // evaluate P_round(t)
+                for (coefficient, products) in products_list.iter() {
+                    let num_mles = products.len();
+                    let mut product = *coefficient;
+                    for &f in products.iter().take(num_mles) {
+                        let table = &flattened_ml_extensions[f]; // f's range is checked in init
+                        product *= table[b << 1] + (table[(b << 1) + 1] - table[b << 1]) * t;
+                    }
+                    *e += product;
+                }
+            }
+        });
+
+        // update prover's state to the partial evaluated polynomial
+        self.poly.flattened_ml_extensions = flattened_ml_extensions
+            .iter()
+            .map(|x| Rc::new(x.clone()))
+            .collect();
+
+        // end_timer!(compute_sum);
+        // end_timer!(start);
+        Ok(IOPProverMessage {
+            evaluations: products_sum,
+        })
+    }
+}

subroutines/src/poly_iop/sum_check/verifier.rs

@@ -0,0 +1,352 @@
+//! Verifier subroutines for a SumCheck protocol.
+
+use super::{SumCheckSubClaim, SumCheckVerifier};
+use crate::poly_iop::{
+    errors::PolyIOPErrors,
+    structs::{IOPProverMessage, IOPVerifierState},
+};
+use arithmetic::VPAuxInfo;
+use ark_ff::PrimeField;
+use ark_std::{end_timer, start_timer};
+use transcript::IOPTranscript;
+
+#[cfg(feature = "parallel")]
+use rayon::iter::{IndexedParallelIterator, IntoParallelIterator, ParallelIterator};
+
+impl<F: PrimeField> SumCheckVerifier<F> for IOPVerifierState<F> {
+    type VPAuxInfo = VPAuxInfo<F>;
+    type ProverMessage = IOPProverMessage<F>;
+    type Challenge = F;
+    type Transcript = IOPTranscript<F>;
+    type SumCheckSubClaim = SumCheckSubClaim<F>;
+
+    /// Initialize the verifier's state.
+    fn verifier_init(index_info: &Self::VPAuxInfo) -> Self {
+        let start = start_timer!(|| "sum check verifier init");
+        let res = Self {
+            round: 1,
+            num_vars: index_info.num_variables,
+            max_degree: index_info.max_degree,
+            finished: false,
+            polynomials_received: Vec::with_capacity(index_info.num_variables),
+            challenges: Vec::with_capacity(index_info.num_variables),
+        };
+        end_timer!(start);
+        res
+    }
+
+    /// Run verifier for the current round, given a prover message.
+    ///
+    /// Note that `verify_round_and_update_state` only samples and stores
+    /// challenges; and update the verifier's state accordingly. The actual
+    /// verifications are deferred (in batch) to `check_and_generate_subclaim`
+    /// at the last step.
+    fn verify_round_and_update_state(
+        &mut self,
+        prover_msg: &Self::ProverMessage,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::Challenge, PolyIOPErrors> {
+        let start =
+            start_timer!(|| format!("sum check verify {}-th round and update state", self.round));
+
+        if self.finished {
+            return Err(PolyIOPErrors::InvalidVerifier(
+                "Incorrect verifier state: Verifier is already finished.".to_string(),
+            ));
+        }
+
+        // In an interactive protocol, the verifier should
+        //
+        // 1. check if the received 'P(0) + P(1) = expected`.
+        // 2. set `expected` to P(r)`
+        //
+        // When we turn the protocol to a non-interactive one, it is sufficient to defer
+        // such checks to `check_and_generate_subclaim` after the last round.
+
+        let challenge = transcript.get_and_append_challenge(b"Internal round")?;
+        self.challenges.push(challenge);
+        self.polynomials_received
+            .push(prover_msg.evaluations.to_vec());
+
+        if self.round == self.num_vars {
+            // accept and close
+            self.finished = true;
+        } else {
+            // proceed to the next round
+            self.round += 1;
+        }
+
+        end_timer!(start);
+        Ok(challenge)
+    }
+
+    /// This function verifies the deferred checks in the interactive version of
+    /// the protocol; and generate the subclaim. Returns an error if the
+    /// proof failed to verify.
+    ///
+    /// If the asserted sum is correct, then the multilinear polynomial
+    /// evaluated at `subclaim.point` will be `subclaim.expected_evaluation`.
+    /// Otherwise, it is highly unlikely that those two will be equal.
+    /// Larger field size guarantees smaller soundness error.
+    fn check_and_generate_subclaim(
+        &self,
+        asserted_sum: &F,
+    ) -> Result<Self::SumCheckSubClaim, PolyIOPErrors> {
+        let start = start_timer!(|| "sum check check and generate subclaim");
+        if !self.finished {
+            return Err(PolyIOPErrors::InvalidVerifier(
+                "Incorrect verifier state: Verifier has not finished.".to_string(),
+            ));
+        }
+
+        if self.polynomials_received.len() != self.num_vars {
+            return Err(PolyIOPErrors::InvalidVerifier(
+                "insufficient rounds".to_string(),
+            ));
+        }
+
+        // the deferred check during the interactive phase:
+        // 2. set `expected` to P(r)`
+        #[cfg(feature = "parallel")]
+        let mut expected_vec = self
+            .polynomials_received
+            .clone()
+            .into_par_iter()
+            .zip(self.challenges.clone().into_par_iter())
+            .map(|(evaluations, challenge)| {
+                if evaluations.len() != self.max_degree + 1 {
+                    return Err(PolyIOPErrors::InvalidVerifier(format!(
+                        "incorrect number of evaluations: {} vs {}",
+                        evaluations.len(),
+                        self.max_degree + 1
+                    )));
+                }
+                interpolate_uni_poly::<F>(&evaluations, challenge)
+            })
+            .collect::<Result<Vec<_>, PolyIOPErrors>>()?;
+
+        #[cfg(not(feature = "parallel"))]
+        let mut expected_vec = self
+            .polynomials_received
+            .clone()
+            .into_iter()
+            .zip(self.challenges.clone().into_iter())
+            .map(|(evaluations, challenge)| {
+                if evaluations.len() != self.max_degree + 1 {
+                    return Err(PolyIOPErrors::InvalidVerifier(format!(
+                        "incorrect number of evaluations: {} vs {}",
+                        evaluations.len(),
+                        self.max_degree + 1
+                    )));
+                }
+                interpolate_uni_poly::<F>(&evaluations, challenge)
+            })
+            .collect::<Result<Vec<_>, PolyIOPErrors>>()?;
+
+        // insert the asserted_sum to the first position of the expected vector
+        expected_vec.insert(0, *asserted_sum);
+
+        for (evaluations, &expected) in self
+            .polynomials_received
+            .iter()
+            .zip(expected_vec.iter())
+            .take(self.num_vars)
+        {
+            // the deferred check during the interactive phase:
+            // 1. check if the received 'P(0) + P(1) = expected`.
+            if evaluations[0] + evaluations[1] != expected {
+                return Err(PolyIOPErrors::InvalidProof(
+                    "Prover message is not consistent with the claim.".to_string(),
+                ));
+            }
+        }
+        end_timer!(start);
+        Ok(SumCheckSubClaim {
+            point: self.challenges.clone(),
+            // the last expected value (not checked within this function) will be included in the
+            // subclaim
+            expected_evaluation: expected_vec[self.num_vars],
+        })
+    }
+}
+
+/// Interpolate a uni-variate degree-`p_i.len()-1` polynomial and evaluate this
+/// polynomial at `eval_at`:
+///
+///   \sum_{i=0}^len p_i * (\prod_{j!=i} (eval_at - j)/(i-j) )
+///
+/// This implementation is linear in number of inputs in terms of field
+/// operations. It also has a quadratic term in primitive operations which is
+/// negligible compared to field operations.
+fn interpolate_uni_poly<F: PrimeField>(p_i: &[F], eval_at: F) -> Result<F, PolyIOPErrors> {
+    let start = start_timer!(|| "sum check interpolate uni poly opt");
+
+    let len = p_i.len();
+    let mut evals = vec![];
+    let mut prod = eval_at;
+    evals.push(eval_at);
+
+    // `prod = \prod_{j} (eval_at - j)`
+    for e in 1..len {
+        let tmp = eval_at - F::from(e as u64);
+        evals.push(tmp);
+        prod *= tmp;
+    }
+    let mut res = F::zero();
+    // we want to compute \prod (j!=i) (i-j) for a given i
+    //
+    // we start from the last step, which is
+    //  denom[len-1] = (len-1) * (len-2) *... * 2 * 1
+    // the step before that is
+    //  denom[len-2] = (len-2) * (len-3) * ... * 2 * 1 * -1
+    // and the step before that is
+    //  denom[len-3] = (len-3) * (len-4) * ... * 2 * 1 * -1 * -2
+    //
+    // i.e., for any i, the one before this will be derived from
+    //  denom[i-1] = denom[i] * (len-i) / i
+    //
+    // that is, we only need to store
+    // - the last denom for i = len-1, and
+    // - the ratio between current step and fhe last step, which is the product of
+    //   (len-i) / i from all previous steps and we store this product as a fraction
+    //   number to reduce field divisions.
+
+    // We know
+    //  - 2^61 < factorial(20) < 2^62
+    //  - 2^122 < factorial(33) < 2^123
+    // so we will be able to compute the ratio
+    //  - for len <= 20 with i64
+    //  - for len <= 33 with i128
+    //  - for len >  33 with BigInt
+    if p_i.len() <= 20 {
+        let last_denominator = F::from(u64_factorial(len - 1));
+        let mut ratio_numerator = 1i64;
+        let mut ratio_denominator = 1u64;
+
+        for i in (0..len).rev() {
+            let ratio_numerator_f = if ratio_numerator < 0 {
+                -F::from((-ratio_numerator) as u64)
+            } else {
+                F::from(ratio_numerator as u64)
+            };
+
+            res += p_i[i] * prod * F::from(ratio_denominator)
+                / (last_denominator * ratio_numerator_f * evals[i]);
+
+            // compute denom for the next step is current_denom * (len-i)/i
+            if i != 0 {
+                ratio_numerator *= -(len as i64 - i as i64);
+                ratio_denominator *= i as u64;
+            }
+        }
+    } else if p_i.len() <= 33 {
+        let last_denominator = F::from(u128_factorial(len - 1));
+        let mut ratio_numerator = 1i128;
+        let mut ratio_denominator = 1u128;
+
+        for i in (0..len).rev() {
+            let ratio_numerator_f = if ratio_numerator < 0 {
+                -F::from((-ratio_numerator) as u128)
+            } else {
+                F::from(ratio_numerator as u128)
+            };
+
+            res += p_i[i] * prod * F::from(ratio_denominator)
+                / (last_denominator * ratio_numerator_f * evals[i]);
+
+            // compute denom for the next step is current_denom * (len-i)/i
+            if i != 0 {
+                ratio_numerator *= -(len as i128 - i as i128);
+                ratio_denominator *= i as u128;
+            }
+        }
+    } else {
+        let mut denom_up = field_factorial::<F>(len - 1);
+        let mut denom_down = F::one();
+
+        for i in (0..len).rev() {
+            res += p_i[i] * prod * denom_down / (denom_up * evals[i]);
+
+            // compute denom for the next step is current_denom * (len-i)/i
+            if i != 0 {
+                denom_up *= -F::from((len - i) as u64);
+                denom_down *= F::from(i as u64);
+            }
+        }
+    }
+    end_timer!(start);
+    Ok(res)
+}
+
+/// compute the factorial(a) = 1 * 2 * ... * a
+#[inline]
+fn field_factorial<F: PrimeField>(a: usize) -> F {
+    let mut res = F::one();
+    for i in 2..=a {
+        res *= F::from(i as u64);
+    }
+    res
+}
+
+/// compute the factorial(a) = 1 * 2 * ... * a
+#[inline]
+fn u128_factorial(a: usize) -> u128 {
+    let mut res = 1u128;
+    for i in 2..=a {
+        res *= i as u128;
+    }
+    res
+}
+
+/// compute the factorial(a) = 1 * 2 * ... * a
+#[inline]
+fn u64_factorial(a: usize) -> u64 {
+    let mut res = 1u64;
+    for i in 2..=a {
+        res *= i as u64;
+    }
+    res
+}
+
+#[cfg(test)]
+mod test {
+    use super::interpolate_uni_poly;
+    use crate::poly_iop::errors::PolyIOPErrors;
+    use ark_bls12_381::Fr;
+    use ark_poly::{univariate::DensePolynomial, Polynomial, UVPolynomial};
+    use ark_std::{vec::Vec, UniformRand};
+
+    #[test]
+    fn test_interpolation() -> Result<(), PolyIOPErrors> {
+        let mut prng = ark_std::test_rng();
+
+        // test a polynomial with 20 known points, i.e., with degree 19
+        let poly = DensePolynomial::<Fr>::rand(20 - 1, &mut prng);
+        let evals = (0..20)
+            .map(|i| poly.evaluate(&Fr::from(i)))
+            .collect::<Vec<Fr>>();
+        let query = Fr::rand(&mut prng);
+
+        assert_eq!(poly.evaluate(&query), interpolate_uni_poly(&evals, query)?);
+
+        // test a polynomial with 33 known points, i.e., with degree 32
+        let poly = DensePolynomial::<Fr>::rand(33 - 1, &mut prng);
+        let evals = (0..33)
+            .map(|i| poly.evaluate(&Fr::from(i)))
+            .collect::<Vec<Fr>>();
+        let query = Fr::rand(&mut prng);
+
+        assert_eq!(poly.evaluate(&query), interpolate_uni_poly(&evals, query)?);
+
+        // test a polynomial with 64 known points, i.e., with degree 63
+        let poly = DensePolynomial::<Fr>::rand(64 - 1, &mut prng);
+        let evals = (0..64)
+            .map(|i| poly.evaluate(&Fr::from(i)))
+            .collect::<Vec<Fr>>();
+        let query = Fr::rand(&mut prng);
+
+        assert_eq!(poly.evaluate(&query), interpolate_uni_poly(&evals, query)?);
+
+        Ok(())
+    }
+}

subroutines/src/poly_iop/utils.rs

@@ -0,0 +1,28 @@
+//! useful macros.
+
+/// Takes as input a struct, and converts them to a series of bytes. All traits
+/// that implement `CanonicalSerialize` can be automatically converted to bytes
+/// in this manner.
+#[macro_export]
+macro_rules! to_bytes {
+    ($x:expr) => {{
+        let mut buf = ark_std::vec![];
+        ark_serialize::CanonicalSerialize::serialize($x, &mut buf).map(|_| buf)
+    }};
+}
+
+#[cfg(test)]
+mod test {
+    use ark_bls12_381::Fr;
+    use ark_serialize::CanonicalSerialize;
+    use ark_std::One;
+
+    #[test]
+    fn test_to_bytes() {
+        let f1 = Fr::one();
+
+        let mut bytes = ark_std::vec![];
+        f1.serialize(&mut bytes).unwrap();
+        assert_eq!(bytes, to_bytes!(&f1).unwrap());
+    }
+}

subroutines/src/poly_iop/zero_check/mod.rs

@@ -0,0 +1,205 @@
+//! Main module for the ZeroCheck protocol.
+
+use std::fmt::Debug;
+
+use crate::poly_iop::{errors::PolyIOPErrors, sum_check::SumCheck, PolyIOP};
+use arithmetic::build_eq_x_r;
+use ark_ff::PrimeField;
+use ark_poly::MultilinearExtension;
+use ark_std::{end_timer, start_timer};
+use transcript::IOPTranscript;
+
+/// A zero check IOP subclaim for `f(x)` consists of the following:
+///   - the initial challenge vector r which is used to build eq(x, r) in
+///     SumCheck
+///   - the random vector `v` to be evaluated
+///   - the claimed evaluation of `f(v)`
+#[derive(Clone, Debug, Default, PartialEq, Eq)]
+pub struct ZeroCheckSubClaim<F: PrimeField> {
+    // the evaluation point
+    pub point: Vec<F>,
+    /// the expected evaluation
+    pub expected_evaluation: F,
+    // the initial challenge r which is used to build eq(x, r)
+    pub init_challenge: Vec<F>,
+}
+
+/// A ZeroCheck for `f(x)` proves that `f(x) = 0` for all `x \in {0,1}^num_vars`
+/// It is derived from SumCheck.
+pub trait ZeroCheck<F: PrimeField>: SumCheck<F> {
+    type ZeroCheckSubClaim: Clone + Debug + Default + PartialEq;
+    type ZeroCheckProof: Clone + Debug + Default + PartialEq;
+
+    /// Initialize the system with a transcript
+    ///
+    /// This function is optional -- in the case where a ZeroCheck is
+    /// an building block for a more complex protocol, the transcript
+    /// may be initialized by this complex protocol, and passed to the
+    /// ZeroCheck prover/verifier.
+    fn init_transcript() -> Self::Transcript;
+
+    /// initialize the prover to argue for the sum of polynomial over
+    /// {0,1}^`num_vars` is zero.
+    fn prove(
+        poly: &Self::VirtualPolynomial,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::ZeroCheckProof, PolyIOPErrors>;
+
+    /// verify the claimed sum using the proof
+    fn verify(
+        proof: &Self::ZeroCheckProof,
+        aux_info: &Self::VPAuxInfo,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::ZeroCheckSubClaim, PolyIOPErrors>;
+}
+
+impl<F: PrimeField> ZeroCheck<F> for PolyIOP<F> {
+    type ZeroCheckSubClaim = ZeroCheckSubClaim<F>;
+    type ZeroCheckProof = Self::SumCheckProof;
+
+    fn init_transcript() -> Self::Transcript {
+        IOPTranscript::<F>::new(b"Initializing ZeroCheck transcript")
+    }
+
+    fn prove(
+        poly: &Self::VirtualPolynomial,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::ZeroCheckProof, PolyIOPErrors> {
+        let start = start_timer!(|| "zero check prove");
+
+        let length = poly.aux_info.num_variables;
+        let r = transcript.get_and_append_challenge_vectors(b"0check r", length)?;
+        let f_hat = poly.build_f_hat(r.as_ref())?;
+        let res = <Self as SumCheck<F>>::prove(&f_hat, transcript);
+
+        end_timer!(start);
+        res
+    }
+
+    fn verify(
+        proof: &Self::ZeroCheckProof,
+        fx_aux_info: &Self::VPAuxInfo,
+        transcript: &mut Self::Transcript,
+    ) -> Result<Self::ZeroCheckSubClaim, PolyIOPErrors> {
+        let start = start_timer!(|| "zero check verify");
+
+        // check that the sum is zero
+        if proof.proofs[0].evaluations[0] + proof.proofs[0].evaluations[1] != F::zero() {
+            return Err(PolyIOPErrors::InvalidProof(format!(
+                "zero check: sum {} is not zero",
+                proof.proofs[0].evaluations[0] + proof.proofs[0].evaluations[1]
+            )));
+        }
+
+        // generate `r` and pass it to the caller for correctness check
+        let length = fx_aux_info.num_variables;
+        let r = transcript.get_and_append_challenge_vectors(b"0check r", length)?;
+
+        // hat_fx's max degree is increased by eq(x, r).degree() which is 1
+        let mut hat_fx_aux_info = fx_aux_info.clone();
+        hat_fx_aux_info.max_degree += 1;
+        let sum_subclaim =
+            <Self as SumCheck<F>>::verify(F::zero(), proof, &hat_fx_aux_info, transcript)?;
+
+        // expected_eval = sumcheck.expect_eval/eq(v, r)
+        // where v = sum_check_sub_claim.point
+        let eq_x_r = build_eq_x_r(&r)?;
+        let expected_evaluation = sum_subclaim.expected_evaluation
+            / eq_x_r.evaluate(&sum_subclaim.point).ok_or_else(|| {
+                PolyIOPErrors::InvalidParameters("evaluation dimension does not match".to_string())
+            })?;
+
+        end_timer!(start);
+        Ok(ZeroCheckSubClaim {
+            point: sum_subclaim.point,
+            expected_evaluation,
+            init_challenge: r,
+        })
+    }
+}
+
+#[cfg(test)]
+mod test {
+
+    use super::ZeroCheck;
+    use crate::poly_iop::{errors::PolyIOPErrors, PolyIOP};
+    use arithmetic::VirtualPolynomial;
+    use ark_bls12_381::Fr;
+    use ark_std::test_rng;
+
+    fn test_zerocheck(
+        nv: usize,
+        num_multiplicands_range: (usize, usize),
+        num_products: usize,
+    ) -> Result<(), PolyIOPErrors> {
+        let mut rng = test_rng();
+
+        {
+            // good path: zero virtual poly
+            let poly =
+                VirtualPolynomial::rand_zero(nv, num_multiplicands_range, num_products, &mut rng)?;
+
+            let mut transcript = <PolyIOP<Fr> as ZeroCheck<Fr>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+            let proof = <PolyIOP<Fr> as ZeroCheck<Fr>>::prove(&poly, &mut transcript)?;
+
+            let poly_info = poly.aux_info.clone();
+            let mut transcript = <PolyIOP<Fr> as ZeroCheck<Fr>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+            let zero_subclaim =
+                <PolyIOP<Fr> as ZeroCheck<Fr>>::verify(&proof, &poly_info, &mut transcript)?;
+            assert!(
+                poly.evaluate(&zero_subclaim.point)? == zero_subclaim.expected_evaluation,
+                "wrong subclaim"
+            );
+        }
+
+        {
+            // bad path: random virtual poly whose sum is not zero
+            let (poly, _sum) =
+                VirtualPolynomial::<Fr>::rand(nv, num_multiplicands_range, num_products, &mut rng)?;
+
+            let mut transcript = <PolyIOP<Fr> as ZeroCheck<Fr>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+            let proof = <PolyIOP<Fr> as ZeroCheck<Fr>>::prove(&poly, &mut transcript)?;
+
+            let poly_info = poly.aux_info.clone();
+            let mut transcript = <PolyIOP<Fr> as ZeroCheck<Fr>>::init_transcript();
+            transcript.append_message(b"testing", b"initializing transcript for testing")?;
+
+            assert!(
+                <PolyIOP<Fr> as ZeroCheck<Fr>>::verify(&proof, &poly_info, &mut transcript)
+                    .is_err()
+            );
+        }
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_trivial_polynomial() -> Result<(), PolyIOPErrors> {
+        let nv = 1;
+        let num_multiplicands_range = (4, 5);
+        let num_products = 1;
+
+        test_zerocheck(nv, num_multiplicands_range, num_products)
+    }
+    #[test]
+    fn test_normal_polynomial() -> Result<(), PolyIOPErrors> {
+        let nv = 5;
+        let num_multiplicands_range = (4, 9);
+        let num_products = 5;
+
+        test_zerocheck(nv, num_multiplicands_range, num_products)
+    }
+
+    #[test]
+    fn zero_polynomial_should_error() -> Result<(), PolyIOPErrors> {
+        let nv = 0;
+        let num_multiplicands_range = (4, 13);
+        let num_products = 5;
+
+        assert!(test_zerocheck(nv, num_multiplicands_range, num_products).is_err());
+        Ok(())
+    }
+}