implemented parallelized kate division, with a lot of TODOs

Cargo.toml

@@ -32,8 +32,8 @@ git = "https://github.com/gtank/blake2-rfc"
 rev = "7a5b5fc99ae483a0043db7547fb79a6fa44b88a9"
 
 [features]
-default = ["multicore"]
+#default = ["multicore"]
-#default = ["multicore", "gm17", "sonic"]
+default = ["multicore", "gm17", "sonic"]
 #default = ["singlecore"]
 multicore = ["futures-cpupool", "num_cpus", "crossbeam"]
 sonic = ["tiny-keccak"]

src/sonic/README.md

@@ -20,8 +20,14 @@ Initial SONIC proof system integration using the code from the [original impleme
 - [x] Basic Ethereum smart-contract
 - [x] Add blinding factors
 - [ ] Implement unhelped version
-  - [x] Implement a part of S poly precomputation (S2)
+  - [x] Implement a part of S poly precomputation (S_2)
   - [x] Implement a "well formed" argument 
   - [x] Implement a coefficients product argument
-  - [ ] Implement a premutation argument
+  - [x] Implement a premutation argument
-  - [ ] Implement synthesizer for proper form of S polynomial
+  - [ ] Implement synthesizer for proper form of S polynomial
+- [ ] Finer tuning
+  - [x] Parallelized evaluation of S polynomial
+  - [ ] Parallelize some polynomial operations in the protocol itself
+  - [x] Parallelized kate division  
+  - [ ] Implement specialized version of polynomial multiplication by degree 1 polynomial
+  - [ ] Non-assigning variant of the "adaptor" (may require a change of SONIC CS trait)

src/sonic/cs/mod.rs

@@ -209,50 +209,6 @@ impl SynthesisDriver for Basic {
 
         circuit.synthesize(&mut tmp)?;
 
-        // let blindings_to_add = 6;
-
-        // for i in 0..blindings_to_add {
-
-        //     match tmp.current_variable.take() {
-        //         Some(index) => {
-        //             let var_a = Variable::A(index);
-        //             let var_b = Variable::B(index);
-        //             let var_c = Variable::C(index);
-
-        //             let mut product = None;
-
-        //             let value_a = tmp.backend.get_var(var_a);
-
-        //             tmp.backend.set_var(var_b, || {
-        //                 let value_b = E::Fr::one();
-        //                 product = Some(value_a.ok_or(SynthesisError::AssignmentMissing)?);
-        //                 product.as_mut().map(|product| product.mul_assign(&value_b));
-
-        //                 Ok(value_b)
-        //             })?;
-
-        //             tmp.backend.set_var(var_c, || {
-        //                 product.ok_or(SynthesisError::AssignmentMissing)
-        //             })?;
-
-        //             tmp.current_variable = None;
-        //         },
-        //         None => {
-        //             self.n += 1;
-        //             let index = self.n ;
-        //             tmp.backend.new_multiplication_gate();
-
-        //             let var_a = Variable::A(index);
-
-        //             tmp.backend.set_var(var_a, value)?;
-
-        //             tmp.current_variable = Some(index);
-        //         }
-        //     }
-        // }
-
-        // TODO: add blinding factors so we actually get zero-knowledge
-
         Ok(())
     }
 }
@@ -357,10 +313,6 @@ impl SynthesisDriver for Nonassigning {
 
             circuit.synthesize(&mut tmp)?;
 
-            // TODO: add blinding factors so we actually get zero-knowledge
-
-            // println!("n = {}", tmp.n);
-
             Ok(())
     }
 }

src/sonic/helped/generator.rs

