From 2173354b1fa2465a48b967e603197a2b714a86df Mon Sep 17 00:00:00 2001 From: Alex Vlasov Date: Mon, 11 Mar 2019 14:57:14 +0100 Subject: [PATCH] implemented parallelized kate division, with a lot of TODOs --- Cargo.toml | 4 +- src/sonic/README.md | 12 +- src/sonic/cs/mod.rs | 48 ------ src/sonic/helped/generator.rs | 1 - src/sonic/helped/helper.rs | 37 +++-- src/sonic/helped/poly.rs | 74 +++++---- src/sonic/helped/prover.rs | 22 --- src/sonic/tests/mod.rs | 1 + src/sonic/tests/sonics.rs | 156 ++++++++++++++--- src/sonic/util.rs | 304 +++++++++++++++++++++++++++++++++- 10 files changed, 507 insertions(+), 152 deletions(-) create mode 100644 src/sonic/tests/mod.rs diff --git a/Cargo.toml b/Cargo.toml index b97cbf9..099343a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -32,8 +32,8 @@ git = "https://github.com/gtank/blake2-rfc" rev = "7a5b5fc99ae483a0043db7547fb79a6fa44b88a9" [features] -default = ["multicore"] -#default = ["multicore", "gm17", "sonic"] +#default = ["multicore"] +default = ["multicore", "gm17", "sonic"] #default = ["singlecore"] multicore = ["futures-cpupool", "num_cpus", "crossbeam"] sonic = ["tiny-keccak"] diff --git a/src/sonic/README.md b/src/sonic/README.md index df0f9c1..9d47757 100644 --- a/src/sonic/README.md +++ b/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 - - [ ] Implement synthesizer for proper form of S polynomial \ No newline at end of file + - [x] Implement a premutation argument + - [ ] 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) \ No newline at end of file diff --git a/src/sonic/cs/mod.rs b/src/sonic/cs/mod.rs index 0a42dfa..464739c 100644 --- a/src/sonic/cs/mod.rs +++ b/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(()) } } \ No newline at end of file diff --git a/src/sonic/helped/generator.rs b/src/sonic/helped/generator.rs index 652c501..1253373 100644 --- a/src/sonic/helped/generator.rs +++ b/src/sonic/helped/generator.rs @@ -383,7 +383,6 @@ pub fn generate_parameters( { let circuit_parameters = get_circuit_parameters::(circuit)?; let min_d = circuit_parameters.n * 4 + 2*NUM_BLINDINGS; - println!{"Mid d = {}", min_d}; let srs = generate_srs(alpha, x, min_d)?; diff --git a/src/sonic/helped/helper.rs b/src/sonic/helped/helper.rs index a8a853a..fb7bf96 100644 --- a/src/sonic/helped/helper.rs +++ b/src/sonic/helped/helper.rs @@ -62,7 +62,7 @@ pub fn create_aggregate_on_srs, 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, S: Sy // Let's open up C to every y. fn compute_value(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; - for &coeff in poly_positive { - let mut coeff = coeff; - coeff.mul_assign(&tmp); - value.add_assign(&coeff); - tmp.mul_assign(&y); - } + let positive_powers_contrib = evaluate_at_consequitive_powers(poly_positive, *y, *y); + let negative_powers_contrib = evaluate_at_consequitive_powers(poly_negative, yinv, yinv); + value.add_assign(&positive_powers_contrib); + value.add_assign(&negative_powers_contrib); + + // 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 } diff --git a/src/sonic/helped/poly.rs b/src/sonic/helped/poly.rs index 22ce7e7..58b046e 100644 --- a/src/sonic/helped/poly.rs +++ b/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 SxEval { 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 { - 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.add_assign(&evaluate_at_consequitive_powers(& self.v[..], tmp, tmp)); + let mut tmp = x.pow(&[(self.v.len()+1) as u64]); + acc.add_assign(&evaluate_at_consequitive_powers(& self.w[..], tmp, x)); + + // let mut tmp = x_inv; + // for mut u in self.u { + // u.mul_assign(&tmp); + // acc.add_assign(&u); + // 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 SyEval { 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 { - coeff.mul_assign(&tmp); - acc.add_assign(&coeff); - tmp.mul_assign(&yinv); - } + + let positive_powers_contrib = evaluate_at_consequitive_powers(& self.positive_coeffs[..], y, y); + let negative_powers_contrib = evaluate_at_consequitive_powers(& self.negative_coeffs[..], yinv, yinv); + acc.add_assign(&positive_powers_contrib); + 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 } diff --git a/src/sonic/helped/prover.rs b/src/sonic/helped/prover.rs index f3f4645..2624bfc 100644 --- a/src/sonic/helped/prover.rs +++ b/src/sonic/helped/prover.rs @@ -310,23 +310,6 @@ pub fn create_proof_on_srs, 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, 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 }) diff --git a/src/sonic/tests/mod.rs b/src/sonic/tests/mod.rs new file mode 100644 index 0000000..075a109 --- /dev/null +++ b/src/sonic/tests/mod.rs @@ -0,0 +1 @@ +mod sonics; \ No newline at end of file diff --git a/src/sonic/tests/sonics.rs b/src/sonic/tests/sonics.rs index 91d6ddd..114a061 100644 --- a/src/sonic/tests/sonics.rs +++ b/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( + 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, + xr: Option, + 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 for MiMCDemo<'a, E> { + fn synthesize>( + 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 bellman::sonic::sonic::AdaptorCircuit; - use bellman::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs}; - use bellman::sonic::helped::{MultiVerifier, get_circuit_parameters}; - use bellman::sonic::helped::helper::{create_aggregate_on_srs}; + use crate::sonic::cs::Basic; + use crate::sonic::sonic::AdaptorCircuit; + use crate::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs}; + use crate::sonic::helped::{MultiVerifier, get_circuit_parameters}; + 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 bellman::sonic::sonic::AdaptorCircuit; - use bellman::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs}; - use bellman::sonic::helped::{MultiVerifier, get_circuit_parameters}; - use bellman::sonic::helped::helper::{create_aggregate_on_srs}; + use crate::sonic::cs::Basic; + use crate::sonic::sonic::AdaptorCircuit; + use crate::sonic::helped::prover::{create_advice_on_srs, create_proof_on_srs}; + use crate::sonic::helped::{MultiVerifier, get_circuit_parameters}; + use crate::sonic::helped::helper::{create_aggregate_on_srs}; let info = get_circuit_parameters::(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, diff --git a/src/sonic/util.rs b/src/sonic/util.rs index 9f3e055..16505a8 100644 --- a/src/sonic/util.rs +++ b/src/sonic/util.rs @@ -123,10 +123,12 @@ pub fn polynomial_commitment_opening< ) -> E::G1Affine where I::IntoIter: DoubleEndedIterator + ExactSizeIterator, { - let poly = kate_divison( - polynomial_coefficients, - point, - ); + let poly = parallel_kate_divison::(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>(a: I, mut b: E::Fr) -> Vec +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 = a.into_iter().map(|el| el.clone()).collect(); + + let (q, _) = kate_divison_inner::(poly, divisor, reciproical, b_squared); + + // assert_eq!(r.len(), 0); + + q +} + +fn kate_divison_inner( + poly: Vec, + divisor: Vec, + reciproical: Vec, + remainder: E::Fr + ) -> (Vec, Vec) { + 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::(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::(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::(divisor, q.clone()); + sub_polynomials(&mut poly[..], &tmp[..]); + + return (q, poly); +} + /// Convenience function to check polynomail commitment pub fn check_polynomial_commitment( commitment: &E::G1Affine, @@ -524,6 +595,90 @@ pub fn multiply_polynomials(a: Vec, b: Vec) -> Vec(a: Vec, b: Vec) -> Vec { + 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> = a.into_iter().map(|e| Scalar::(e)).collect(); + let mut scalars_b: Vec> = b.into_iter().map(|e| Scalar::(e)).collect(); + scalars_a.resize(m, Scalar::(E::Fr::zero())); + scalars_b.resize(m, Scalar::(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 = scalars_a.into_iter().map(|e| e.0).collect(); + + mul_result.truncate(result_len); + + mul_result +} + pub fn multiply_polynomials_serial(mut a: Vec, mut b: Vec) -> Vec { let result_len = a.len() + b.len() - 1; @@ -574,6 +729,7 @@ pub fn multiply_polynomials_serial(mut a: Vec, mut b: Vec(a: &mut [F], b: &[F]) { use crate::multicore::Worker; use crate::domain::{EvaluationDomain, Scalar}; @@ -594,6 +750,28 @@ pub fn add_polynomials(a: &mut [F], b: &[F]) { }); } +// subtract polynomails in coefficient form +pub fn sub_polynomials(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(a: &mut [F], b: F) { use crate::multicore::Worker; use crate::domain::{EvaluationDomain, Scalar}; @@ -612,6 +790,8 @@ pub fn mul_polynomial_by_scalar(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(a: &mut [F], b: &[F], c: F) { use crate::multicore::Worker; use crate::domain::{EvaluationDomain, Scalar}; @@ -693,7 +873,6 @@ impl OptionExt for Option { } #[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::(&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::(&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::(&poly, z); + + assert_eq!(quotient_poly, parallel_q_poly); } \ No newline at end of file