diff --git a/Cargo.lock b/Cargo.lock index e4de0a4894..fc48089f77 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -7619,6 +7619,7 @@ dependencies = [ "aes-prng", "anyhow", "bc2wrap", + "criterion", "error-utils", "g2p", "itertools 0.14.0", diff --git a/core/threshold-algebra/Cargo.toml b/core/threshold-algebra/Cargo.toml index 06c1de5846..7cc55522b2 100644 --- a/core/threshold-algebra/Cargo.toml +++ b/core/threshold-algebra/Cargo.toml @@ -24,7 +24,12 @@ tfhe-versionable.workspace = true [dev-dependencies] aes-prng.workspace = true bc2wrap.workspace = true +criterion.workspace = true num-traits.workspace = true paste.workspace = true proptest.workspace = true -rstest.workspace = true \ No newline at end of file +rstest.workspace = true + +[[bench]] +name = "bivariate" +harness = false diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs new file mode 100644 index 0000000000..3c3a681b3b --- /dev/null +++ b/core/threshold-algebra/benches/bivariate.rs @@ -0,0 +1,74 @@ +use aes_prng::AesRng; +use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; +use rand::SeedableRng; +use std::hint::black_box; +use threshold_algebra::{ + bivariate::BivariatePoly, galois_rings::degree_4::ResiduePolyF4Z128, structure_traits::Sample, +}; + +const DEGREES: [usize; 3] = [1, 4, 13]; + +fn bivariate_setup(degree: usize) -> (BivariatePoly, ResiduePolyF4Z128) { + let mut rng = AesRng::seed_from_u64(degree as u64); + let secret = ResiduePolyF4Z128::sample(&mut rng); + let point = ResiduePolyF4Z128::sample(&mut rng); + let poly = BivariatePoly::from_secret(&mut rng, secret, degree); + (poly, point) +} + +fn bench_bivariate_sampling(c: &mut Criterion) { + let mut group = c.benchmark_group("bivariate/from_secret"); + + for degree in DEGREES { + group.bench_function(BenchmarkId::from_parameter(degree), |b| { + let mut rng = AesRng::seed_from_u64(degree as u64); + let secret = ResiduePolyF4Z128::sample(&mut rng); + b.iter(|| { + black_box(BivariatePoly::from_secret( + &mut rng, + black_box(secret), + black_box(degree), + )) + }); + }); + } + + group.finish(); +} + +fn bench_bivariate_evaluation(c: &mut Criterion) { + let mut group = c.benchmark_group("bivariate/evaluation"); + + for degree in DEGREES { + let (poly, point) = bivariate_setup(degree); + + group.bench_function(BenchmarkId::new("partial_x", degree), |b| { + b.iter(|| black_box(poly.partial_x_eval(black_box(point)))); + }); + group.bench_function(BenchmarkId::new("partial_y", degree), |b| { + b.iter(|| black_box(poly.partial_y_eval(black_box(point)))); + }); + // Baseline: two separate partial evals, mirroring the code before PR 576. + group.bench_function(BenchmarkId::new("partial_x_then_y", degree), |b| { + b.iter(|| { + let p = black_box(point); + ( + black_box(poly.partial_x_eval(p)), + black_box(poly.partial_y_eval(p)), + ) + }); + }); + group.bench_function(BenchmarkId::new("full", degree), |b| { + b.iter(|| black_box(poly.full_eval(black_box(point), black_box(point)))); + }); + } + + group.finish(); +} + +criterion_group!( + bivariate, + bench_bivariate_sampling, + bench_bivariate_evaluation, +); +criterion_main!(bivariate); diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 4aa3268352..5471657f2e 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -1,194 +1,115 @@ use super::poly::Poly; -use super::structure_traits::One; use super::structure_traits::Ring; use super::structure_traits::Sample; -use super::structure_traits::Zero; - -use anyhow::Result; -use error_utils::anyhow_error_and_log; -use itertools::Itertools; -use ndarray::Array; -use ndarray::ArrayD; -use ndarray::IxDyn; + use rand::{CryptoRng, Rng}; -use std::ops::Mul; -/// Bivariate polynomial is a matrix of coefficients of ResiduePolynomials -/// The row view of the polynomials is the following: -/// [[a_{00}, a_{01}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}]] -#[derive(Clone, Default, Debug)] +/// Bivariate polynomial represented as a square, row-major coefficient matrix. +/// +/// The row view of the coefficients is: +/// [[a_{00}, a_{01}, ..., a_{0t}], ..., [a_{t0}, ..., a_{tt}]] +/// +/// The `degree` field is the common degree in each variable, not the usual +/// total degree of a bivariate polynomial. This type assumes +/// `deg_X == deg_Y == degree`, so `coefs.len() == (degree + 1)^2`. +#[derive(Clone, Debug)] pub struct BivariatePoly { - pub coefs: ArrayD, + /// Row-major; `(degree + 1)^2` elements. + coefs: Vec, degree: usize, } impl BivariatePoly { - /// method for sampling random bivariate polynomial where free term is the secret - pub fn from_secret(rng: &mut R, secret: Z, degree: usize) -> Result + /// Samples a random bivariate polynomial where the free term is the secret. + pub fn from_secret(rng: &mut R, secret: Z, degree: usize) -> Self where - Z: Sample + Zero + Copy, + Z: Sample, { let d = degree + 1; - let coefs: Vec<_> = (0..d * d) - .map(|i| if i == 0 { secret } else { Z::sample(rng) }) - .collect(); - Ok(BivariatePoly { - coefs: ArrayD::from_shape_vec(IxDyn(&[d, d]), coefs)?.into_dyn(), - degree, - }) + let n = d * d; + let mut coefs = Vec::with_capacity(n); + coefs.push(secret); + coefs.extend((1..n).map(|_| Z::sample(rng))); + + BivariatePoly { coefs, degree } } -} -/// computes powers of a list of points in F up to a given maximal exponent -pub fn compute_powers_list + Copy>( - points: &[F], - max_exponent: usize, -) -> Vec> { - let mut alpha_powers = Vec::new(); - for p in points { - alpha_powers.push(compute_powers(*p, max_exponent)); + // Test-only constructor to allow for externally generated test-vectors (e.g. Sage). + #[cfg(test)] + fn from_coeffs(coefs: Vec, degree: usize) -> Self { + let d = degree + 1; + assert_eq!(coefs.len(), d * d); + BivariatePoly { coefs, degree } } - alpha_powers -} -/// Computes powers of a specific point up to degree: p^0, p^1,...,p^degree -pub fn compute_powers + Copy>(point: Z, degree: usize) -> Vec { - let mut powers_of_point = Vec::new(); - powers_of_point.push(Z::ONE); // start with - for i in 1..=degree { - powers_of_point.push(powers_of_point[i - 1] * point); + fn dim(&self) -> usize { + self.degree + 1 } - powers_of_point -} -pub trait MatrixMul { - fn matmul(&self, rhs: &ArrayD) -> Result>; + fn coeffs(&self) -> &[Z] { + &self.coefs + } } -impl MatrixMul for ArrayD { - fn matmul(&self, rhs: &ArrayD) -> Result> { - match (self.ndim(), rhs.ndim()) { - (1, 1) => { - if self.dim() != rhs.dim() { - return Err(anyhow_error_and_log(format!( - "Cannot compute multiplication between rank 1 tensor where dimension of lhs {:?} and rhs {:?}", - self.dim(), - rhs.dim() - ))); - } - if self.len() != rhs.len() { - return Err(anyhow_error_and_log(format!( - "Cannot multiply lhs of {:?} elements and rhs of {:?} elements for rank 1 tensors", - self.len(), - rhs.len() - ))); - } - let res = self - .iter() - .zip_eq(rhs) - .fold(Z::ZERO, |acc, (a, b)| acc + *a * *b); - Ok(Array::from_elem(IxDyn(&[1]), res).into_dyn()) +impl BivariatePoly { + /// Given a bivariate polynomial + /// F(X, Y) = \sum_{i=0}^{d-1} \sum_{j=0}^{d-1} a_{ij} X^i Y^j, + /// evaluate the X variable at \alpha. + /// + /// Returns G(Y) = F(\alpha, Y), whose Y^j coefficient is + /// \sum_{i=0}^{d-1} a_{ij} \alpha^i. + pub fn partial_x_eval(&self, alpha: Z) -> Poly { + let d = self.dim(); + let coefs = self.coeffs(); + let mut res = coefs[(d - 1) * d..d * d].to_vec(); + // Every row has exactly `d` coefficients. + for row in coefs[..(d - 1) * d].chunks_exact(d).rev() { + for (res_coef, coef) in res.iter_mut().zip(row) { + *res_coef = *res_coef * alpha + *coef; } - (1, 2) => { - if self.dim()[0] != rhs.dim()[0] { - Err(anyhow_error_and_log(format!( - "Cannot compute multiplication between rank 1 tensor and rank 2 tensor where dimension of lhs {:?} and rhs {:?}", - self.dim(), - rhs.dim() - ))) - } else { - let mut res = Vec::new(); - for col in rhs.columns() { - if col.len() != self.len() { - return Err(anyhow_error_and_log(format!( - "Cannot multiply lhs of {:?} elements and rhs of {:?} elements for rank 1 tensors and rank 2 tensors", - self.len(), - rhs.len() - ))); - } - let s = col - .iter() - .zip_eq(self) - .fold(Z::ZERO, |acc, (a, b)| acc + *b * *a); - res.push(s); - } - Ok(Array::from_vec(res).into_dyn()) - } - } - (2, 1) => { - if self.dim()[1] != rhs.dim()[0] { - Err(anyhow_error_and_log(format!( - "Cannot compute multiplication between rank 2 tensor and rank 1 tensor where dimension of lhs {:?} and rhs {:?}", - self.dim(), - rhs.dim() - ))) - } else { - let mut res = Vec::new(); - for row in self.rows() { - if row.len() != rhs.len() { - return Err(anyhow_error_and_log(format!( - "Cannot multiply lhs of {:?} elements and rhs of {:?} elements for rank 2 tensors and rank 1 tensors", - self.len(), - rhs.len() - ))); - } - let s = row - .iter() - .zip_eq(rhs) - .fold(Z::ZERO, |acc, (a, b)| acc + *b * *a); - res.push(s); - } - Ok(Array::from_vec(res).into_dyn()) - } - } - (l_rank, r_rank) => Err(anyhow_error_and_log(format!( - "Matmul not implemented for tensors of rank {l_rank:?}, {r_rank:?}", - ))), } - } -} - -pub trait BivariateEval { - /// Given a degree T bivariate poly F(X,Y) and a point \alpha, we compute - /// G(X) = F(X, \alpha) as - /// [\alpha^0, ..., \alpha_d].matmul([a_{00}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}] = - /// [sum(alpha^j * a_{j0}), ..., sum(alpha^j * a_{jd})] - fn partial_x_evaluation(&self, alpha: Z) -> Result>; - - /// Given a degree T bivariate poly F(X,Y) and a point \alpha, we compute - /// G(Y) := F(\alpha, Y) as - /// [a_{00}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}].matmul([\alpha^0, ..., \alpha_d]) - /// [sum(alpha^j * a_{0j}), ..., sum(alpha^j * a_{dj})] - fn partial_y_evaluation(&self, alpha: Z) -> Result>; - - /// Given a degree T bivariate poly F(X,Y) and two points \alpha_x, \alpha_y, we compute - /// F(\alpha_x, \alpha_y) - fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result; -} - -impl BivariateEval for BivariatePoly -where - ArrayD: MatrixMul, -{ - fn partial_x_evaluation(&self, alpha: Z) -> Result> { - let powers_array = Array::from(compute_powers(alpha, self.degree)).into_dyn(); - let res_vector = powers_array.matmul(&self.coefs)?; - Ok(Poly::from_coefs(res_vector.into_raw_vec_and_offset().0)) + Poly::from_coefs(res) } - fn partial_y_evaluation(&self, alpha: Z) -> Result> { - let powers_array = Array::from(compute_powers(alpha, self.degree)).into_dyn(); - let res_vector = self.coefs.matmul(&powers_array)?; - Ok(Poly::from_coefs(res_vector.into_raw_vec_and_offset().0)) + /// Given a bivariate polynomial + /// F(X, Y) = \sum_{i=0}^{d-1} \sum_{j=0}^{d-1} a_{ij} X^i Y^j, + /// evaluate the Y variable at \alpha. + /// + /// Returns G(X) = F(X, \alpha), whose X^i coefficient is + /// \sum_{j=0}^{d-1} a_{ij} \alpha^j. + pub fn partial_y_eval(&self, alpha: Z) -> Poly { + let d = self.dim(); + let coefs = self.coeffs(); + let mut res = Vec::with_capacity(d); + // Each pass consumes exactly one row of length `d`. + for row in coefs.chunks_exact(d) { + let mut acc = row[d - 1]; + for coef in row[..d - 1].iter().rev() { + acc = acc * alpha + *coef; + } + res.push(acc); + } + Poly::from_coefs(res) } - fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result { - let powers_array_x = Array::from(compute_powers(alpha_x, self.degree)).into_dyn(); - let powers_array_y = Array::from(compute_powers(alpha_y, self.degree)).into_dyn(); - - let lhs = powers_array_x.matmul(&self.coefs)?; - let res = lhs.matmul(&powers_array_y)?; - Ok(res[0]) + /// Given a bivariate poly F(X,Y) and two points \alpha_x, \alpha_y, compute F(\alpha_x, \alpha_y) + pub fn full_eval(&self, alpha_x: Z, alpha_y: Z) -> Z { + let d = self.dim(); + let coefs = self.coeffs(); + let mut rows = coefs.chunks_exact(d).rev(); + let last = rows.next().expect("d >= 1"); + let mut acc = last[d - 1]; + for coef in last[..d - 1].iter().rev() { + acc = acc * alpha_y + *coef; + } + for row in rows { + let mut row_acc = row[d - 1]; + for coef in row[..d - 1].iter().rev() { + row_acc = row_acc * alpha_y + *coef; + } + acc = acc * alpha_x + row_acc; + } + acc } } @@ -196,12 +117,13 @@ where mod tests { use super::*; use crate::galois_rings::degree_8::ResiduePolyF8Z128; + use crate::structure_traits::{FromU128, ZConsts}; use crate::{ galois_rings::{ common::ResiduePoly, degree_4::{ResiduePolyF4Z64, ResiduePolyF4Z128}, }, - structure_traits::One, + structure_traits::{One, Zero}, }; use aes_prng::AesRng; @@ -209,91 +131,107 @@ mod tests { use rstest::rstest; use std::num::Wrapping; - //Checks that we error on incompatible sizes and dimensions + //Checks the hot coefficient shapes used by bivariate evaluation. #[test] - fn test_matmul_bounds() { - let x11 = ArrayD::from_elem(IxDyn(&[1, 1]), ResiduePolyF4Z128::ONE); - let y2 = ArrayD::from_elem(IxDyn(&[2]), ResiduePolyF4Z128::ONE); - // test (1, 1) X (2) mul error - assert!(x11.matmul(&y2).is_err()); - assert!(y2.matmul(&x11).is_err()); - - // we do not support mul between two 2d matrices - assert!(x11.matmul(&x11).is_err()); - - let z22 = ArrayD::from_elem(IxDyn(&[2, 2]), ResiduePolyF4Z128::ONE); - // test vec-matrix bound check returns ok - assert!(y2.matmul(&z22).is_ok()); - - // test matrix-vec bound check returns ok - assert!(z22.matmul(&y2).is_ok()); - - let y4 = ArrayD::from_elem(IxDyn(&[4]), ResiduePolyF4Z128::ONE); + fn test_bivariate_supported_shapes() { + let two = ResiduePolyF4Z128::TWO; + let bpoly2 = BivariatePoly::from_coeffs(vec![ResiduePolyF4Z128::ONE; 4], 1); + let expected2 = Poly::from_coefs(vec![two; 2]); + assert_eq!(bpoly2.partial_x_eval(ResiduePolyF4Z128::ONE), expected2); + assert_eq!(bpoly2.partial_y_eval(ResiduePolyF4Z128::ONE), expected2); + assert_eq!( + bpoly2.full_eval(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), + two + two + ); + + let five = ResiduePolyF4Z128::from_u128(5); + let bpoly5 = BivariatePoly::from_coeffs(vec![ResiduePolyF4Z128::ONE; 25], 4); + let expected5 = Poly::from_coefs(vec![five; 5]); + assert_eq!(bpoly5.partial_x_eval(ResiduePolyF4Z128::ONE), expected5); + assert_eq!(bpoly5.partial_y_eval(ResiduePolyF4Z128::ONE), expected5); + assert_eq!( + bpoly5.full_eval(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), + ResiduePolyF4Z128::from_u128(25) + ); + } - // test 1x1 vector mul errors - assert!(y4.matmul(&y2).is_err()); - assert!(y2.matmul(&y4).is_err()); + #[test] + fn bivariate_partial_evals_match_full() { + // F(X, Y) = 1 + 2*Y + 3*X + 4*X*Y, stored row-major as [a_00, a_01, a_10, a_11]. + let one = ResiduePolyF4Z128::ONE; + let two = ResiduePolyF4Z128::TWO; + let three = ResiduePolyF4Z128::THREE; + let four = ResiduePolyF4Z128::from_u128(4); + let five = ResiduePolyF4Z128::from_u128(5); + let seven = ResiduePolyF4Z128::from_u128(7); + let bpoly = BivariatePoly::from_coeffs(vec![one, two, three, four], 1); + + // F(5, Y) = (1 + 3*5) + (2 + 4*5)*Y = 16 + 22*Y + let sixteen = ResiduePolyF4Z128::from_u128(16); + let twentytwo = ResiduePolyF4Z128::from_u128(22); + assert_eq!( + bpoly.partial_x_eval(five), + Poly::from_coefs(vec![sixteen, twentytwo]) + ); + + // F(X, 5) = (1 + 2*5) + (3 + 4*5)*X = 11 + 23*X + let eleven = ResiduePolyF4Z128::from_u128(11); + let twentythree = ResiduePolyF4Z128::from_u128(23); + assert_eq!( + bpoly.partial_y_eval(five), + Poly::from_coefs(vec![eleven, twentythree]) + ); + + // F(5, 7) = 1 + 2*7 + 3*5 + 4*5*7 = 170 + let expected_full = ResiduePolyF4Z128::from_u128(170); + assert_eq!(bpoly.full_eval(five, seven), expected_full); } //Test that eval at 0 return the secret for ResiduePolyF4Z128 #[rstest] - #[case(4)] - #[case(10)] - #[case(20)] - fn test_bivariate_zero_128(#[case] degree: usize) { + fn test_bivariate_zero_128(#[values(4, 10, 20)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z128::sample(&mut rng); - let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree).unwrap(); - let ev_zero = bpoly - .full_evaluation(ResiduePolyF4Z128::ZERO, ResiduePolyF4Z128::ZERO) - .unwrap(); + let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); + let ev_zero = bpoly.full_eval(ResiduePolyF4Z128::ZERO, ResiduePolyF4Z128::ZERO); assert_eq!(ev_zero, secret); } //Test that eval at 0 return the secret for ResiduePolyF4Z64 #[rstest] - #[case(4)] - #[case(10)] - #[case(20)] - fn test_bivariate_zero_64(#[case] degree: usize) { + fn test_bivariate_zero_64(#[values(4, 10, 20)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z64::sample(&mut rng); - let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree).unwrap(); - let ev_zero = bpoly - .full_evaluation(ResiduePolyF4Z64::ZERO, ResiduePolyF4Z64::ZERO) - .unwrap(); + let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); + let ev_zero = bpoly.full_eval(ResiduePolyF4Z64::ZERO, ResiduePolyF4Z64::ZERO); assert_eq!(ev_zero, secret); } //Test that eval at 1 return the sum of all coefs of the poly for ResiduePolyF4Z128 #[rstest] - #[case(4)] - #[case(10)] - #[case(20)] - fn test_bivariate_one_128(#[case] degree: usize) { + fn test_bivariate_one_128(#[values(4, 10, 20)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z128::sample(&mut rng); - let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree).unwrap(); - let ev_one = bpoly - .full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE) - .unwrap(); - let sum_coefs = bpoly.coefs.iter().fold(ResiduePoly::ZERO, |acc, x| acc + x); + let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); + let ev_one = bpoly.full_eval(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE); + let sum_coefs = bpoly + .coeffs() + .iter() + .fold(ResiduePoly::ZERO, |acc, x| acc + x); assert_eq!(ev_one, sum_coefs); } //Test that eval at 1 return the sum of all coefs of the poly for ResiduePolyF4Z64 #[rstest] - #[case(4)] - #[case(10)] - #[case(20)] - fn test_bivariate_one_64(#[case] degree: usize) { + fn test_bivariate_one_64(#[values(4, 10, 20)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z64::sample(&mut rng); - let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree).unwrap(); - let ev_one = bpoly - .full_evaluation(ResiduePolyF4Z64::ONE, ResiduePolyF4Z64::ONE) - .unwrap(); - let sum_coefs = bpoly.coefs.iter().fold(ResiduePoly::ZERO, |acc, x| acc + x); + let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); + let ev_one = bpoly.full_eval(ResiduePolyF4Z64::ONE, ResiduePolyF4Z64::ONE); + let sum_coefs = bpoly + .coeffs() + .iter() + .fold(ResiduePoly::ZERO, |acc, x| acc + x); assert_eq!(ev_one, sum_coefs); } @@ -602,12 +540,7 @@ mod tests { }, //x4y4 ]; - let bpoly = BivariatePoly { - coefs: ArrayD::from_shape_vec(IxDyn(&[5, 5]), coefs) - .unwrap() - .to_owned(), - degree: 4, - }; + let bpoly = BivariatePoly::from_coeffs(coefs, 4); let point = ResiduePoly { coefs: [ @@ -629,7 +562,7 @@ mod tests { #[test] fn test_bivariate_partial_eval_x() { let (bpoly, point) = poly_setup(); - let res = bpoly.partial_x_evaluation(point).unwrap(); + let res = bpoly.partial_x_eval(point); let expected_result = Poly::::from_coefs(vec![ ResiduePoly { @@ -702,7 +635,7 @@ mod tests { fn test_bivariate_partial_eval_y() { //Taking Sage as reference let (bpoly, point) = poly_setup(); - let res = bpoly.partial_y_evaluation(point).unwrap(); + let res = bpoly.partial_y_eval(point); let expected_result = Poly::::from_coefs(vec![ ResiduePoly { @@ -774,9 +707,9 @@ mod tests { let (bpoly, point) = poly_setup(); let point_x = point; let point_y = point_x + point_x; - let res = bpoly.full_evaluation(point_x, point_y).unwrap(); + let res = bpoly.full_eval(point_x, point_y); - let expected_res = bpoly.partial_x_evaluation(point_x).unwrap().eval(&point_y); + let expected_res = bpoly.partial_x_eval(point_x).eval(&point_y); assert_eq!(res, expected_res); } diff --git a/core/threshold-algebra/src/galois_rings/common.rs b/core/threshold-algebra/src/galois_rings/common.rs index 4f893428d5..8bd79fdc9d 100644 --- a/core/threshold-algebra/src/galois_rings/common.rs +++ b/core/threshold-algebra/src/galois_rings/common.rs @@ -1,7 +1,7 @@ use super::utils::ArrayVisitor; use crate::{ base_ring::ToZ64, - bivariate::compute_powers_list, + matrix::compute_powers_list, poly::Poly, structure_traits::{ BaseRing, Derive, FromU128, Invert, One, QuotientMaximalIdeal, Ring, diff --git a/core/threshold-algebra/src/galois_rings/degree_3.rs b/core/threshold-algebra/src/galois_rings/degree_3.rs index 66aadd216f..5d945c6cb2 100644 --- a/core/threshold-algebra/src/galois_rings/degree_3.rs +++ b/core/threshold-algebra/src/galois_rings/degree_3.rs @@ -11,9 +11,9 @@ use std::{ use crate::galois_rings::ExceptionalSetMap; use crate::{ base_ring::{Z64, Z128}, - bivariate::compute_powers, error_correction::MemoizedExceptionals, galois_fields::gf8::{GF8, GF8_FROM_GENERATOR}, + matrix::compute_powers, poly::{BitWiseEval, BitwisePoly}, structure_traits::{ BaseRing, One, QuotientMaximalIdeal, Ring, RingWithExceptionalSequence, Solve1, ZConsts, diff --git a/core/threshold-algebra/src/galois_rings/degree_4.rs b/core/threshold-algebra/src/galois_rings/degree_4.rs index ce20c16e5e..c09b754e81 100644 --- a/core/threshold-algebra/src/galois_rings/degree_4.rs +++ b/core/threshold-algebra/src/galois_rings/degree_4.rs @@ -11,9 +11,9 @@ use std::{ use crate::galois_rings::ExceptionalSetMap; use crate::{ base_ring::{Z64, Z128}, - bivariate::compute_powers, error_correction::MemoizedExceptionals, galois_fields::gf16::{GF16, GF16_FROM_GENERATOR, GF16_NEWTON_INNER_LOOP, two_powers}, + matrix::compute_powers, poly::{BitWiseEval, BitwisePoly}, structure_traits::{ BaseRing, One, QuotientMaximalIdeal, Ring, RingWithExceptionalSequence, Solve1, ZConsts, diff --git a/core/threshold-algebra/src/galois_rings/degree_5.rs b/core/threshold-algebra/src/galois_rings/degree_5.rs index a7d0c38580..09107916dc 100644 --- a/core/threshold-algebra/src/galois_rings/degree_5.rs +++ b/core/threshold-algebra/src/galois_rings/degree_5.rs @@ -11,9 +11,9 @@ use std::{ use crate::galois_rings::ExceptionalSetMap; use crate::{ base_ring::{Z64, Z128}, - bivariate::compute_powers, error_correction::MemoizedExceptionals, galois_fields::gf32::{GF32, GF32_FROM_GENERATOR}, + matrix::compute_powers, poly::{BitWiseEval, BitwisePoly}, structure_traits::{ BaseRing, One, QuotientMaximalIdeal, Ring, RingWithExceptionalSequence, Solve1, ZConsts, diff --git a/core/threshold-algebra/src/galois_rings/degree_6.rs b/core/threshold-algebra/src/galois_rings/degree_6.rs index 33fe12acd8..f7142445d8 100644 --- a/core/threshold-algebra/src/galois_rings/degree_6.rs +++ b/core/threshold-algebra/src/galois_rings/degree_6.rs @@ -11,9 +11,9 @@ use std::{ use crate::galois_rings::ExceptionalSetMap; use crate::{ base_ring::{Z64, Z128}, - bivariate::compute_powers, error_correction::MemoizedExceptionals, galois_fields::gf64::{GF64, GF64_FROM_GENERATOR, GF64_NEWTON_INNER_LOOP, two_powers}, + matrix::compute_powers, poly::{BitWiseEval, BitwisePoly}, structure_traits::{ BaseRing, One, QuotientMaximalIdeal, Ring, RingWithExceptionalSequence, Solve1, ZConsts, diff --git a/core/threshold-algebra/src/galois_rings/degree_7.rs b/core/threshold-algebra/src/galois_rings/degree_7.rs index 1725c4ccdb..55db561f4a 100644 --- a/core/threshold-algebra/src/galois_rings/degree_7.rs +++ b/core/threshold-algebra/src/galois_rings/degree_7.rs @@ -11,9 +11,9 @@ use std::{ use crate::galois_rings::ExceptionalSetMap; use crate::{ base_ring::{Z64, Z128}, - bivariate::compute_powers, error_correction::MemoizedExceptionals, galois_fields::gf128::{GF128, GF128_FROM_GENERATOR}, + matrix::compute_powers, poly::{BitWiseEval, BitwisePoly}, structure_traits::{ BaseRing, One, QuotientMaximalIdeal, Ring, RingWithExceptionalSequence, Solve1, ZConsts, diff --git a/core/threshold-algebra/src/galois_rings/degree_8.rs b/core/threshold-algebra/src/galois_rings/degree_8.rs index fca9802bfc..149b510a77 100644 --- a/core/threshold-algebra/src/galois_rings/degree_8.rs +++ b/core/threshold-algebra/src/galois_rings/degree_8.rs @@ -11,9 +11,9 @@ use std::{ use crate::galois_rings::ExceptionalSetMap; use crate::{ base_ring::{Z64, Z128}, - bivariate::compute_powers, error_correction::MemoizedExceptionals, galois_fields::gf256::{GF256, GF256_FROM_GENERATOR, GF256_NEWTON_INNER_LOOP, two_powers}, + matrix::compute_powers, poly::{BitWiseEval, BitwisePoly}, structure_traits::{ BaseRing, One, QuotientMaximalIdeal, Ring, RingWithExceptionalSequence, Solve1, ZConsts, diff --git a/core/threshold-algebra/src/lib.rs b/core/threshold-algebra/src/lib.rs index b1e3d8baa9..027e941624 100644 --- a/core/threshold-algebra/src/lib.rs +++ b/core/threshold-algebra/src/lib.rs @@ -4,6 +4,7 @@ pub mod commitment; pub mod error_correction; pub mod galois_fields; pub mod galois_rings; +pub mod matrix; pub mod poly; pub mod randomness_check; pub mod sharing; diff --git a/core/threshold-algebra/src/matrix.rs b/core/threshold-algebra/src/matrix.rs new file mode 100644 index 0000000000..bef220e234 --- /dev/null +++ b/core/threshold-algebra/src/matrix.rs @@ -0,0 +1,205 @@ +use crate::structure_traits::{One, Ring}; + +use anyhow::Result; +use error_utils::anyhow_error_and_log; +use itertools::Itertools; +use ndarray::{Array, ArrayD, IxDyn}; +use std::ops::Mul; + +pub trait MatrixMul: Sized { + type Output; + fn matmul(&self, rhs: &Rhs) -> Result; +} + +/// Computes powers of a specific point up to degree: p^0, p^1,...,p^degree +pub fn compute_powers + Copy>(point: Z, degree: usize) -> Vec { + let mut powers_of_point = Vec::with_capacity(degree + 1); + powers_of_point.push(Z::ONE); + for i in 1..=degree { + powers_of_point.push(powers_of_point[i - 1] * point); + } + powers_of_point +} + +/// Computes powers of a list of points in F up to a given maximal exponent. +pub fn compute_powers_list + Copy>( + points: &[F], + max_exponent: usize, +) -> Vec> { + let mut alpha_powers = Vec::with_capacity(points.len()); + for p in points { + alpha_powers.push(compute_powers(*p, max_exponent)); + } + alpha_powers +} + +impl MatrixMul> for ArrayD { + type Output = ArrayD; + + fn matmul(&self, rhs: &ArrayD) -> Result { + match (self.ndim(), rhs.ndim()) { + (1, 1) => { + if self.dim() != rhs.dim() { + return Err(anyhow_error_and_log(format!( + "Cannot compute multiplication between rank 1 tensor where dimension of lhs {:?} and rhs {:?}", + self.dim(), + rhs.dim() + ))); + } + if self.len() != rhs.len() { + return Err(anyhow_error_and_log(format!( + "Cannot multiply lhs of {:?} elements and rhs of {:?} elements for rank 1 tensors", + self.len(), + rhs.len() + ))); + } + let res = self + .iter() + .zip_eq(rhs) + .fold(Z::ZERO, |acc, (a, b)| acc + *a * *b); + Ok(Array::from_elem(IxDyn(&[1]), res).into_dyn()) + } + (1, 2) => { + if self.dim()[0] != rhs.dim()[0] { + Err(anyhow_error_and_log(format!( + "Cannot compute multiplication between rank 1 tensor and rank 2 tensor where dimension of lhs {:?} and rhs {:?}", + self.dim(), + rhs.dim() + ))) + } else { + let mut res = Vec::with_capacity(rhs.shape()[1]); + for col in rhs.columns() { + if col.len() != self.len() { + return Err(anyhow_error_and_log(format!( + "Cannot multiply lhs of {:?} elements and rhs of {:?} elements for rank 1 tensors and rank 2 tensors", + self.len(), + rhs.len() + ))); + } + let s = col + .iter() + .zip_eq(self) + .fold(Z::ZERO, |acc, (a, b)| acc + *b * *a); + res.push(s); + } + Ok(Array::from_vec(res).into_dyn()) + } + } + (2, 1) => { + if self.dim()[1] != rhs.dim()[0] { + Err(anyhow_error_and_log(format!( + "Cannot compute multiplication between rank 2 tensor and rank 1 tensor where dimension of lhs {:?} and rhs {:?}", + self.dim(), + rhs.dim() + ))) + } else { + let mut res = Vec::with_capacity(self.shape()[0]); + for row in self.rows() { + if row.len() != rhs.len() { + return Err(anyhow_error_and_log(format!( + "Cannot multiply lhs of {:?} elements and rhs of {:?} elements for rank 2 tensors and rank 1 tensors", + self.len(), + rhs.len() + ))); + } + let s = row + .iter() + .zip_eq(rhs) + .fold(Z::ZERO, |acc, (a, b)| acc + *b * *a); + res.push(s); + } + Ok(Array::from_vec(res).into_dyn()) + } + } + (l_rank, r_rank) => Err(anyhow_error_and_log(format!( + "Matmul not implemented for tensors of rank {l_rank:?}, {r_rank:?}", + ))), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::galois_rings::degree_4::ResiduePolyF4Z128; + use crate::structure_traits::{One, Zero}; + + fn rp(n: u8) -> ResiduePolyF4Z128 { + let mut acc = ResiduePolyF4Z128::ZERO; + for _ in 0..n { + acc += ResiduePolyF4Z128::ONE; + } + acc + } + + #[test] + fn compute_powers_basic() { + let powers = compute_powers(rp(2), 4); + assert_eq!(powers, vec![rp(1), rp(2), rp(4), rp(8), rp(16)]); + } + + #[test] + fn compute_powers_zero_degree() { + let powers = compute_powers(rp(7), 0); + assert_eq!(powers, vec![rp(1)]); + } + + #[test] + fn compute_powers_list_basic() { + let lists = compute_powers_list(&[rp(2), rp(3)], 2); + assert_eq!(lists.len(), 2); + assert_eq!(lists[0], vec![rp(1), rp(2), rp(4)]); + assert_eq!(lists[1], vec![rp(1), rp(3), rp(9)]); + } + + #[test] + fn matmul_rank1_dot() { + let a = ArrayD::from_shape_vec(IxDyn(&[3]), vec![rp(1), rp(2), rp(3)]).unwrap(); + let b = ArrayD::from_shape_vec(IxDyn(&[3]), vec![rp(4), rp(5), rp(6)]).unwrap(); + let res = a.matmul(&b).unwrap(); + // 1*4 + 2*5 + 3*6 = 32 + assert_eq!(res.into_raw_vec_and_offset().0, vec![rp(32)]); + } + + #[test] + fn matmul_vector_matrix() { + // [1,2] * [[1,2,3],[4,5,6]] = [9, 12, 15] + let v = ArrayD::from_shape_vec(IxDyn(&[2]), vec![rp(1), rp(2)]).unwrap(); + let m = ArrayD::from_shape_vec( + IxDyn(&[2, 3]), + vec![rp(1), rp(2), rp(3), rp(4), rp(5), rp(6)], + ) + .unwrap(); + let res = v.matmul(&m).unwrap(); + assert_eq!(res.into_raw_vec_and_offset().0, vec![rp(9), rp(12), rp(15)]); + } + + #[test] + fn matmul_matrix_vector() { + // [[1,2,3],[4,5,6]] * [7,8,9] = [50, 122] + let m = ArrayD::from_shape_vec( + IxDyn(&[2, 3]), + vec![rp(1), rp(2), rp(3), rp(4), rp(5), rp(6)], + ) + .unwrap(); + let v = ArrayD::from_shape_vec(IxDyn(&[3]), vec![rp(7), rp(8), rp(9)]).unwrap(); + let res = m.matmul(&v).unwrap(); + assert_eq!(res.into_raw_vec_and_offset().0, vec![rp(50), rp(122)]); + } + + #[test] + fn matmul_dim_mismatches_err() { + let v3 = ArrayD::from_shape_vec(IxDyn(&[3]), vec![rp(1); 3]).unwrap(); + let v2 = ArrayD::from_shape_vec(IxDyn(&[2]), vec![rp(1); 2]).unwrap(); + // rank-1 length mismatch + assert!(v3.matmul(&v2).is_err()); + // rank mismatch on rank-2 sides + let m22 = ArrayD::from_shape_vec(IxDyn(&[2, 2]), vec![rp(1); 4]).unwrap(); + assert!(m22.matmul(&m22).is_err()); + // (1, 2) shape mismatch + let m23 = ArrayD::from_shape_vec(IxDyn(&[2, 3]), vec![rp(1); 6]).unwrap(); + assert!(v3.matmul(&m23).is_err()); + // (2, 1) shape mismatch + assert!(m23.matmul(&v2).is_err()); + } +} diff --git a/core/threshold-algebra/src/syndrome.rs b/core/threshold-algebra/src/syndrome.rs index 3f284a312f..41aa8884b3 100644 --- a/core/threshold-algebra/src/syndrome.rs +++ b/core/threshold-algebra/src/syndrome.rs @@ -1,5 +1,5 @@ use super::{ - bivariate::compute_powers_list, + matrix::compute_powers_list, poly::Poly, sharing::shamir::ShamirFieldPoly, structure_traits::{Field, Ring}, diff --git a/core/threshold-execution/src/large_execution/double_sharing.rs b/core/threshold-execution/src/large_execution/double_sharing.rs index 8c14ebc89f..fe3d6012a2 100644 --- a/core/threshold-execution/src/large_execution/double_sharing.rs +++ b/core/threshold-execution/src/large_execution/double_sharing.rs @@ -4,7 +4,7 @@ use super::{ }; use crate::runtime::sessions::large_session::LargeSessionHandles; use algebra::{ - bivariate::MatrixMul, + matrix::MatrixMul, structure_traits::{Derive, ErrorCorrect, Invert, Ring}, }; use async_trait::async_trait; diff --git a/core/threshold-execution/src/large_execution/single_sharing.rs b/core/threshold-execution/src/large_execution/single_sharing.rs index 1797880f1e..cf6b033be8 100644 --- a/core/threshold-execution/src/large_execution/single_sharing.rs +++ b/core/threshold-execution/src/large_execution/single_sharing.rs @@ -1,7 +1,7 @@ use super::local_single_share::{LocalSingleShare, SecureLocalSingleShare}; use crate::runtime::sessions::large_session::LargeSessionHandles; use algebra::{ - bivariate::{MatrixMul, compute_powers}, + matrix::MatrixMul, structure_traits::{Derive, ErrorCorrect, Invert, Ring, RingWithExceptionalSequence}, }; use async_trait::async_trait; @@ -142,11 +142,14 @@ pub fn init_vdm( .map(|idx| Z::get_from_exceptional_sequence(idx + 1)) .try_collect()?; - let powers_of_exceptional_sequence: Vec = exceptional_sequence - .into_iter() - .fold(Vec::::new(), |acc, point| { - [acc, compute_powers(point, width - 1)].concat() - }); + let mut powers_of_exceptional_sequence = Vec::with_capacity(height * width); + for point in exceptional_sequence { + let mut power = Z::ONE; + for _ in 0..width { + powers_of_exceptional_sequence.push(power); + power *= point; + } + } Ok(ArrayD::from_shape_vec(IxDyn(&[height, width]), powers_of_exceptional_sequence)?.into_dyn()) } diff --git a/core/threshold-execution/src/large_execution/vss.rs b/core/threshold-execution/src/large_execution/vss.rs index 19d6c38e48..7b686ea3fd 100644 --- a/core/threshold-execution/src/large_execution/vss.rs +++ b/core/threshold-execution/src/large_execution/vss.rs @@ -15,7 +15,7 @@ use crate::{ }, }; use algebra::{ - bivariate::{BivariateEval, BivariatePoly}, + bivariate::BivariatePoly, poly::Poly, structure_traits::{Ring, RingWithExceptionalSequence}, }; @@ -102,10 +102,21 @@ pub(crate) struct DoublePoly where Poly: Eq, { + /// The share polynomial in variable X, F(X, \alpha) pub(crate) share_in_x: Poly, + /// The share polynomial in variable Y, F(\alpha, Y) pub(crate) share_in_y: Poly, } +impl DoublePoly { + pub(crate) fn from_bivariate(poly: &BivariatePoly, point: Z) -> Self { + Self { + share_in_x: poly.partial_y_eval(point), // F(X, alpha) + share_in_y: poly.partial_x_eval(point), // F(alpha, Y) + } + } +} + /// Struct to hold data sent during round 1 of VSS, composed of /// - double_poly is my share in a single VSS instance /// - we need n challenges sent and n challenges received (one from every party) @@ -262,10 +273,10 @@ pub(crate) fn sample_secret_polys anyhow::Result<(Vec>, MapRoleDoublePoly)> { let degree = session.threshold() as usize; //Sample the bivariate polynomials Vec - let bivariate_poly = secrets + let bivariate_poly: Vec<_> = secrets .iter() .map(|s| BivariatePoly::from_secret(session.rng(), *s, degree)) - .collect::, _>>()?; + .collect(); //Evaluate the bivariate poly in its first and second variables //to create a mapping role -> Vec<(F(X,alpha_role), F(alpha_role,Y))> let map_double_shares: MapRoleDoublePoly = session @@ -275,12 +286,7 @@ pub(crate) fn sample_secret_polys>), anyhow::Error>((*r, vec_map)) }) @@ -906,8 +912,8 @@ where (*own_role, *pi_role, *pj_role), vss.my_poly .iter() - .map(|poly| poly.full_evaluation(point_x, point_y)) - .collect::, _>>()?, + .map(|poly| poly.full_eval(point_x, point_y)) + .collect(), ); } //If im a Pi send Fi(alpha_j) @@ -1030,8 +1036,8 @@ fn round_4_conflict_resolution( true => ValueOrPoly::Poly( vss.my_poly .iter() - .map(|poly| poly.partial_y_evaluation(point_pi)) - .collect::, _>>()?, + .map(|poly| poly.partial_y_eval(point_pi)) + .collect(), ), //As P_j external from the conflict, resolve conflict with P_i by sending F(alpha_j,alpha_i) false => ValueOrPoly::Value( @@ -1177,7 +1183,6 @@ pub(crate) mod tests { }; use crate::tests::helper::tests_and_benches::execute_protocol_small; use algebra::{ - bivariate::BivariateEval, galois_rings::degree_4::{ResiduePolyF4, ResiduePolyF4Z64, ResiduePolyF4Z128}, sharing::{ shamir::{RevealOp, ShamirSharings}, @@ -1332,7 +1337,7 @@ pub(crate) mod tests { &result .my_poly .iter() - .map(|p| p.full_evaluation(x_0, y_0).unwrap()) + .map(|p| p.full_eval(x_0, y_0)) .collect_vec(), expected_secret, ); @@ -1344,12 +1349,12 @@ pub(crate) mod tests { let expected_result_x = result .my_poly .iter() - .map(|p| p.partial_y_evaluation(embedded_pn).unwrap()) + .map(|p| p.partial_y_eval(embedded_pn)) .collect_vec(); let expected_result_y = result .my_poly .iter() - .map(|p| p.partial_x_evaluation(embedded_pn).unwrap()) + .map(|p| p.partial_x_eval(embedded_pn)) .collect_vec(); assert_eq!( diff --git a/core/threshold-execution/src/malicious_execution/large_execution/malicious_vss.rs b/core/threshold-execution/src/malicious_execution/large_execution/malicious_vss.rs index cd4a3902dd..7a918e2124 100644 --- a/core/threshold-execution/src/malicious_execution/large_execution/malicious_vss.rs +++ b/core/threshold-execution/src/malicious_execution/large_execution/malicious_vss.rs @@ -13,7 +13,7 @@ use crate::{ runtime::sessions::base_session::BaseSessionHandles, }; use algebra::{ - bivariate::{BivariateEval, BivariatePoly}, + bivariate::BivariatePoly, poly::Poly, structure_traits::{Ring, RingWithExceptionalSequence}, }; @@ -197,9 +197,7 @@ async fn malicious_round_1 = session .roles() @@ -207,14 +205,7 @@ async fn malicious_round_1 WrongDegreeSharingVss { // to an honest behaviour. let degree = session.threshold() as usize + 1; //Sample the bivariate polynomials Vec - let bivariate_poly = secrets + let bivariate_poly: Vec<_> = secrets .iter() .map(|s| BivariatePoly::from_secret(session.rng(), *s, degree)) - .collect::, _>>()?; + .collect(); //Evaluate the bivariate poly in its first and second variables //to create a mapping role -> Vec<(F(X,alpha_role), F(alpha_role,Y))> let map_double_shares: MapRoleDoublePoly = session @@ -338,12 +329,7 @@ impl WrongDegreeSharingVss { let embedded_role = Z::embed_role_to_exceptional_sequence(r)?; let mut vec_map = Vec::with_capacity(bivariate_poly.len()); for p in &bivariate_poly { - let share_in_x = p.partial_y_evaluation(embedded_role)?; - let share_in_y = p.partial_x_evaluation(embedded_role)?; - vec_map.push(DoublePoly { - share_in_x, - share_in_y, - }); + vec_map.push(DoublePoly::from_bivariate(p, embedded_role)); } Ok::<(Role, Vec>), anyhow::Error>((*r, vec_map)) }) diff --git a/core/threshold-execution/src/small_execution/prss.rs b/core/threshold-execution/src/small_execution/prss.rs index 2817c7c891..2ae2dbd10b 100644 --- a/core/threshold-execution/src/small_execution/prss.rs +++ b/core/threshold-execution/src/small_execution/prss.rs @@ -21,7 +21,7 @@ use crate::{ }, }; use algebra::{ - bivariate::{MatrixMul, compute_powers_list}, + matrix::{MatrixMul, compute_powers_list}, poly::Poly, structure_traits::{ErrorCorrect, Invert, Ring, RingWithExceptionalSequence}, }; diff --git a/dylint.toml b/dylint.toml index 0297b0ecd4..81bbf2da9c 100644 --- a/dylint.toml +++ b/dylint.toml @@ -5,4 +5,3 @@ libraries = [ [non_local_effect_before_error_return] work_limit = 100000000 -