@@ -383,7 +383,6 @@ pub fn generate_parameters<E, C>(
 {
     let circuit_parameters = get_circuit_parameters::<E, C>(circuit)?; 
     let min_d = circuit_parameters.n * 4 + 2*NUM_BLINDINGS;
-    println!{"Mid d = {}", min_d};
 
     let srs = generate_srs(alpha, x, min_d)?;
 

src/sonic/helped/helper.rs

@@ -62,7 +62,7 @@ pub fn create_aggregate_on_srs<E: Engine, C: Circuit<E>, S: SynthesisDriver>(
             }
         }
 
-        let mut tmp = CountN{n:0,q:0};
+        let mut tmp = CountN{n:0, q:0};
         S::synthesize(&mut tmp, circuit).unwrap(); // TODO
 
         (tmp.n, tmp.q)
@@ -139,23 +139,28 @@ pub fn create_aggregate_on_srs_using_information<E: Engine, C: Circuit<E>, S: Sy
     // Let's open up C to every y.
     fn compute_value<E: Engine>(y: &E::Fr, poly_positive: &[E::Fr], poly_negative: &[E::Fr]) -> E::Fr {
         let mut value = E::Fr::zero();
-
         let yinv = y.inverse().unwrap(); // TODO
-        let mut tmp = yinv;
-        for &coeff in poly_negative {
-            let mut coeff = coeff;
-            coeff.mul_assign(&tmp);
-            value.add_assign(&coeff);
-            tmp.mul_assign(&yinv);
-        }
 
-        let mut tmp = *y;
+        let positive_powers_contrib = evaluate_at_consequitive_powers(poly_positive, *y, *y);
-        for &coeff in poly_positive {
+        let negative_powers_contrib = evaluate_at_consequitive_powers(poly_negative, yinv, yinv);
-            let mut coeff = coeff;
+        value.add_assign(&positive_powers_contrib);
-            coeff.mul_assign(&tmp);
+        value.add_assign(&negative_powers_contrib);
-            value.add_assign(&coeff);
+
-            tmp.mul_assign(&y);
+        // let mut tmp = yinv;
-        }
+        // for &coeff in poly_negative {
+        //     let mut coeff = coeff;
+        //     coeff.mul_assign(&tmp);
+        //     value.add_assign(&coeff);
+        //     tmp.mul_assign(&yinv);
+        // }
+
+        // let mut tmp = *y;
+        // for &coeff in poly_positive {
+        //     let mut coeff = coeff;
+        //     coeff.mul_assign(&tmp);
+        //     value.add_assign(&coeff);
+        //     tmp.mul_assign(&y);
+        // }
 
         value
     }

src/sonic/helped/poly.rs

@@ -4,6 +4,7 @@ use std::marker::PhantomData;
 
 use crate::sonic::cs::{Backend};
 use crate::sonic::cs::{Coeff, Variable, LinearCombination};
+use crate::sonic::util::*;
 
 /*
 s(X, Y) =   \sum\limits_{i=1}^N u_i(Y) X^{-i}
@@ -71,26 +72,33 @@ impl<E: Engine> SxEval<E> {
     pub fn finalize(self, x: E::Fr) -> E::Fr {
         let x_inv = x.inverse().unwrap(); // TODO
 
-        let mut tmp = x_inv;
-
         let mut acc = E::Fr::zero();
-        for mut u in self.u {
-            u.mul_assign(&tmp);
-            acc.add_assign(&u);
-            tmp.mul_assign(&x_inv);
-        }
 
+        let mut tmp = x_inv;
+        acc.add_assign(&evaluate_at_consequitive_powers(& self.u[..], tmp, tmp));
         let mut tmp = x;
-        for mut v in self.v {
+        acc.add_assign(&evaluate_at_consequitive_powers(& self.v[..], tmp, tmp));
-            v.mul_assign(&tmp);
+        let mut tmp = x.pow(&[(self.v.len()+1) as u64]);
-            acc.add_assign(&v);
+        acc.add_assign(&evaluate_at_consequitive_powers(& self.w[..], tmp, x));
-            tmp.mul_assign(&x);
+
-        }
+        // let mut tmp = x_inv;
-        for mut w in self.w {
+        // for mut u in self.u {
-            w.mul_assign(&tmp);
+        //     u.mul_assign(&tmp);
-            acc.add_assign(&w);
+        //     acc.add_assign(&u);
-            tmp.mul_assign(&x);
+        //     tmp.mul_assign(&x_inv);
-        }
+        // }
+
+        // let mut tmp = x;
+        // for mut v in self.v {
+        //     v.mul_assign(&tmp);
+        //     acc.add_assign(&v);
+        //     tmp.mul_assign(&x);
+        // }
+        // for mut w in self.w {
+        //     w.mul_assign(&tmp);
+        //     acc.add_assign(&w);
+        //     tmp.mul_assign(&x);
+        // }
 
         acc
     }
@@ -209,20 +217,26 @@ impl<E: Engine> SyEval<E> {
 
     pub fn finalize(self, y: E::Fr) -> E::Fr {
         let mut acc = E::Fr::zero();
-
-        let mut tmp = y;
-        for mut coeff in self.positive_coeffs {
-            coeff.mul_assign(&tmp);
-            acc.add_assign(&coeff);
-            tmp.mul_assign(&y);
-        }
         let yinv = y.inverse().unwrap(); // TODO
-        let mut tmp = yinv;
+
-        for mut coeff in self.negative_coeffs {
+        let positive_powers_contrib = evaluate_at_consequitive_powers(& self.positive_coeffs[..], y, y);
-            coeff.mul_assign(&tmp);
+        let negative_powers_contrib = evaluate_at_consequitive_powers(& self.negative_coeffs[..], yinv, yinv);
-            acc.add_assign(&coeff);
+        acc.add_assign(&positive_powers_contrib);
-            tmp.mul_assign(&yinv);
+        acc.add_assign(&negative_powers_contrib);
-        }
+
+        // let mut tmp = y;
+        // for mut coeff in self.positive_coeffs {
+        //     coeff.mul_assign(&tmp);
+        //     acc.add_assign(&coeff);
+        //     tmp.mul_assign(&y);
+        // }
+
+        // let mut tmp = yinv;
+        // for mut coeff in self.negative_coeffs {
+        //     coeff.mul_assign(&tmp);
+        //     acc.add_assign(&coeff);
+        //     tmp.mul_assign(&yinv);
+        // }
 
         acc
     }

src/sonic/helped/prover.rs

@@ -310,23 +310,6 @@ pub fn create_proof_on_srs<E: Engine, C: Circuit<E>, S: SynthesisDriver>(
         evaluate_at_consequitive_powers(&rx1, tmp, z)
     };
 
-    // {
-    //     rx1[(2 * n + NUM_BLINDINGS)].sub_assign(&rz);
-
-    //     let opening = polynomial_commitment_opening(
-    //         2 * n + NUM_BLINDINGS,
-    //         n, 
-    //         rx1.iter(),
-    //         z,
-    //         srs
-    //     );
-
-    //     let valid_rz_commitment = check_polynomial_commitment(&r, &z, &rz, &opening, n, &srs);
-    //     assert!(valid_rz_commitment);
-
-    //     rx1[(2 * n + NUM_BLINDINGS)].add_assign(&rz);
-    // }
-
     // rzy is evaluation of r(X, Y) at z, y
     let rzy = {
         let tmp = z_inv.pow(&[(2*n + NUM_BLINDINGS) as u64]);
@@ -385,11 +368,6 @@ pub fn create_proof_on_srs<E: Engine, C: Circuit<E>, S: SynthesisDriver>(
             srs)
     };
 
-    // let mut zy = z;
-    // zy.mul_assign(&y);
-    // let valid_rzy_commitment = check_polynomial_commitment(&r, &zy, &rzy, &zy_opening, n, &srs);
-    // assert!(valid_rzy_commitment);
-
     Ok(Proof {
         r, rz, rzy, t, z_opening, zy_opening
     })

src/sonic/tests/mod.rs

@@ -0,0 +1 @@
+mod sonics;

src/sonic/tests/sonics.rs

@@ -1,5 +1,3 @@
-extern crate bellman;
-extern crate pairing;
 extern crate rand;
 
 // For randomness (during paramgen and proof generation)
@@ -27,23 +25,131 @@ use pairing::bn256::{
 };
 
 // We'll use these interfaces to construct our circuit.
-use bellman::{
+use crate::{
     Circuit,
     ConstraintSystem,
     SynthesisError
 };
 
-// We're going to use the Groth16 proving system.
-use bellman::groth16::{
-    Proof,
-    generate_random_parameters,
-    prepare_verifying_key,
-    create_random_proof,
-    verify_proof,
-};
-
 const MIMC_ROUNDS: usize = 322;
 
+fn mimc<E: Engine>(
+    mut xl: E::Fr,
+    mut xr: E::Fr,
+    constants: &[E::Fr]
+) -> E::Fr
+{
+    assert_eq!(constants.len(), MIMC_ROUNDS);
+
+    for i in 0..MIMC_ROUNDS {
+        let mut tmp1 = xl;
+        tmp1.add_assign(&constants[i]);
+        let mut tmp2 = tmp1;
+        tmp2.square();
+        tmp2.mul_assign(&tmp1);
+        tmp2.add_assign(&xr);
+        xr = xl;
+        xl = tmp2;
+    }
+
+    xl
+}
+
+/// This is our demo circuit for proving knowledge of the
+/// preimage of a MiMC hash invocation.
+#[derive(Clone)]
+struct MiMCDemo<'a, E: Engine> {
+    xl: Option<E::Fr>,
+    xr: Option<E::Fr>,
+    constants: &'a [E::Fr]
+}
+
+/// Our demo circuit implements this `Circuit` trait which
+/// is used during paramgen and proving in order to
+/// synthesize the constraint system.
+impl<'a, E: Engine> Circuit<E> for MiMCDemo<'a, E> {
+    fn synthesize<CS: ConstraintSystem<E>>(
+        self,
+        cs: &mut CS
+    ) -> Result<(), SynthesisError>
+    {
+        assert_eq!(self.constants.len(), MIMC_ROUNDS);
+
+        // Allocate the first component of the preimage.
+        let mut xl_value = self.xl;
+        let mut xl = cs.alloc(|| "preimage xl", || {
+            xl_value.ok_or(SynthesisError::AssignmentMissing)
+        })?;
+
+        // Allocate the second component of the preimage.
+        let mut xr_value = self.xr;
+        let mut xr = cs.alloc(|| "preimage xr", || {
+            xr_value.ok_or(SynthesisError::AssignmentMissing)
+        })?;
+
+        for i in 0..MIMC_ROUNDS {
+            // xL, xR := xR + (xL + Ci)^3, xL
+            let cs = &mut cs.namespace(|| format!("round {}", i));
+
+            // tmp = (xL + Ci)^2
+            let tmp_value = xl_value.map(|mut e| {
+                e.add_assign(&self.constants[i]);
+                e.square();
+                e
+            });
+            let tmp = cs.alloc(|| "tmp", || {
+                tmp_value.ok_or(SynthesisError::AssignmentMissing)
+            })?;
+
+            cs.enforce(
+                || "tmp = (xL + Ci)^2",
+                |lc| lc + xl + (self.constants[i], CS::one()),
+                |lc| lc + xl + (self.constants[i], CS::one()),
+                |lc| lc + tmp
+            );
+
+            // new_xL = xR + (xL + Ci)^3
+            // new_xL = xR + tmp * (xL + Ci)
+            // new_xL - xR = tmp * (xL + Ci)
+            let new_xl_value = xl_value.map(|mut e| {
+                e.add_assign(&self.constants[i]);
+                e.mul_assign(&tmp_value.unwrap());
+                e.add_assign(&xr_value.unwrap());
+                e
+            });
+
+            let new_xl = if i == (MIMC_ROUNDS-1) {
+                // This is the last round, xL is our image and so
+                // we allocate a public input.
+                cs.alloc_input(|| "image", || {
+                    new_xl_value.ok_or(SynthesisError::AssignmentMissing)
+                })?
+            } else {
+                cs.alloc(|| "new_xl", || {
+                    new_xl_value.ok_or(SynthesisError::AssignmentMissing)
+                })?
+            };
+
+            cs.enforce(
+                || "new_xL = xR + (xL + Ci)^3",
+                |lc| lc + tmp,
+                |lc| lc + xl + (self.constants[i], CS::one()),
+                |lc| lc + new_xl - xr
+            );
+
+            // xR = xL
+            xr = xl;
+            xr_value = xl_value;
+
+            // xL = new_xL
+            xl = new_xl;
+            xl_value = new_xl_value;
+        }
+
+        Ok(())
+    }
+}
+
 /// This is our demo circuit for proving knowledge of the
 /// preimage of a MiMC hash invocation.
 #[derive(Clone)]
@@ -147,7 +253,7 @@ fn test_sonic_mimc() {
     use pairing::{Engine, CurveAffine, CurveProjective};
     use pairing::bls12_381::{Bls12, Fr};
     use std::time::{Instant};
-    use bellman::sonic::srs::SRS;
+    use crate::sonic::srs::SRS;
 
     let srs_x = Fr::from_str("23923").unwrap();
     let srs_alpha = Fr::from_str("23728792").unwrap();
@@ -178,11 +284,11 @@ fn test_sonic_mimc() {
             constants: &constants
         };
 
-        use bellman::sonic::cs::Basic;
+        use crate::sonic::cs::Basic;
-        use bellman::sonic::sonic::AdaptorCircuit;
+        use crate::sonic::sonic::AdaptorCircuit;
-        use bellman::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs};
+        use crate::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs};
-        use bellman::sonic::helped::{MultiVerifier, get_circuit_parameters};
+        use crate::sonic::helped::{MultiVerifier, get_circuit_parameters};
-        use bellman::sonic::helped::helper::{create_aggregate_on_srs};
+        use crate::sonic::helped::helper::{create_aggregate_on_srs};
 
         println!("creating proof");
         let start = Instant::now();
@@ -252,7 +358,7 @@ fn test_inputs_into_sonic_mimc() {
     use pairing::bn256::{Bn256, Fr};
     // use pairing::bls12_381::{Bls12, Fr};
     use std::time::{Instant};
-    use bellman::sonic::srs::SRS;
+    use crate::sonic::srs::SRS;
 
     let srs_x = Fr::from_str("23923").unwrap();
     let srs_alpha = Fr::from_str("23728792").unwrap();
@@ -282,11 +388,11 @@ fn test_inputs_into_sonic_mimc() {
             constants: &constants
         };
 
-        use bellman::sonic::cs::Basic;
+        use crate::sonic::cs::Basic;
-        use bellman::sonic::sonic::AdaptorCircuit;
+        use crate::sonic::sonic::AdaptorCircuit;
-        use bellman::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs};
+        use crate::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs};
-        use bellman::sonic::helped::{MultiVerifier, get_circuit_parameters};
+        use crate::sonic::helped::{MultiVerifier, get_circuit_parameters};
-        use bellman::sonic::helped::helper::{create_aggregate_on_srs};
+        use crate::sonic::helped::helper::{create_aggregate_on_srs};
 
         let info = get_circuit_parameters::<Bn256, _>(circuit.clone()).expect("Must get circuit info");
         println!("{:?}", info);
@@ -356,7 +462,7 @@ fn test_inputs_into_sonic_mimc() {
 fn test_high_level_sonic_api() {
     use pairing::bn256::{Bn256};
     use std::time::{Instant};
-    use bellman::sonic::helped::{
+    use crate::sonic::helped::{
         generate_random_parameters, 
         verify_aggregate, 
         verify_proofs, 

src/sonic/util.rs

@@ -123,10 +123,12 @@ pub fn polynomial_commitment_opening<
     ) -> E::G1Affine
         where I::IntoIter: DoubleEndedIterator + ExactSizeIterator,
     {
-        let poly = kate_divison(
+        let poly = parallel_kate_divison::<E, _>(polynomial_coefficients, point);
-            polynomial_coefficients,
+
-            point,
+        // let poly = kate_divison(
-        );
+        //     polynomial_coefficients,
+        //     point,
+        // );
 
         let negative_poly = poly[0..largest_negative_power].iter().rev();
         let positive_poly = poly[largest_negative_power..].iter();
@@ -400,6 +402,75 @@ where
     q
 }
 
+/// Divides polynomial `a` in `x` by `x - b` with
+/// no remainder using fft.
+pub fn parallel_kate_divison<'a, E: Engine, I: IntoIterator<Item = &'a E::Fr>>(a: I, mut b: E::Fr) -> Vec<E::Fr>
+where
+    I::IntoIter: DoubleEndedIterator + ExactSizeIterator,
+{
+    // this implementation is only for division by `x - b` form polynomail,
+    // so we can manuall calculate the reciproical poly of the form `x^2/(x-b)`
+    // and the reminder
+
+    // x^2 /(x - b) = x + b*x/(x - b) = (x + b) + b^2/(x - b)
+
+    let reciproical = vec![b, E::Fr::one()]; // x + b
+
+    // and remainder b^2
+    let mut b_squared = b;
+    b_squared.square();
+
+    let mut b_neg = b;
+    b_neg.negate();
+
+    let divisor = vec![b_neg, E::Fr::one()];
+
+    let poly: Vec<E::Fr> = a.into_iter().map(|el| el.clone()).collect();
+
+    let (q, _) = kate_divison_inner::<E>(poly, divisor, reciproical, b_squared);
+
+    // assert_eq!(r.len(), 0);
+
+    q
+}
+
+fn kate_divison_inner<E: Engine>(
+        poly: Vec<E::Fr>, 
+        divisor: Vec<E::Fr>, 
+        reciproical: Vec<E::Fr>, 
+        remainder: E::Fr
+    ) -> (Vec<E::Fr>, Vec<E::Fr>) {
+    if poly.len() == 1 {
+        return (vec![], poly);
+    }
+    // TODO: Change generic multiplications by multiplications by degree 1 polynomial
+    let poly_degree = poly.len() - 1;
+    let mut q = multiply_polynomials::<E>(poly.clone(), reciproical.clone());
+    q.drain(0..2);
+    // recursion step
+    if poly_degree > 2 {
+        let mut rec_step = poly.clone();
+        mul_polynomial_by_scalar(&mut rec_step[..], remainder);
+        // truncate low order terms 
+        rec_step.drain(0..2);
+        let (q2, _) = kate_divison_inner::<E>(rec_step, divisor.clone(), reciproical, remainder);
+        // length of q2 is smaller
+        add_polynomials(&mut q[..q2.len()], &q2[..]);
+    }        
+
+    // although r must be zero, calculate it for now
+    if q.len() == 0 {
+        return (q, poly);
+    }
+
+    // r = u - v*q
+    let mut poly = poly;
+    let tmp = multiply_polynomials::<E>(divisor, q.clone());
+    sub_polynomials(&mut poly[..], &tmp[..]);
+
+    return (q, poly);
+}
+
 /// Convenience function to check polynomail commitment
 pub fn check_polynomial_commitment<E: Engine>(
     commitment: &E::G1Affine,
@@ -524,6 +595,90 @@ pub fn multiply_polynomials<E: Engine>(a: Vec<E::Fr>, b: Vec<E::Fr>) -> Vec<E::F
     mul_result
 }
 
+
+// alternative implementation that does not require an `Evaluation domain` struct
+pub fn multiply_polynomials_fft<E: Engine>(a: Vec<E::Fr>, b: Vec<E::Fr>) -> Vec<E::Fr> {
+    use crate::multicore::Worker;
+    use crate::domain::{best_fft, Scalar};
+    use crate::group::Group;
+
+    let result_len = a.len() + b.len() - 1;
+
+    // m is a size of domain where Z polynomial does NOT vanish
+    // in normal domain Z is in a form of (X-1)(X-2)...(X-N)
+    let mut m = 1;
+    let mut exp = 0;
+    let mut omega = E::Fr::root_of_unity();
+    let max_degree = (1 << E::Fr::S) - 1;
+
+    if result_len > max_degree {
+        panic!("multiplication result degree is too large");
+    }
+
+    while m < result_len {
+        m *= 2;
+        exp += 1;
+
+        // The pairing-friendly curve may not be able to support
+        // large enough (radix2) evaluation domains.
+        if exp > E::Fr::S {
+            panic!("multiplication result degree is too large");
+        }
+    }
+
+    // If full domain is not needed - limit it,
+    // e.g. if (2^N)th power is not required, just double omega and get 2^(N-1)th
+    // Compute omega, the 2^exp primitive root of unity
+    for _ in exp..E::Fr::S {
+        omega.square();
+    }
+
+    let omegainv = omega.inverse().unwrap();
+    let minv = E::Fr::from_str(&format!("{}", m)).unwrap().inverse().unwrap();
+
+    let worker = Worker::new();
+
+    let mut scalars_a: Vec<Scalar<E>> = a.into_iter().map(|e| Scalar::<E>(e)).collect();
+    let mut scalars_b: Vec<Scalar<E>> = b.into_iter().map(|e| Scalar::<E>(e)).collect();
+    scalars_a.resize(m, Scalar::<E>(E::Fr::zero()));
+    scalars_b.resize(m, Scalar::<E>(E::Fr::zero()));
+
+
+    best_fft(&mut scalars_a[..], &worker, &omega, exp);
+    best_fft(&mut scalars_b[..], &worker, &omega, exp);
+
+    // do the convolution
+    worker.scope(scalars_a.len(), |scope, chunk| {
+        for (a, b) in scalars_a.chunks_mut(chunk).zip(scalars_b.chunks(chunk)) {
+            scope.spawn(move |_| {
+                for (a, b) in a.iter_mut().zip(b.iter()) {
+                    a.group_mul_assign(&b.0);
+                }
+            });
+        }
+    });
+
+    // no longer need it
+    drop(scalars_b);
+
+    best_fft(&mut scalars_a[..], &worker, &omegainv, exp);
+    worker.scope(scalars_a.len(), |scope, chunk| {
+        for v in scalars_a.chunks_mut(chunk) {
+            scope.spawn(move |_| {
+                for v in v {
+                    v.group_mul_assign(&minv);
+                }
+            });
+        }
+    });
+
+    let mut mul_result: Vec<E::Fr> = scalars_a.into_iter().map(|e| e.0).collect();
+
+    mul_result.truncate(result_len);
+
+    mul_result
+}
+
 pub fn multiply_polynomials_serial<E: Engine>(mut a: Vec<E::Fr>, mut b: Vec<E::Fr>) -> Vec<E::Fr> {
     let result_len = a.len() + b.len() - 1;
 
@@ -574,6 +729,7 @@ pub fn multiply_polynomials_serial<E: Engine>(mut a: Vec<E::Fr>, mut b: Vec<E::F
     a
 }
 
+// add polynomails in coefficient form
 pub fn add_polynomials<F: Field>(a: &mut [F], b: &[F]) {
         use crate::multicore::Worker;
         use crate::domain::{EvaluationDomain, Scalar};
@@ -594,6 +750,28 @@ pub fn add_polynomials<F: Field>(a: &mut [F], b: &[F]) {
         });
 }
 
+// subtract polynomails in coefficient form
+pub fn sub_polynomials<F: Field>(a: &mut [F], b: &[F]) {
+    use crate::multicore::Worker;
+    use crate::domain::{EvaluationDomain, Scalar};
+
+    let worker = Worker::new();
+
+    assert_eq!(a.len(), b.len());
+
+    worker.scope(a.len(), |scope, chunk| {
+        for (a, b) in a.chunks_mut(chunk).zip(b.chunks(chunk))
+        {
+            scope.spawn(move |_| {
+                for (a, b) in a.iter_mut().zip(b.iter()) {
+                    a.sub_assign(b);
+                }
+            });
+        }
+    });
+}
+
+// multiply coefficients of the polynomial by the scalar
 pub fn mul_polynomial_by_scalar<F: Field>(a: &mut [F], b: F) {
         use crate::multicore::Worker;
         use crate::domain::{EvaluationDomain, Scalar};
@@ -612,6 +790,8 @@ pub fn mul_polynomial_by_scalar<F: Field>(a: &mut [F], b: F) {
         });
 }
 
+// elementwise add coeffs of one polynomial with coeffs of other, that are 
+// first multiplied by a scalar 
 pub fn mul_add_polynomials<F: Field>(a: &mut [F], b: &[F], c: F) {
         use crate::multicore::Worker;
         use crate::domain::{EvaluationDomain, Scalar};
@@ -693,7 +873,6 @@ impl<T> OptionExt<T> for Option<T> {
 }
 
 #[test]
-
 fn test_mul() {
     use rand::{self, Rand};
     use pairing::bls12_381::Bls12;
@@ -804,4 +983,119 @@ fn test_mut_distribute_powers() {
     mut_distribute_consequitive_powers(&mut b[..], first_power, x);
 
     assert!(a == b);
+}
+
+
+#[test]
+fn test_trivial_parallel_kate_division() {
+    use pairing::ff::PrimeField;
+    use pairing::bls12_381::{Bls12, Fr};
+
+    let mut minus_one = Fr::one();
+    minus_one.negate();
+
+    let z = Fr::one();
+
+    // this is x^2 - 1
+    let poly = vec![
+        minus_one,
+        Fr::from_str("0").unwrap(),
+        Fr::from_str("1").unwrap(),
+    ];
+
+    let quotient_poly = kate_divison(&poly, z);
+
+    let parallel_q_poly = parallel_kate_divison::<Bls12, _>(&poly, z);
+
+    assert_eq!(quotient_poly, parallel_q_poly);
+}
+
+#[test]
+fn test_less_trivial_parallel_kate_division() {
+    use pairing::ff::PrimeField;
+    use pairing::bls12_381::{Bls12, Fr};
+
+    let z = Fr::one();
+
+    let mut poly = vec![
+        Fr::from_str("328947234").unwrap(),
+        Fr::from_str("3545623451111").unwrap(),
+        Fr::from_str("5").unwrap(),
+        Fr::from_str("55555").unwrap(),
+        Fr::from_str("1235685").unwrap(),
+    ];
+
+    fn eval(poly: &[Fr], point: Fr) -> Fr {
+        let mut acc = Fr::zero();
+        let mut tmp = Fr::one();
+        for p in &poly[..] {
+            let mut t = *p;
+            t.mul_assign(&tmp);
+            acc.add_assign(&t);
+            tmp.mul_assign(&point);
+        }
+    
+        acc
+    }
+
+    let p_at_z = eval(&poly, z);
+
+    // poly = poly(X) - poly(z)
+    poly[0].sub_assign(&p_at_z);
+
+    let quotient_poly = kate_divison(&poly, z);
+
+    let parallel_q_poly = parallel_kate_divison::<Bls12, _>(&poly, z);
+
+    assert_eq!(quotient_poly, parallel_q_poly);
+}
+
+
+#[test]
+fn test_parallel_kate_division() {
+    use pairing::ff::PrimeField;
+    use pairing::bls12_381::{Bls12, Fr};
+
+    let mut poly = vec![
+        Fr::from_str("328947234").unwrap(),
+        Fr::from_str("3545623451111").unwrap(),
+        Fr::from_str("0").unwrap(),
+        Fr::from_str("55555").unwrap(),
+        Fr::from_str("1235685").unwrap(),
+    ];
+
+    fn eval(poly: &[Fr], point: Fr) -> Fr {
+        let point_inv = point.inverse().unwrap();
+
+        let mut acc = Fr::zero();
+        let mut tmp = Fr::one();
+        for p in &poly[2..] {
+            let mut t = *p;
+            t.mul_assign(&tmp);
+            acc.add_assign(&t);
+            tmp.mul_assign(&point);
+        }
+        let mut tmp = point_inv;
+        for p in poly[0..2].iter().rev() {
+            let mut t = *p;
+            t.mul_assign(&tmp);
+            acc.add_assign(&t);
+            tmp.mul_assign(&point_inv);
+        }
+
+        acc
+    }
+
+    let z = Fr::from_str("2000").unwrap();
+
+    let p_at_z = eval(&poly, z);
+
+    // poly = poly(X) - poly(z)
+    poly[2].sub_assign(&p_at_z);
+
+    let quotient_poly = kate_divison(&poly, z);
+
+    let parallel_q_poly = parallel_kate_divison::<Bls12, _>(&poly, z);
+
+    assert_eq!(quotient_poly, parallel_q_poly);
 }