From 3b5d83cda8dc8aecad5b0e3de59f5a40d1270f53 Mon Sep 17 00:00:00 2001 From: David Palm Date: Wed, 6 May 2026 15:34:09 +0200 Subject: [PATCH 01/20] add: benches for bivariate ops and matmul --- Cargo.lock | 1 + core/threshold-algebra/Cargo.toml | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/Cargo.lock b/Cargo.lock index 12aa83646d..f3f16119aa 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -7579,6 +7579,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 a0f9ca8dbb..da0d54a038 100644 --- a/core/threshold-algebra/Cargo.toml +++ b/core/threshold-algebra/Cargo.toml @@ -24,11 +24,16 @@ 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 +[[bench]] +name = "bivariate" +harness = false + [features] default = ["extension_degree_4"] extension_degree_3 =[] From c3e724deee2a162adc1f57c1ae03d74794dcae48 Mon Sep 17 00:00:00 2001 From: David Palm Date: Wed, 6 May 2026 15:34:21 +0200 Subject: [PATCH 02/20] add: benches for bivariate ops and matmul --- core/threshold-algebra/benches/bivariate.rs | 133 ++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 core/threshold-algebra/benches/bivariate.rs diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs new file mode 100644 index 0000000000..434949bae0 --- /dev/null +++ b/core/threshold-algebra/benches/bivariate.rs @@ -0,0 +1,133 @@ +use aes_prng::AesRng; +use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; +use ndarray::{ArrayD, IxDyn}; +use rand::SeedableRng; +use std::hint::black_box; +use threshold_algebra::{ + bivariate::{BivariateEval, BivariatePoly, MatrixMul, compute_powers}, + galois_rings::degree_4::ResiduePolyF4Z128, + structure_traits::{Ring, Sample}, +}; + +const DEGREES: [usize; 5] = [1, 3, 4, 13, 20]; + +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).unwrap(); + (poly, point) +} + +fn matrix_setup(degree: usize) -> (ArrayD, ArrayD) { + let mut rng = AesRng::seed_from_u64(100 + degree as u64); + let d = degree + 1; + let vector = (0..d) + .map(|_| ResiduePolyF4Z128::sample(&mut rng)) + .collect(); + let matrix = (0..d * d) + .map(|_| ResiduePolyF4Z128::sample(&mut rng)) + .collect(); + + ( + ArrayD::from_shape_vec(IxDyn(&[d]), vector).unwrap(), + ArrayD::from_shape_vec(IxDyn(&[d, d]), matrix).unwrap(), + ) +} + +fn full_evaluation_direct( + poly: &BivariatePoly, + degree: usize, + alpha_x: Z, + alpha_y: Z, +) -> Z { + let powers_x = compute_powers(alpha_x, degree); + let powers_y = compute_powers(alpha_y, degree); + let mut acc = Z::ZERO; + + for (row_idx, row) in poly.coefs.rows().into_iter().enumerate() { + for (col_idx, coef) in row.iter().enumerate() { + acc += powers_x[row_idx] * *coef * powers_y[col_idx]; + } + } + + acc +} + +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)) + .unwrap(), + ) + }); + }); + } + + 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); + let current = poly.full_evaluation(point, point).unwrap(); + let direct = full_evaluation_direct(&poly, degree, point, point); + assert_eq!(current, direct); + + group.bench_function(BenchmarkId::new("partial_x", degree), |b| { + b.iter(|| black_box(poly.partial_x_evaluation(black_box(point)).unwrap())); + }); + group.bench_function(BenchmarkId::new("partial_y", degree), |b| { + b.iter(|| black_box(poly.partial_y_evaluation(black_box(point)).unwrap())); + }); + group.bench_function(BenchmarkId::new("full", degree), |b| { + b.iter(|| { + black_box( + poly.full_evaluation(black_box(point), black_box(point)) + .unwrap(), + ) + }); + }); + group.bench_function(BenchmarkId::new("full_direct", degree), |b| { + b.iter(|| black_box(full_evaluation_direct(&poly, degree, point, point))); + }); + } + + group.finish(); +} + +fn bench_matrix_mul(c: &mut Criterion) { + let mut group = c.benchmark_group("bivariate/matrix_mul"); + + for degree in DEGREES { + let (vector, matrix) = matrix_setup(degree); + + group.bench_function(BenchmarkId::new("vector_matrix", degree), |b| { + b.iter(|| black_box(vector.matmul(black_box(&matrix)).unwrap())); + }); + group.bench_function(BenchmarkId::new("matrix_vector", degree), |b| { + b.iter(|| black_box(matrix.matmul(black_box(&vector)).unwrap())); + }); + group.bench_function(BenchmarkId::new("vector_dot", degree), |b| { + b.iter(|| black_box(vector.matmul(black_box(&vector)).unwrap())); + }); + } + + group.finish(); +} + +criterion_group!( + bivariate, + bench_bivariate_sampling, + bench_bivariate_evaluation, + bench_matrix_mul +); +criterion_main!(bivariate); From 9e799164c51acf1553d98e423ed529f099dd45a1 Mon Sep 17 00:00:00 2001 From: David Palm Date: Wed, 6 May 2026 17:04:12 +0200 Subject: [PATCH 03/20] wip: optimizations sweep 1 --- core/threshold-algebra/benches/bivariate.rs | 13 +- core/threshold-algebra/src/bivariate.rs | 323 ++++++++++++-------- 2 files changed, 209 insertions(+), 127 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 434949bae0..71e0291af9 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -1,6 +1,6 @@ use aes_prng::AesRng; use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; -use ndarray::{ArrayD, IxDyn}; +use ndarray::{Array1, Array2}; use rand::SeedableRng; use std::hint::black_box; use threshold_algebra::{ @@ -10,6 +10,8 @@ use threshold_algebra::{ }; const DEGREES: [usize; 5] = [1, 3, 4, 13, 20]; +// Use this shorter sweep for fast local iteration. +// const DEGREES: [usize; 2] = [4, 20]; fn bivariate_setup(degree: usize) -> (BivariatePoly, ResiduePolyF4Z128) { let mut rng = AesRng::seed_from_u64(degree as u64); @@ -19,7 +21,7 @@ fn bivariate_setup(degree: usize) -> (BivariatePoly, ResidueP (poly, point) } -fn matrix_setup(degree: usize) -> (ArrayD, ArrayD) { +fn matrix_setup(degree: usize) -> (Array1, Array2) { let mut rng = AesRng::seed_from_u64(100 + degree as u64); let d = degree + 1; let vector = (0..d) @@ -30,8 +32,8 @@ fn matrix_setup(degree: usize) -> (ArrayD, ArrayD { - pub coefs: ArrayD, + pub coefs: Array2, degree: usize, } @@ -29,11 +26,13 @@ impl BivariatePoly { Z: Sample + Zero + Copy, { let d = degree + 1; - let coefs: Vec<_> = (0..d * d) - .map(|i| if i == 0 { secret } else { Z::sample(rng) }) - .collect(); + let n = d * d; + let mut coefs = Vec::with_capacity(n); + coefs.push(secret); + coefs.extend((1..n).map(|_| Z::sample(rng))); + Ok(BivariatePoly { - coefs: ArrayD::from_shape_vec(IxDyn(&[d, d]), coefs)?.into_dyn(), + coefs: Array2::from_shape_vec((d, d), coefs)?, degree, }) } @@ -61,93 +60,178 @@ pub fn compute_powers + Copy>(point: Z, degree: usize) powers_of_point } -pub trait MatrixMul { - fn matmul(&self, rhs: &ArrayD) -> Result>; +pub trait MatrixMul: Sized { + type Output; + fn matmul(&self, rhs: &Rhs) -> Result; } -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()) - } - (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); +macro_rules! dot_col { + ($lhs:expr, $rhs:expr, $col:expr, $first:expr $(, $idx:expr)+) => {{ + let mut acc = $lhs[$first] * $rhs[($first, $col)]; + $(acc += $lhs[$idx] * $rhs[($idx, $col)];)+ + acc + }}; +} + +macro_rules! dot_row { + ($lhs:expr, $rhs:expr, $row:expr, $first:expr $(, $idx:expr)+) => {{ + let mut acc = $lhs[($row, $first)] * $rhs[$first]; + $(acc += $lhs[($row, $idx)] * $rhs[$idx];)+ + acc + }}; +} + +macro_rules! dot_vec { + ($lhs:expr, $rhs:expr, $first:expr $(, $idx:expr)+) => {{ + let mut acc = $lhs[$first] * $rhs[$first]; + $(acc += $lhs[$idx] * $rhs[$idx];)+ + acc + }}; +} + +impl MatrixMul> for Array1 { + type Output = Array1; + fn matmul(&self, rhs: &Array2) -> Result { + Ok(match rhs.dim() { + (2, 2) => Array1::from_vec(vec![ + dot_col!(self, rhs, 0, 0, 1), + dot_col!(self, rhs, 1, 0, 1), + ]), + (5, 5) => Array1::from_vec(vec![ + dot_col!(self, rhs, 0, 0, 1, 2, 3, 4), + dot_col!(self, rhs, 1, 0, 1, 2, 3, 4), + dot_col!(self, rhs, 2, 0, 1, 2, 3, 4), + dot_col!(self, rhs, 3, 0, 1, 2, 3, 4), + dot_col!(self, rhs, 4, 0, 1, 2, 3, 4), + ]), + (_, cols) => { + let mut res = Vec::with_capacity(cols); + for col_idx in 0..cols { + let mut acc = Z::ZERO; + for row_idx in 0..self.len() { + acc += self[row_idx] * rhs[(row_idx, col_idx)]; } - Ok(Array::from_vec(res).into_dyn()) + res.push(acc); } + Array1::from_vec(res) } - (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); + }) + } +} + +impl MatrixMul> for Array2 { + type Output = Array1; + fn matmul(&self, rhs: &Array1) -> Result { + Ok(match self.dim() { + (2, 2) => Array1::from_vec(vec![ + dot_row!(self, rhs, 0, 0, 1), + dot_row!(self, rhs, 1, 0, 1), + ]), + (5, 5) => Array1::from_vec(vec![ + dot_row!(self, rhs, 0, 0, 1, 2, 3, 4), + dot_row!(self, rhs, 1, 0, 1, 2, 3, 4), + dot_row!(self, rhs, 2, 0, 1, 2, 3, 4), + dot_row!(self, rhs, 3, 0, 1, 2, 3, 4), + dot_row!(self, rhs, 4, 0, 1, 2, 3, 4), + ]), + (rows, cols) => { + let mut res = Vec::with_capacity(rows); + for row_idx in 0..rows { + let mut acc = Z::ZERO; + for col_idx in 0..cols { + acc += self[(row_idx, col_idx)] * rhs[col_idx]; } - Ok(Array::from_vec(res).into_dyn()) + res.push(acc); } + Array1::from_vec(res) } - (l_rank, r_rank) => Err(anyhow_error_and_log(format!( - "Matmul not implemented for tensors of rank {l_rank:?}, {r_rank:?}", - ))), - } + }) } } +// 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()) +// } +// (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 @@ -167,28 +251,37 @@ pub trait BivariateEval { } impl BivariateEval for BivariatePoly -where - ArrayD: MatrixMul, +// where +// ArrayD: MatrixMul, { fn partial_x_evaluation(&self, alpha: Z) -> Result> { - let powers_array = Array::from(compute_powers(alpha, self.degree)).into_dyn(); + let powers_array = Array1::from(compute_powers(alpha, self.degree)); let res_vector = powers_array.matmul(&self.coefs)?; Ok(Poly::from_coefs(res_vector.into_raw_vec_and_offset().0)) } fn partial_y_evaluation(&self, alpha: Z) -> Result> { - let powers_array = Array::from(compute_powers(alpha, self.degree)).into_dyn(); + let powers_array = Array1::from(compute_powers(alpha, self.degree)); let res_vector = self.coefs.matmul(&powers_array)?; Ok(Poly::from_coefs(res_vector.into_raw_vec_and_offset().0)) } 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 powers_array_x = Array1::from(compute_powers(alpha_x, self.degree)); + let powers_array_y = Array1::from(compute_powers(alpha_y, self.degree)); let lhs = powers_array_x.matmul(&self.coefs)?; - let res = lhs.matmul(&powers_array_y)?; - Ok(res[0]) + Ok(match lhs.len() { + 2 => dot_vec!(lhs, powers_array_y, 0, 1), + 5 => dot_vec!(lhs, powers_array_y, 0, 1, 2, 3, 4), + len => { + let mut acc = Z::ZERO; + for idx in 0..len { + acc += lhs[idx] * powers_array_y[idx]; + } + acc + } + }) } } @@ -211,30 +304,22 @@ mod tests { #[cfg(feature = "extension_degree_8")] use std::num::Wrapping; - //Checks that we error on incompatible sizes and dimensions + //Checks the hot matrix/vector 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); - - // test 1x1 vector mul errors - assert!(y4.matmul(&y2).is_err()); - assert!(y2.matmul(&y4).is_err()); + fn test_matmul_supported_shapes() { + let y2 = Array1::from_elem(2, ResiduePolyF4Z128::ONE); + let z22 = Array2::from_elem((2, 2), ResiduePolyF4Z128::ONE); + let two = ResiduePolyF4Z128::ONE + ResiduePolyF4Z128::ONE; + let expected2 = Array1::from_elem(2, two); + assert_eq!(y2.matmul(&z22).unwrap(), expected2); + assert_eq!(z22.matmul(&y2).unwrap(), expected2); + + let y5 = Array1::from_elem(5, ResiduePolyF4Z128::ONE); + let z55 = Array2::from_elem((5, 5), ResiduePolyF4Z128::ONE); + let five = two + two + ResiduePolyF4Z128::ONE; + let expected5 = Array1::from_elem(5, five); + assert_eq!(y5.matmul(&z55).unwrap(), expected5); + assert_eq!(z55.matmul(&y5).unwrap(), expected5); } //Test that eval at 0 return the secret for ResiduePolyF4Z128 @@ -606,9 +691,7 @@ mod tests { ]; let bpoly = BivariatePoly { - coefs: ArrayD::from_shape_vec(IxDyn(&[5, 5]), coefs) - .unwrap() - .to_owned(), + coefs: Array2::from_shape_vec((5, 5), coefs).unwrap(), degree: 4, }; From dbc2767ae16dc8bb30f1e332279d14ed944da903 Mon Sep 17 00:00:00 2001 From: David Palm Date: Wed, 6 May 2026 18:08:58 +0200 Subject: [PATCH 04/20] wip: aggressive unrolling and removal of ndarray (almost) completely. --- core/threshold-algebra/benches/bivariate.rs | 100 ++++- core/threshold-algebra/src/bivariate.rs | 347 +++++++----------- core/threshold-algebra/src/lib.rs | 1 + core/threshold-algebra/src/matrix.rs | 96 +++++ .../src/large_execution/double_sharing.rs | 2 +- .../src/large_execution/single_sharing.rs | 3 +- .../src/small_execution/prss.rs | 3 +- 7 files changed, 331 insertions(+), 221 deletions(-) create mode 100644 core/threshold-algebra/src/matrix.rs diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 71e0291af9..6d5af11a2b 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -1,10 +1,9 @@ use aes_prng::AesRng; use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; -use ndarray::{Array1, Array2}; use rand::SeedableRng; use std::hint::black_box; use threshold_algebra::{ - bivariate::{BivariateEval, BivariatePoly, MatrixMul, compute_powers}, + bivariate::{BivariateEval, BivariatePoly, compute_powers}, galois_rings::degree_4::ResiduePolyF4Z128, structure_traits::{Ring, Sample}, }; @@ -21,7 +20,7 @@ fn bivariate_setup(degree: usize) -> (BivariatePoly, ResidueP (poly, point) } -fn matrix_setup(degree: usize) -> (Array1, Array2) { +fn matrix_setup(degree: usize) -> (Vec, Vec) { let mut rng = AesRng::seed_from_u64(100 + degree as u64); let d = degree + 1; let vector = (0..d) @@ -31,10 +30,7 @@ fn matrix_setup(degree: usize) -> (Array1, Array2( @@ -46,8 +42,9 @@ fn full_evaluation_direct( let powers_x = compute_powers(alpha_x, degree); let powers_y = compute_powers(alpha_y, degree); let mut acc = Z::ZERO; + let d = degree + 1; - for (row_idx, row) in poly.coefs.rows().into_iter().enumerate() { + for (row_idx, row) in poly.coefs.chunks_exact(d).enumerate() { for (col_idx, coef) in row.iter().enumerate() { acc += powers_x[row_idx] * *coef * powers_y[col_idx]; } @@ -56,6 +53,79 @@ fn full_evaluation_direct( acc } +macro_rules! dot_vector_matrix_col { + ($vector:expr, $matrix:expr, $d:expr, $col:expr, $first:expr $(, $row:expr)+) => {{ + let mut acc = $vector[$first] * $matrix[$first * $d + $col]; + $(acc += $vector[$row] * $matrix[$row * $d + $col];)+ + acc + }}; +} + +macro_rules! dot_matrix_vector_row { + ($matrix:expr, $vector:expr, $d:expr, $row:expr, $first:expr $(, $col:expr)+) => {{ + let row_start = $row * $d; + let mut acc = $matrix[row_start + $first] * $vector[$first]; + $(acc += $matrix[row_start + $col] * $vector[$col];)+ + acc + }}; +} + +fn vector_matrix_product(vector: &[Z], matrix: &[Z]) -> Vec { + let d = vector.len(); + match d { + 2 => vec![ + dot_vector_matrix_col!(vector, matrix, d, 0, 0, 1), + dot_vector_matrix_col!(vector, matrix, d, 1, 0, 1), + ], + 5 => vec![ + dot_vector_matrix_col!(vector, matrix, d, 0, 0, 1, 2, 3, 4), + dot_vector_matrix_col!(vector, matrix, d, 1, 0, 1, 2, 3, 4), + dot_vector_matrix_col!(vector, matrix, d, 2, 0, 1, 2, 3, 4), + dot_vector_matrix_col!(vector, matrix, d, 3, 0, 1, 2, 3, 4), + dot_vector_matrix_col!(vector, matrix, d, 4, 0, 1, 2, 3, 4), + ], + _ => { + let mut res = vec![Z::ZERO; d]; + for (row_idx, vector_value) in vector.iter().enumerate() { + let row_start = row_idx * d; + for col_idx in 0..d { + res[col_idx] += *vector_value * matrix[row_start + col_idx]; + } + } + res + } + } +} + +fn matrix_vector_product(matrix: &[Z], vector: &[Z]) -> Vec { + let d = vector.len(); + match d { + 2 => vec![ + dot_matrix_vector_row!(matrix, vector, d, 0, 0, 1), + dot_matrix_vector_row!(matrix, vector, d, 1, 0, 1), + ], + 5 => vec![ + dot_matrix_vector_row!(matrix, vector, d, 0, 0, 1, 2, 3, 4), + dot_matrix_vector_row!(matrix, vector, d, 1, 0, 1, 2, 3, 4), + dot_matrix_vector_row!(matrix, vector, d, 2, 0, 1, 2, 3, 4), + dot_matrix_vector_row!(matrix, vector, d, 3, 0, 1, 2, 3, 4), + dot_matrix_vector_row!(matrix, vector, d, 4, 0, 1, 2, 3, 4), + ], + _ => { + let mut res = Vec::with_capacity(d); + for row_idx in 0..d { + let row_start = row_idx * d; + let mut acc = Z::ZERO; + for col_idx in 0..d { + acc += matrix[row_start + col_idx] * vector[col_idx]; + } + res.push(acc); + } + res + } + } +} + fn bench_bivariate_sampling(c: &mut Criterion) { let mut group = c.benchmark_group("bivariate/from_secret"); @@ -113,10 +183,20 @@ fn bench_matrix_mul(c: &mut Criterion) { let (vector, matrix) = matrix_setup(degree); group.bench_function(BenchmarkId::new("vector_matrix", degree), |b| { - b.iter(|| black_box(vector.matmul(black_box(&matrix)).unwrap())); + b.iter(|| { + black_box(vector_matrix_product( + black_box(&vector), + black_box(&matrix), + )) + }); }); group.bench_function(BenchmarkId::new("matrix_vector", degree), |b| { - b.iter(|| black_box(matrix.matmul(black_box(&vector)).unwrap())); + b.iter(|| { + black_box(matrix_vector_product( + black_box(&matrix), + black_box(&vector), + )) + }); }); } diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 3dced71115..2cc96914c4 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -5,17 +5,15 @@ use super::structure_traits::Sample; use super::structure_traits::Zero; use anyhow::Result; -use ndarray::Array1; -use ndarray::Array2; use rand::{CryptoRng, Rng}; use std::ops::Mul; -/// Bivariate polynomial is a matrix of coefficients of ResiduePolynomials +/// Bivariate polynomial is a row-major 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)] pub struct BivariatePoly { - pub coefs: Array2, + pub coefs: Vec, degree: usize, } @@ -31,10 +29,7 @@ impl BivariatePoly { coefs.push(secret); coefs.extend((1..n).map(|_| Z::sample(rng))); - Ok(BivariatePoly { - coefs: Array2::from_shape_vec((d, d), coefs)?, - degree, - }) + Ok(BivariatePoly { coefs, degree }) } } @@ -60,188 +55,36 @@ pub fn compute_powers + Copy>(point: Z, degree: usize) powers_of_point } -pub trait MatrixMul: Sized { - type Output; - fn matmul(&self, rhs: &Rhs) -> Result; -} - -macro_rules! dot_col { - ($lhs:expr, $rhs:expr, $col:expr, $first:expr $(, $idx:expr)+) => {{ - let mut acc = $lhs[$first] * $rhs[($first, $col)]; - $(acc += $lhs[$idx] * $rhs[($idx, $col)];)+ +macro_rules! eval_row_y { + ($coefs:expr, $start:expr, $alpha:expr, $last:expr $(, $idx:expr)+) => {{ + let mut acc = $coefs[$start + $last]; + $(acc = acc * $alpha + $coefs[$start + $idx];)+ acc }}; } -macro_rules! dot_row { - ($lhs:expr, $rhs:expr, $row:expr, $first:expr $(, $idx:expr)+) => {{ - let mut acc = $lhs[($row, $first)] * $rhs[$first]; - $(acc += $lhs[($row, $idx)] * $rhs[$idx];)+ - acc +macro_rules! eval_col_x { + ($coefs:expr, $col:expr, $d:expr, $alpha:expr, $last:expr $(, $idx:expr)+) => {{ + let mut acc = $alpha * $coefs[$last * $d + $col]; + $(acc = $alpha * ($coefs[$idx * $d + $col] + acc);)+ + $coefs[$col] + acc }}; } -macro_rules! dot_vec { - ($lhs:expr, $rhs:expr, $first:expr $(, $idx:expr)+) => {{ - let mut acc = $lhs[$first] * $rhs[$first]; - $(acc += $lhs[$idx] * $rhs[$idx];)+ - acc - }}; -} - -impl MatrixMul> for Array1 { - type Output = Array1; - fn matmul(&self, rhs: &Array2) -> Result { - Ok(match rhs.dim() { - (2, 2) => Array1::from_vec(vec![ - dot_col!(self, rhs, 0, 0, 1), - dot_col!(self, rhs, 1, 0, 1), - ]), - (5, 5) => Array1::from_vec(vec![ - dot_col!(self, rhs, 0, 0, 1, 2, 3, 4), - dot_col!(self, rhs, 1, 0, 1, 2, 3, 4), - dot_col!(self, rhs, 2, 0, 1, 2, 3, 4), - dot_col!(self, rhs, 3, 0, 1, 2, 3, 4), - dot_col!(self, rhs, 4, 0, 1, 2, 3, 4), - ]), - (_, cols) => { - let mut res = Vec::with_capacity(cols); - for col_idx in 0..cols { - let mut acc = Z::ZERO; - for row_idx in 0..self.len() { - acc += self[row_idx] * rhs[(row_idx, col_idx)]; - } - res.push(acc); - } - Array1::from_vec(res) - } - }) - } -} - -impl MatrixMul> for Array2 { - type Output = Array1; - fn matmul(&self, rhs: &Array1) -> Result { - Ok(match self.dim() { - (2, 2) => Array1::from_vec(vec![ - dot_row!(self, rhs, 0, 0, 1), - dot_row!(self, rhs, 1, 0, 1), - ]), - (5, 5) => Array1::from_vec(vec![ - dot_row!(self, rhs, 0, 0, 1, 2, 3, 4), - dot_row!(self, rhs, 1, 0, 1, 2, 3, 4), - dot_row!(self, rhs, 2, 0, 1, 2, 3, 4), - dot_row!(self, rhs, 3, 0, 1, 2, 3, 4), - dot_row!(self, rhs, 4, 0, 1, 2, 3, 4), - ]), - (rows, cols) => { - let mut res = Vec::with_capacity(rows); - for row_idx in 0..rows { - let mut acc = Z::ZERO; - for col_idx in 0..cols { - acc += self[(row_idx, col_idx)] * rhs[col_idx]; - } - res.push(acc); - } - Array1::from_vec(res) - } - }) +impl BivariatePoly { + fn dimension(&self) -> usize { + self.degree + 1 } } -// 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()) -// } -// (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>; @@ -250,34 +93,97 @@ pub trait BivariateEval { fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result; } -impl BivariateEval for BivariatePoly -// where -// ArrayD: MatrixMul, -{ +impl BivariateEval for BivariatePoly { fn partial_x_evaluation(&self, alpha: Z) -> Result> { - let powers_array = Array1::from(compute_powers(alpha, self.degree)); - let res_vector = powers_array.matmul(&self.coefs)?; - Ok(Poly::from_coefs(res_vector.into_raw_vec_and_offset().0)) + let d = self.dimension(); + let coefs = self.coefs.as_slice(); + let res = match d { + 2 => vec![coefs[0] + alpha * coefs[2], coefs[1] + alpha * coefs[3]], + 5 => vec![ + eval_col_x!(coefs, 0, d, alpha, 4, 3, 2, 1), + eval_col_x!(coefs, 1, d, alpha, 4, 3, 2, 1), + eval_col_x!(coefs, 2, d, alpha, 4, 3, 2, 1), + eval_col_x!(coefs, 3, d, alpha, 4, 3, 2, 1), + eval_col_x!(coefs, 4, d, alpha, 4, 3, 2, 1), + ], + _ => { + let mut res = vec![Z::ZERO; d]; + let mut alpha_power = Z::ONE; + for row_idx in 0..d { + let row_start = row_idx * d; + for col_idx in 0..d { + res[col_idx] += alpha_power * coefs[row_start + col_idx]; + } + alpha_power *= alpha; + } + res + } + }; + Ok(Poly::from_coefs(res)) } fn partial_y_evaluation(&self, alpha: Z) -> Result> { - let powers_array = Array1::from(compute_powers(alpha, self.degree)); - let res_vector = self.coefs.matmul(&powers_array)?; - Ok(Poly::from_coefs(res_vector.into_raw_vec_and_offset().0)) + let d = self.dimension(); + let coefs = self.coefs.as_slice(); + let res = match d { + 2 => vec![ + eval_row_y!(coefs, 0, alpha, 1, 0), + eval_row_y!(coefs, 2, alpha, 1, 0), + ], + 5 => vec![ + eval_row_y!(coefs, 0, alpha, 4, 3, 2, 1, 0), + eval_row_y!(coefs, 5, alpha, 4, 3, 2, 1, 0), + eval_row_y!(coefs, 10, alpha, 4, 3, 2, 1, 0), + eval_row_y!(coefs, 15, alpha, 4, 3, 2, 1, 0), + eval_row_y!(coefs, 20, alpha, 4, 3, 2, 1, 0), + ], + _ => { + let mut res = Vec::with_capacity(d); + for row_idx in 0..d { + let row_start = row_idx * d; + let mut acc = coefs[row_start + d - 1]; + for col_idx in (0..d - 1).rev() { + acc = acc * alpha + coefs[row_start + col_idx]; + } + res.push(acc); + } + res + } + }; + Ok(Poly::from_coefs(res)) } fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result { - let powers_array_x = Array1::from(compute_powers(alpha_x, self.degree)); - let powers_array_y = Array1::from(compute_powers(alpha_y, self.degree)); - - let lhs = powers_array_x.matmul(&self.coefs)?; - Ok(match lhs.len() { - 2 => dot_vec!(lhs, powers_array_y, 0, 1), - 5 => dot_vec!(lhs, powers_array_y, 0, 1, 2, 3, 4), - len => { + let d = self.dimension(); + let coefs = self.coefs.as_slice(); + Ok(match d { + 2 => { + let row0 = eval_row_y!(coefs, 0, alpha_y, 1, 0); + let row1 = eval_row_y!(coefs, 2, alpha_y, 1, 0); + row0 + alpha_x * row1 + } + 5 => { + let row0 = eval_row_y!(coefs, 0, alpha_y, 4, 3, 2, 1, 0); + let row1 = eval_row_y!(coefs, 5, alpha_y, 4, 3, 2, 1, 0); + let row2 = eval_row_y!(coefs, 10, alpha_y, 4, 3, 2, 1, 0); + let row3 = eval_row_y!(coefs, 15, alpha_y, 4, 3, 2, 1, 0); + let row4 = eval_row_y!(coefs, 20, alpha_y, 4, 3, 2, 1, 0); + let alpha_x2 = alpha_x * alpha_x; + let alpha_x3 = alpha_x2 * alpha_x; + let alpha_x4 = alpha_x3 * alpha_x; + row0 + alpha_x * row1 + alpha_x2 * row2 + alpha_x3 * row3 + alpha_x4 * row4 + } + _ => { let mut acc = Z::ZERO; - for idx in 0..len { - acc += lhs[idx] * powers_array_y[idx]; + let mut alpha_x_power = Z::ONE; + for row_idx in 0..d { + let row_start = row_idx * d; + let mut row_acc = coefs[row_start + d - 1]; + for col_idx in (0..d - 1).rev() { + row_acc = row_acc * alpha_y + coefs[row_start + col_idx]; + } + acc += alpha_x_power * row_acc; + alpha_x_power *= alpha_x; } acc } @@ -304,22 +210,50 @@ mod tests { #[cfg(feature = "extension_degree_8")] use std::num::Wrapping; - //Checks the hot matrix/vector shapes used by bivariate evaluation. + //Checks the hot coefficient shapes used by bivariate evaluation. #[test] - fn test_matmul_supported_shapes() { - let y2 = Array1::from_elem(2, ResiduePolyF4Z128::ONE); - let z22 = Array2::from_elem((2, 2), ResiduePolyF4Z128::ONE); + fn test_bivariate_supported_shapes() { let two = ResiduePolyF4Z128::ONE + ResiduePolyF4Z128::ONE; - let expected2 = Array1::from_elem(2, two); - assert_eq!(y2.matmul(&z22).unwrap(), expected2); - assert_eq!(z22.matmul(&y2).unwrap(), expected2); + let bpoly2 = BivariatePoly { + coefs: vec![ResiduePolyF4Z128::ONE; 4], + degree: 1, + }; + let expected2 = Poly::from_coefs(vec![two; 2]); + assert_eq!( + bpoly2.partial_x_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + expected2 + ); + assert_eq!( + bpoly2.partial_y_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + expected2 + ); + assert_eq!( + bpoly2 + .full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE) + .unwrap(), + two + two + ); - let y5 = Array1::from_elem(5, ResiduePolyF4Z128::ONE); - let z55 = Array2::from_elem((5, 5), ResiduePolyF4Z128::ONE); let five = two + two + ResiduePolyF4Z128::ONE; - let expected5 = Array1::from_elem(5, five); - assert_eq!(y5.matmul(&z55).unwrap(), expected5); - assert_eq!(z55.matmul(&y5).unwrap(), expected5); + let bpoly5 = BivariatePoly { + coefs: vec![ResiduePolyF4Z128::ONE; 25], + degree: 4, + }; + let expected5 = Poly::from_coefs(vec![five; 5]); + assert_eq!( + bpoly5.partial_x_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + expected5 + ); + assert_eq!( + bpoly5.partial_y_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + expected5 + ); + assert_eq!( + bpoly5 + .full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE) + .unwrap(), + five + five + five + five + five + ); } //Test that eval at 0 return the secret for ResiduePolyF4Z128 @@ -690,10 +624,7 @@ mod tests { }, //x4y4 ]; - let bpoly = BivariatePoly { - coefs: Array2::from_shape_vec((5, 5), coefs).unwrap(), - degree: 4, - }; + let bpoly = BivariatePoly { coefs, degree: 4 }; let point = ResiduePoly { coefs: [ 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..47809c2e36 --- /dev/null +++ b/core/threshold-algebra/src/matrix.rs @@ -0,0 +1,96 @@ +use crate::structure_traits::Ring; + +use anyhow::Result; +use error_utils::anyhow_error_and_log; +use itertools::Itertools; +use ndarray::{Array, ArrayD, IxDyn}; + +pub trait MatrixMul: Sized { + type Output; + fn matmul(&self, rhs: &Rhs) -> Result; +} + +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:?}", + ))), + } + } +} 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 b52182eb08..570f15b1dd 100644 --- a/core/threshold-execution/src/large_execution/single_sharing.rs +++ b/core/threshold-execution/src/large_execution/single_sharing.rs @@ -1,7 +1,8 @@ use super::local_single_share::{LocalSingleShare, SecureLocalSingleShare}; use crate::runtime::sessions::large_session::LargeSessionHandles; use algebra::{ - bivariate::{MatrixMul, compute_powers}, + bivariate::compute_powers, + matrix::MatrixMul, structure_traits::{Derive, ErrorCorrect, Invert, Ring, RingWithExceptionalSequence}, }; use async_trait::async_trait; diff --git a/core/threshold-execution/src/small_execution/prss.rs b/core/threshold-execution/src/small_execution/prss.rs index f2594edaae..2d6e578c36 100644 --- a/core/threshold-execution/src/small_execution/prss.rs +++ b/core/threshold-execution/src/small_execution/prss.rs @@ -21,7 +21,8 @@ use crate::{ }, }; use algebra::{ - bivariate::{MatrixMul, compute_powers_list}, + bivariate::compute_powers_list, + matrix::MatrixMul, poly::Poly, structure_traits::{ErrorCorrect, Invert, Ring, RingWithExceptionalSequence}, }; From cb554ce7cfe3b92ef5dbf8ed3b1e7ea7593b0499 Mon Sep 17 00:00:00 2001 From: David Palm Date: Thu, 7 May 2026 07:53:59 +0200 Subject: [PATCH 05/20] chore: remove aggressive unrolling and reformulate loops for speed. --- core/threshold-algebra/benches/bivariate.rs | 75 ++------ core/threshold-algebra/src/bivariate.rs | 189 +++++++++----------- 2 files changed, 101 insertions(+), 163 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 6d5af11a2b..6370428b2f 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -53,77 +53,30 @@ fn full_evaluation_direct( acc } -macro_rules! dot_vector_matrix_col { - ($vector:expr, $matrix:expr, $d:expr, $col:expr, $first:expr $(, $row:expr)+) => {{ - let mut acc = $vector[$first] * $matrix[$first * $d + $col]; - $(acc += $vector[$row] * $matrix[$row * $d + $col];)+ - acc - }}; -} - -macro_rules! dot_matrix_vector_row { - ($matrix:expr, $vector:expr, $d:expr, $row:expr, $first:expr $(, $col:expr)+) => {{ - let row_start = $row * $d; - let mut acc = $matrix[row_start + $first] * $vector[$first]; - $(acc += $matrix[row_start + $col] * $vector[$col];)+ - acc - }}; -} - fn vector_matrix_product(vector: &[Z], matrix: &[Z]) -> Vec { let d = vector.len(); - match d { - 2 => vec![ - dot_vector_matrix_col!(vector, matrix, d, 0, 0, 1), - dot_vector_matrix_col!(vector, matrix, d, 1, 0, 1), - ], - 5 => vec![ - dot_vector_matrix_col!(vector, matrix, d, 0, 0, 1, 2, 3, 4), - dot_vector_matrix_col!(vector, matrix, d, 1, 0, 1, 2, 3, 4), - dot_vector_matrix_col!(vector, matrix, d, 2, 0, 1, 2, 3, 4), - dot_vector_matrix_col!(vector, matrix, d, 3, 0, 1, 2, 3, 4), - dot_vector_matrix_col!(vector, matrix, d, 4, 0, 1, 2, 3, 4), - ], - _ => { - let mut res = vec![Z::ZERO; d]; - for (row_idx, vector_value) in vector.iter().enumerate() { - let row_start = row_idx * d; - for col_idx in 0..d { - res[col_idx] += *vector_value * matrix[row_start + col_idx]; - } - } - res + let mut res = vec![Z::ZERO; d]; + for (row_idx, vector_value) in vector.iter().enumerate() { + let row_start = row_idx * d; + for col_idx in 0..d { + res[col_idx] += *vector_value * matrix[row_start + col_idx]; } } + res } fn matrix_vector_product(matrix: &[Z], vector: &[Z]) -> Vec { let d = vector.len(); - match d { - 2 => vec![ - dot_matrix_vector_row!(matrix, vector, d, 0, 0, 1), - dot_matrix_vector_row!(matrix, vector, d, 1, 0, 1), - ], - 5 => vec![ - dot_matrix_vector_row!(matrix, vector, d, 0, 0, 1, 2, 3, 4), - dot_matrix_vector_row!(matrix, vector, d, 1, 0, 1, 2, 3, 4), - dot_matrix_vector_row!(matrix, vector, d, 2, 0, 1, 2, 3, 4), - dot_matrix_vector_row!(matrix, vector, d, 3, 0, 1, 2, 3, 4), - dot_matrix_vector_row!(matrix, vector, d, 4, 0, 1, 2, 3, 4), - ], - _ => { - let mut res = Vec::with_capacity(d); - for row_idx in 0..d { - let row_start = row_idx * d; - let mut acc = Z::ZERO; - for col_idx in 0..d { - acc += matrix[row_start + col_idx] * vector[col_idx]; - } - res.push(acc); - } - res + let mut res = Vec::with_capacity(d); + for row_idx in 0..d { + let row_start = row_idx * d; + let mut acc = Z::ZERO; + for col_idx in 0..d { + acc += matrix[row_start + col_idx] * vector[col_idx]; } + res.push(acc); } + res } fn bench_bivariate_sampling(c: &mut Criterion) { diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 2cc96914c4..3100295b00 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -38,7 +38,7 @@ pub fn compute_powers_list + Copy>( points: &[F], max_exponent: usize, ) -> Vec> { - let mut alpha_powers = Vec::new(); + let mut alpha_powers = Vec::with_capacity(points.len()); for p in points { alpha_powers.push(compute_powers(*p, max_exponent)); } @@ -47,7 +47,7 @@ pub fn compute_powers_list + Copy>( /// 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(); + let mut powers_of_point = Vec::with_capacity(degree + 1); powers_of_point.push(Z::ONE); // start with for i in 1..=degree { powers_of_point.push(powers_of_point[i - 1] * point); @@ -55,22 +55,6 @@ pub fn compute_powers + Copy>(point: Z, degree: usize) powers_of_point } -macro_rules! eval_row_y { - ($coefs:expr, $start:expr, $alpha:expr, $last:expr $(, $idx:expr)+) => {{ - let mut acc = $coefs[$start + $last]; - $(acc = acc * $alpha + $coefs[$start + $idx];)+ - acc - }}; -} - -macro_rules! eval_col_x { - ($coefs:expr, $col:expr, $d:expr, $alpha:expr, $last:expr $(, $idx:expr)+) => {{ - let mut acc = $alpha * $coefs[$last * $d + $col]; - $(acc = $alpha * ($coefs[$idx * $d + $col] + acc);)+ - $coefs[$col] + acc - }}; -} - impl BivariatePoly { fn dimension(&self) -> usize { self.degree + 1 @@ -78,16 +62,25 @@ impl BivariatePoly { } pub trait BivariateEval { - /// Given a degree T bivariate poly F(X,Y) and a point \alpha, we compute - /// G(X) = F(X, \alpha) as - /// [sum(alpha^j * a_{j0}), ..., sum(alpha^j * a_{jd})] + /// Given a degree T bivariate poly F(X,Y) and a point \alpha, evaluate the X variable: + /// G(Y) = F(\alpha, Y). 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 - /// [sum(alpha^j * a_{0j}), ..., sum(alpha^j * a_{dj})] + /// Given a degree T bivariate poly F(X,Y) and a point \alpha, evaluate the Y variable: + /// G(X) = F(X, \alpha). fn partial_y_evaluation(&self, alpha: Z) -> Result>; + /// Compute both partial evaluations at the same point. + /// + /// Returns `(F(\alpha, Y), F(X, \alpha))`, matching + /// [`Self::partial_x_evaluation`] and [`Self::partial_y_evaluation`]. + fn partial_evaluations(&self, alpha: Z) -> Result<(Poly, Poly)> { + Ok(( + self.partial_x_evaluation(alpha)?, + self.partial_y_evaluation(alpha)?, + )) + } + /// 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; @@ -96,98 +89,74 @@ pub trait BivariateEval { impl BivariateEval for BivariatePoly { fn partial_x_evaluation(&self, alpha: Z) -> Result> { let d = self.dimension(); - let coefs = self.coefs.as_slice(); - let res = match d { - 2 => vec![coefs[0] + alpha * coefs[2], coefs[1] + alpha * coefs[3]], - 5 => vec![ - eval_col_x!(coefs, 0, d, alpha, 4, 3, 2, 1), - eval_col_x!(coefs, 1, d, alpha, 4, 3, 2, 1), - eval_col_x!(coefs, 2, d, alpha, 4, 3, 2, 1), - eval_col_x!(coefs, 3, d, alpha, 4, 3, 2, 1), - eval_col_x!(coefs, 4, d, alpha, 4, 3, 2, 1), - ], - _ => { - let mut res = vec![Z::ZERO; d]; - let mut alpha_power = Z::ONE; - for row_idx in 0..d { - let row_start = row_idx * d; - for col_idx in 0..d { - res[col_idx] += alpha_power * coefs[row_start + col_idx]; - } - alpha_power *= alpha; - } - res + let coefs = &self.coefs[..d * d]; + let mut res = coefs[(d - 1) * d..d * d].to_vec(); + 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; } - }; + } Ok(Poly::from_coefs(res)) } fn partial_y_evaluation(&self, alpha: Z) -> Result> { let d = self.dimension(); - let coefs = self.coefs.as_slice(); - let res = match d { - 2 => vec![ - eval_row_y!(coefs, 0, alpha, 1, 0), - eval_row_y!(coefs, 2, alpha, 1, 0), - ], - 5 => vec![ - eval_row_y!(coefs, 0, alpha, 4, 3, 2, 1, 0), - eval_row_y!(coefs, 5, alpha, 4, 3, 2, 1, 0), - eval_row_y!(coefs, 10, alpha, 4, 3, 2, 1, 0), - eval_row_y!(coefs, 15, alpha, 4, 3, 2, 1, 0), - eval_row_y!(coefs, 20, alpha, 4, 3, 2, 1, 0), - ], - _ => { - let mut res = Vec::with_capacity(d); - for row_idx in 0..d { - let row_start = row_idx * d; - let mut acc = coefs[row_start + d - 1]; - for col_idx in (0..d - 1).rev() { - acc = acc * alpha + coefs[row_start + col_idx]; - } - res.push(acc); - } - res + let coefs = &self.coefs[..d * d]; + let mut res = Vec::with_capacity(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); + } Ok(Poly::from_coefs(res)) } - fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result { + fn partial_evaluations(&self, alpha: Z) -> Result<(Poly, Poly)> { let d = self.dimension(); - let coefs = self.coefs.as_slice(); - Ok(match d { - 2 => { - let row0 = eval_row_y!(coefs, 0, alpha_y, 1, 0); - let row1 = eval_row_y!(coefs, 2, alpha_y, 1, 0); - row0 + alpha_x * row1 + let coefs = &self.coefs[..d * d]; + let last_row_start = (d - 1) * d; + let last_row = &coefs[last_row_start..last_row_start + d]; + + let mut partial_x = last_row.to_vec(); + let mut partial_y = vec![Z::ZERO; d]; + + let mut acc = last_row[d - 1]; + for coef in last_row[..d - 1].iter().rev() { + acc = acc * alpha + *coef; + } + partial_y[d - 1] = acc; + + for (row_idx, row) in coefs[..last_row_start].chunks_exact(d).enumerate().rev() { + let mut acc = row[d - 1]; + for coef in row[..d - 1].iter().rev() { + acc = acc * alpha + *coef; } - 5 => { - let row0 = eval_row_y!(coefs, 0, alpha_y, 4, 3, 2, 1, 0); - let row1 = eval_row_y!(coefs, 5, alpha_y, 4, 3, 2, 1, 0); - let row2 = eval_row_y!(coefs, 10, alpha_y, 4, 3, 2, 1, 0); - let row3 = eval_row_y!(coefs, 15, alpha_y, 4, 3, 2, 1, 0); - let row4 = eval_row_y!(coefs, 20, alpha_y, 4, 3, 2, 1, 0); - let alpha_x2 = alpha_x * alpha_x; - let alpha_x3 = alpha_x2 * alpha_x; - let alpha_x4 = alpha_x3 * alpha_x; - row0 + alpha_x * row1 + alpha_x2 * row2 + alpha_x3 * row3 + alpha_x4 * row4 + partial_y[row_idx] = acc; + + for (res_coef, coef) in partial_x.iter_mut().zip(row) { + *res_coef = *res_coef * alpha + *coef; } - _ => { - let mut acc = Z::ZERO; - let mut alpha_x_power = Z::ONE; - for row_idx in 0..d { - let row_start = row_idx * d; - let mut row_acc = coefs[row_start + d - 1]; - for col_idx in (0..d - 1).rev() { - row_acc = row_acc * alpha_y + coefs[row_start + col_idx]; - } - acc += alpha_x_power * row_acc; - alpha_x_power *= alpha_x; - } - acc + } + + Ok((Poly::from_coefs(partial_x), Poly::from_coefs(partial_y))) + } + + fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result { + let d = self.dimension(); + let coefs = &self.coefs[..d * d]; + let mut acc = Z::ZERO; + let mut alpha_x_power = Z::ONE; + for row in coefs.chunks_exact(d) { + let mut row_acc = row[d - 1]; + for coef in row[..d - 1].iter().rev() { + row_acc = row_acc * alpha_y + *coef; } - }) + acc += alpha_x_power * row_acc; + alpha_x_power *= alpha_x; + } + Ok(acc) } } @@ -256,6 +225,22 @@ mod tests { ); } + #[rstest] + #[case(1)] + #[case(4)] + #[case(10)] + fn test_bivariate_partial_evaluations(#[case] degree: usize) { + let mut rng = AesRng::seed_from_u64(0); + let secret = ResiduePolyF4Z128::sample(&mut rng); + let point = ResiduePolyF4Z128::sample(&mut rng); + let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree).unwrap(); + + let (partial_x, partial_y) = bpoly.partial_evaluations(point).unwrap(); + + assert_eq!(partial_x, bpoly.partial_x_evaluation(point).unwrap()); + assert_eq!(partial_y, bpoly.partial_y_evaluation(point).unwrap()); + } + //Test that eval at 0 return the secret for ResiduePolyF4Z128 #[rstest] #[case(4)] From 2382591816571b1cf9416125542764765adc8770 Mon Sep 17 00:00:00 2001 From: David Palm Date: Thu, 7 May 2026 08:56:24 +0200 Subject: [PATCH 06/20] chore: better docs, sprinkled a few quesitions/notices for reviewers, added some validation where strictly needed, stop returning Result when code is infallible --- core/threshold-algebra/src/bivariate.rs | 185 +++++++++++++----------- 1 file changed, 99 insertions(+), 86 deletions(-) diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 3100295b00..a33c114ba9 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -2,26 +2,26 @@ 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 rand::{CryptoRng, Rng}; use std::ops::Mul; -/// Bivariate polynomial is a row-major matrix of coefficients of ResiduePolynomials +/// Bivariate polynomial represented as a row-major square matrix of coefficients. +/// /// The row view of the polynomials is the following: /// [[a_{00}, a_{01}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}]] #[derive(Clone, Default, Debug)] pub struct BivariatePoly { + /// Row-major coefficients. Contains `(degree + 1)^2` elements. pub 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 n = d * d; @@ -29,7 +29,7 @@ impl BivariatePoly { coefs.push(secret); coefs.extend((1..n).map(|_| Z::sample(rng))); - Ok(BivariatePoly { coefs, degree }) + BivariatePoly { coefs, degree } } } @@ -59,50 +59,73 @@ impl BivariatePoly { fn dimension(&self) -> usize { self.degree + 1 } + + fn coefficients(&self) -> &[Z] { + let d = self.dimension(); + let n = d * d; + // @reviewers The old ndarray-based matmul checked ranks and dimensions on + // every evaluation. The only supported bivariate shape is the square coefficient + // matrix implied by `degree`, so we check that invariant once here and let the hot + // loops below use `d` directly. + assert_eq!( + self.coefs.len(), + n, + "BivariatePoly expected {n} coefficients for degree {}, got {}", + self.degree, + self.coefs.len() + ); + &self.coefs + } } pub trait BivariateEval { - /// Given a degree T bivariate poly F(X,Y) and a point \alpha, evaluate the X variable: - /// G(Y) = F(\alpha, Y). - fn partial_x_evaluation(&self, alpha: Z) -> Result>; + /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the X variable: + /// G(Y) = F(\alpha, Y) with coefficients [sum_i a_i0 \alpha^i, ..., sum_i a_id \alpha^i]. + fn partial_x_evaluation(&self, alpha: Z) -> Poly; - /// Given a degree T bivariate poly F(X,Y) and a point \alpha, evaluate the Y variable: - /// G(X) = F(X, \alpha). - fn partial_y_evaluation(&self, alpha: Z) -> Result>; + /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the Y variable: + /// G(X) = F(X, \alpha) with coefficients [sum_j a_0j \alpha^j, ..., sum_j a_dj \alpha^j]. + fn partial_y_evaluation(&self, alpha: Z) -> Poly; /// Compute both partial evaluations at the same point. /// /// Returns `(F(\alpha, Y), F(X, \alpha))`, matching /// [`Self::partial_x_evaluation`] and [`Self::partial_y_evaluation`]. - fn partial_evaluations(&self, alpha: Z) -> Result<(Poly, Poly)> { - Ok(( - self.partial_x_evaluation(alpha)?, - self.partial_y_evaluation(alpha)?, - )) + fn partial_evaluations(&self, alpha: Z) -> (Poly, Poly) { + ( + self.partial_x_evaluation(alpha), + self.partial_y_evaluation(alpha), + ) } /// 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; + fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Z; } impl BivariateEval for BivariatePoly { - fn partial_x_evaluation(&self, alpha: Z) -> Result> { + fn partial_x_evaluation(&self, alpha: Z) -> Poly { let d = self.dimension(); - let coefs = &self.coefs[..d * d]; + let coefs = self.coefficients(); let mut res = coefs[(d - 1) * d..d * d].to_vec(); + // @reviewers The old implementation checked each column length inside generic + // `matmul`. After the invariant check in `coefficients`, every row has exactly `d` + // coefficients, and `chunks_exact(d)` encodes the only shape this evaluator supports. 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; } } - Ok(Poly::from_coefs(res)) + Poly::from_coefs(res) } - fn partial_y_evaluation(&self, alpha: Z) -> Result> { + fn partial_y_evaluation(&self, alpha: Z) -> Poly { let d = self.dimension(); - let coefs = &self.coefs[..d * d]; + let coefs = self.coefficients(); let mut res = Vec::with_capacity(d); + // @reviewers The previous generic matrix path checked row/vector dimensions before + // multiplication. Here `coefficients` asserts the only valid bivariate shape once, + // and each Horner 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() { @@ -110,12 +133,12 @@ impl BivariateEval for BivariatePoly { } res.push(acc); } - Ok(Poly::from_coefs(res)) + Poly::from_coefs(res) } - fn partial_evaluations(&self, alpha: Z) -> Result<(Poly, Poly)> { + fn partial_evaluations(&self, alpha: Z) -> (Poly, Poly) { let d = self.dimension(); - let coefs = &self.coefs[..d * d]; + let coefs = self.coefficients(); let last_row_start = (d - 1) * d; let last_row = &coefs[last_row_start..last_row_start + d]; @@ -128,6 +151,9 @@ impl BivariateEval for BivariatePoly { } partial_y[d - 1] = acc; + // @reviewers This fused path relies on the same `d * d` coefficient invariant as + // the separate evaluators. The only difference is that it computes both directions + // in one pass to avoid repeating the row scan. for (row_idx, row) in coefs[..last_row_start].chunks_exact(d).enumerate().rev() { let mut acc = row[d - 1]; for coef in row[..d - 1].iter().rev() { @@ -140,12 +166,12 @@ impl BivariateEval for BivariatePoly { } } - Ok((Poly::from_coefs(partial_x), Poly::from_coefs(partial_y))) + (Poly::from_coefs(partial_x), Poly::from_coefs(partial_y)) } - fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Result { + fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Z { let d = self.dimension(); - let coefs = &self.coefs[..d * d]; + let coefs = self.coefficients(); let mut acc = Z::ZERO; let mut alpha_x_power = Z::ONE; for row in coefs.chunks_exact(d) { @@ -156,7 +182,10 @@ impl BivariateEval for BivariatePoly { acc += alpha_x_power * row_acc; alpha_x_power *= alpha_x; } - Ok(acc) + // @reviewers This replaces the old vector-matrix-then-dot-product path. Both power + // vectors had length `d` by construction, so the final rank-1 length check bought no + // extra validation beyond the coefficient invariant checked above. + acc } } @@ -170,7 +199,7 @@ mod tests { common::ResiduePoly, degree_4::{ResiduePolyF4Z64, ResiduePolyF4Z128}, }, - structure_traits::One, + structure_traits::{One, Zero}, }; use aes_prng::AesRng; @@ -189,17 +218,15 @@ mod tests { }; let expected2 = Poly::from_coefs(vec![two; 2]); assert_eq!( - bpoly2.partial_x_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + bpoly2.partial_x_evaluation(ResiduePolyF4Z128::ONE), expected2 ); assert_eq!( - bpoly2.partial_y_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + bpoly2.partial_y_evaluation(ResiduePolyF4Z128::ONE), expected2 ); assert_eq!( - bpoly2 - .full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE) - .unwrap(), + bpoly2.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), two + two ); @@ -210,95 +237,81 @@ mod tests { }; let expected5 = Poly::from_coefs(vec![five; 5]); assert_eq!( - bpoly5.partial_x_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + bpoly5.partial_x_evaluation(ResiduePolyF4Z128::ONE), expected5 ); assert_eq!( - bpoly5.partial_y_evaluation(ResiduePolyF4Z128::ONE).unwrap(), + bpoly5.partial_y_evaluation(ResiduePolyF4Z128::ONE), expected5 ); assert_eq!( - bpoly5 - .full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE) - .unwrap(), + bpoly5.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), five + five + five + five + five ); } + #[test] + #[should_panic(expected = "BivariatePoly expected 4 coefficients for degree 1")] + fn test_bivariate_invalid_coefficient_count() { + let bpoly = BivariatePoly { + coefs: vec![ResiduePolyF4Z128::ONE; 3], + degree: 1, + }; + + let _ = bpoly.partial_x_evaluation(ResiduePolyF4Z128::ONE); + } + #[rstest] - #[case(1)] - #[case(4)] - #[case(10)] - fn test_bivariate_partial_evaluations(#[case] degree: usize) { + fn test_bivariate_partial_evaluations(#[values(1, 4, 10)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z128::sample(&mut rng); let point = ResiduePolyF4Z128::sample(&mut rng); - let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree).unwrap(); + let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - let (partial_x, partial_y) = bpoly.partial_evaluations(point).unwrap(); + let (partial_x, partial_y) = bpoly.partial_evaluations(point); - assert_eq!(partial_x, bpoly.partial_x_evaluation(point).unwrap()); - assert_eq!(partial_y, bpoly.partial_y_evaluation(point).unwrap()); + assert_eq!(partial_x, bpoly.partial_x_evaluation(point)); + assert_eq!(partial_y, bpoly.partial_y_evaluation(point)); } //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_evaluation(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_evaluation(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 bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); + let ev_one = bpoly.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE); let sum_coefs = bpoly.coefs.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 bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); + let ev_one = bpoly.full_evaluation(ResiduePolyF4Z64::ONE, ResiduePolyF4Z64::ONE); let sum_coefs = bpoly.coefs.iter().fold(ResiduePoly::ZERO, |acc, x| acc + x); assert_eq!(ev_one, sum_coefs); } @@ -632,7 +645,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_evaluation(point); let expected_result = Poly::::from_coefs(vec![ ResiduePoly { @@ -706,7 +719,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_evaluation(point); let expected_result = Poly::::from_coefs(vec![ ResiduePoly { @@ -779,9 +792,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_evaluation(point_x, point_y); - let expected_res = bpoly.partial_x_evaluation(point_x).unwrap().eval(&point_y); + let expected_res = bpoly.partial_x_evaluation(point_x).eval(&point_y); assert_eq!(res, expected_res); } From 11871c9ce175cf2dd686a54da9b4de29b5533059 Mon Sep 17 00:00:00 2001 From: David Palm Date: Thu, 7 May 2026 08:56:56 +0200 Subject: [PATCH 07/20] chore: update bivariate benches to new api --- core/threshold-algebra/benches/bivariate.rs | 24 +++++++++------------ 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 6370428b2f..78c3e52a63 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -16,7 +16,7 @@ fn bivariate_setup(degree: usize) -> (BivariatePoly, ResidueP 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).unwrap(); + let poly = BivariatePoly::from_secret(&mut rng, secret, degree); (poly, point) } @@ -87,10 +87,11 @@ fn bench_bivariate_sampling(c: &mut Criterion) { 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)) - .unwrap(), - ) + black_box(BivariatePoly::from_secret( + &mut rng, + black_box(secret), + black_box(degree), + )) }); }); } @@ -103,23 +104,18 @@ fn bench_bivariate_evaluation(c: &mut Criterion) { for degree in DEGREES { let (poly, point) = bivariate_setup(degree); - let current = poly.full_evaluation(point, point).unwrap(); + let current = poly.full_evaluation(point, point); let direct = full_evaluation_direct(&poly, degree, point, point); assert_eq!(current, direct); group.bench_function(BenchmarkId::new("partial_x", degree), |b| { - b.iter(|| black_box(poly.partial_x_evaluation(black_box(point)).unwrap())); + b.iter(|| black_box(poly.partial_x_evaluation(black_box(point)))); }); group.bench_function(BenchmarkId::new("partial_y", degree), |b| { - b.iter(|| black_box(poly.partial_y_evaluation(black_box(point)).unwrap())); + b.iter(|| black_box(poly.partial_y_evaluation(black_box(point)))); }); group.bench_function(BenchmarkId::new("full", degree), |b| { - b.iter(|| { - black_box( - poly.full_evaluation(black_box(point), black_box(point)) - .unwrap(), - ) - }); + b.iter(|| black_box(poly.full_evaluation(black_box(point), black_box(point)))); }); group.bench_function(BenchmarkId::new("full_direct", degree), |b| { b.iter(|| black_box(full_evaluation_direct(&poly, degree, point, point))); From 3853a18be05b84a9f490434095037fef0571b4df Mon Sep 17 00:00:00 2001 From: David Palm Date: Thu, 7 May 2026 09:18:06 +0200 Subject: [PATCH 08/20] perf: avoid a few allocations --- .../src/large_execution/single_sharing.rs | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/core/threshold-execution/src/large_execution/single_sharing.rs b/core/threshold-execution/src/large_execution/single_sharing.rs index 570f15b1dd..5ddc78b31b 100644 --- a/core/threshold-execution/src/large_execution/single_sharing.rs +++ b/core/threshold-execution/src/large_execution/single_sharing.rs @@ -1,7 +1,6 @@ use super::local_single_share::{LocalSingleShare, SecureLocalSingleShare}; use crate::runtime::sessions::large_session::LargeSessionHandles; use algebra::{ - bivariate::compute_powers, matrix::MatrixMul, structure_traits::{Derive, ErrorCorrect, Invert, Ring, RingWithExceptionalSequence}, }; @@ -143,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()) } From d424382b80ab888020cbc9c410e59a75fd0745a9 Mon Sep 17 00:00:00 2001 From: David Palm Date: Thu, 7 May 2026 09:18:48 +0200 Subject: [PATCH 09/20] chore: update to new api add: DoublePoly::from_bivariate helper --- .../src/large_execution/vss.rs | 33 +++++++++++-------- .../large_execution/malicious_vss.rs | 26 ++++----------- 2 files changed, 26 insertions(+), 33 deletions(-) diff --git a/core/threshold-execution/src/large_execution/vss.rs b/core/threshold-execution/src/large_execution/vss.rs index e04a33b4cc..cb54be3c98 100644 --- a/core/threshold-execution/src/large_execution/vss.rs +++ b/core/threshold-execution/src/large_execution/vss.rs @@ -106,6 +106,18 @@ where pub(crate) share_in_y: Poly, } +impl DoublePoly { + pub(crate) fn from_bivariate(poly: &BivariatePoly, point: Z) -> Self { + // `partial_evaluations` returns `(F(alpha, Y), F(X, alpha))`, which are the + // recipient's Y-share and X-share respectively. + let (share_in_y, share_in_x) = poly.partial_evaluations(point); + Self { + share_in_x, // F(X, alpha) + share_in_y, // 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 +274,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 +287,7 @@ pub(crate) fn sample_secret_polys>), anyhow::Error>((*r, vec_map)) }) @@ -907,7 +914,7 @@ where vss.my_poly .iter() .map(|poly| poly.full_evaluation(point_x, point_y)) - .collect::, _>>()?, + .collect(), ); } //If im a Pi send Fi(alpha_j) @@ -1031,7 +1038,7 @@ fn round_4_conflict_resolution( vss.my_poly .iter() .map(|poly| poly.partial_y_evaluation(point_pi)) - .collect::, _>>()?, + .collect(), ), //As P_j external from the conflict, resolve conflict with P_i by sending F(alpha_j,alpha_i) false => ValueOrPoly::Value( @@ -1332,7 +1339,7 @@ pub(crate) mod tests { &result .my_poly .iter() - .map(|p| p.full_evaluation(x_0, y_0).unwrap()) + .map(|p| p.full_evaluation(x_0, y_0)) .collect_vec(), expected_secret, ); @@ -1344,12 +1351,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_evaluation(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_evaluation(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)) }) From 5af93188c82e01d8291f49d16d3a37ba6754b081 Mon Sep 17 00:00:00 2001 From: David Palm Date: Thu, 7 May 2026 10:03:35 +0200 Subject: [PATCH 10/20] chore: dial down the benchmarks to more relevant sizes and cases. --- core/threshold-algebra/benches/bivariate.rs | 74 +-------------------- 1 file changed, 2 insertions(+), 72 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 78c3e52a63..d8ef7a39a3 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -8,9 +8,9 @@ use threshold_algebra::{ structure_traits::{Ring, Sample}, }; -const DEGREES: [usize; 5] = [1, 3, 4, 13, 20]; +const DEGREES: [usize; 3] = [1, 4, 13]; // Use this shorter sweep for fast local iteration. -// const DEGREES: [usize; 2] = [4, 20]; +// const DEGREES: [usize; 2] = [4, 13]; fn bivariate_setup(degree: usize) -> (BivariatePoly, ResiduePolyF4Z128) { let mut rng = AesRng::seed_from_u64(degree as u64); @@ -20,19 +20,6 @@ fn bivariate_setup(degree: usize) -> (BivariatePoly, ResidueP (poly, point) } -fn matrix_setup(degree: usize) -> (Vec, Vec) { - let mut rng = AesRng::seed_from_u64(100 + degree as u64); - let d = degree + 1; - let vector = (0..d) - .map(|_| ResiduePolyF4Z128::sample(&mut rng)) - .collect(); - let matrix = (0..d * d) - .map(|_| ResiduePolyF4Z128::sample(&mut rng)) - .collect(); - - (vector, matrix) -} - fn full_evaluation_direct( poly: &BivariatePoly, degree: usize, @@ -53,32 +40,6 @@ fn full_evaluation_direct( acc } -fn vector_matrix_product(vector: &[Z], matrix: &[Z]) -> Vec { - let d = vector.len(); - let mut res = vec![Z::ZERO; d]; - for (row_idx, vector_value) in vector.iter().enumerate() { - let row_start = row_idx * d; - for col_idx in 0..d { - res[col_idx] += *vector_value * matrix[row_start + col_idx]; - } - } - res -} - -fn matrix_vector_product(matrix: &[Z], vector: &[Z]) -> Vec { - let d = vector.len(); - let mut res = Vec::with_capacity(d); - for row_idx in 0..d { - let row_start = row_idx * d; - let mut acc = Z::ZERO; - for col_idx in 0..d { - acc += matrix[row_start + col_idx] * vector[col_idx]; - } - res.push(acc); - } - res -} - fn bench_bivariate_sampling(c: &mut Criterion) { let mut group = c.benchmark_group("bivariate/from_secret"); @@ -117,36 +78,6 @@ fn bench_bivariate_evaluation(c: &mut Criterion) { group.bench_function(BenchmarkId::new("full", degree), |b| { b.iter(|| black_box(poly.full_evaluation(black_box(point), black_box(point)))); }); - group.bench_function(BenchmarkId::new("full_direct", degree), |b| { - b.iter(|| black_box(full_evaluation_direct(&poly, degree, point, point))); - }); - } - - group.finish(); -} - -fn bench_matrix_mul(c: &mut Criterion) { - let mut group = c.benchmark_group("bivariate/matrix_mul"); - - for degree in DEGREES { - let (vector, matrix) = matrix_setup(degree); - - group.bench_function(BenchmarkId::new("vector_matrix", degree), |b| { - b.iter(|| { - black_box(vector_matrix_product( - black_box(&vector), - black_box(&matrix), - )) - }); - }); - group.bench_function(BenchmarkId::new("matrix_vector", degree), |b| { - b.iter(|| { - black_box(matrix_vector_product( - black_box(&matrix), - black_box(&vector), - )) - }); - }); } group.finish(); @@ -156,6 +87,5 @@ criterion_group!( bivariate, bench_bivariate_sampling, bench_bivariate_evaluation, - bench_matrix_mul ); criterion_main!(bivariate); From 0c7b8b42e7569965914553911708ac8af0224edf Mon Sep 17 00:00:00 2001 From: David Palm Date: Fri, 8 May 2026 15:52:28 +0200 Subject: [PATCH 11/20] chore: remove trait, mak`from_secret` the only constructor, remove `Default` derivation --- core/threshold-algebra/benches/bivariate.rs | 29 +--- core/threshold-algebra/src/bivariate.rs | 156 ++++++++---------- .../src/large_execution/vss.rs | 3 +- 3 files changed, 67 insertions(+), 121 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index d8ef7a39a3..1be5276d1d 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -3,14 +3,10 @@ use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; use rand::SeedableRng; use std::hint::black_box; use threshold_algebra::{ - bivariate::{BivariateEval, BivariatePoly, compute_powers}, - galois_rings::degree_4::ResiduePolyF4Z128, - structure_traits::{Ring, Sample}, + bivariate::BivariatePoly, galois_rings::degree_4::ResiduePolyF4Z128, structure_traits::Sample, }; const DEGREES: [usize; 3] = [1, 4, 13]; -// Use this shorter sweep for fast local iteration. -// const DEGREES: [usize; 2] = [4, 13]; fn bivariate_setup(degree: usize) -> (BivariatePoly, ResiduePolyF4Z128) { let mut rng = AesRng::seed_from_u64(degree as u64); @@ -20,26 +16,6 @@ fn bivariate_setup(degree: usize) -> (BivariatePoly, ResidueP (poly, point) } -fn full_evaluation_direct( - poly: &BivariatePoly, - degree: usize, - alpha_x: Z, - alpha_y: Z, -) -> Z { - let powers_x = compute_powers(alpha_x, degree); - let powers_y = compute_powers(alpha_y, degree); - let mut acc = Z::ZERO; - let d = degree + 1; - - for (row_idx, row) in poly.coefs.chunks_exact(d).enumerate() { - for (col_idx, coef) in row.iter().enumerate() { - acc += powers_x[row_idx] * *coef * powers_y[col_idx]; - } - } - - acc -} - fn bench_bivariate_sampling(c: &mut Criterion) { let mut group = c.benchmark_group("bivariate/from_secret"); @@ -65,9 +41,6 @@ fn bench_bivariate_evaluation(c: &mut Criterion) { for degree in DEGREES { let (poly, point) = bivariate_setup(degree); - let current = poly.full_evaluation(point, point); - let direct = full_evaluation_direct(&poly, degree, point, point); - assert_eq!(current, direct); group.bench_function(BenchmarkId::new("partial_x", degree), |b| { b.iter(|| black_box(poly.partial_x_evaluation(black_box(point)))); diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index a33c114ba9..2ecf88551f 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -10,10 +10,10 @@ use std::ops::Mul; /// /// The row view of the polynomials is the following: /// [[a_{00}, a_{01}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}]] -#[derive(Clone, Default, Debug)] +#[derive(Clone, Debug)] pub struct BivariatePoly { /// Row-major coefficients. Contains `(degree + 1)^2` elements. - pub coefs: Vec, + coefs: Vec, degree: usize, } @@ -31,6 +31,29 @@ impl BivariatePoly { BivariatePoly { coefs, degree } } + + #[cfg(test)] + fn from_coeffs(coefs: Vec, degree: usize) -> Self { + let d = degree + 1; + let n = d * d; + assert_eq!( + coefs.len(), + n, + "BivariatePoly expected {n} coefficients for degree {}, got {}", + degree, + coefs.len() + ); + + BivariatePoly { coefs, degree } + } + + fn dim(&self) -> usize { + self.degree + 1 + } + + fn coeffs(&self) -> &[Z] { + &self.coefs + } } /// computes powers of a list of points in F up to a given maximal exponent @@ -55,62 +78,16 @@ pub fn compute_powers + Copy>(point: Z, degree: usize) powers_of_point } -impl BivariatePoly { - fn dimension(&self) -> usize { - self.degree + 1 - } - - fn coefficients(&self) -> &[Z] { - let d = self.dimension(); - let n = d * d; - // @reviewers The old ndarray-based matmul checked ranks and dimensions on - // every evaluation. The only supported bivariate shape is the square coefficient - // matrix implied by `degree`, so we check that invariant once here and let the hot - // loops below use `d` directly. - assert_eq!( - self.coefs.len(), - n, - "BivariatePoly expected {n} coefficients for degree {}, got {}", - self.degree, - self.coefs.len() - ); - &self.coefs - } -} - -pub trait BivariateEval { +impl BivariatePoly { /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the X variable: /// G(Y) = F(\alpha, Y) with coefficients [sum_i a_i0 \alpha^i, ..., sum_i a_id \alpha^i]. - fn partial_x_evaluation(&self, alpha: Z) -> Poly; - - /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the Y variable: - /// G(X) = F(X, \alpha) with coefficients [sum_j a_0j \alpha^j, ..., sum_j a_dj \alpha^j]. - fn partial_y_evaluation(&self, alpha: Z) -> Poly; - - /// Compute both partial evaluations at the same point. - /// - /// Returns `(F(\alpha, Y), F(X, \alpha))`, matching - /// [`Self::partial_x_evaluation`] and [`Self::partial_y_evaluation`]. - fn partial_evaluations(&self, alpha: Z) -> (Poly, Poly) { - ( - self.partial_x_evaluation(alpha), - self.partial_y_evaluation(alpha), - ) - } - - /// 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) -> Z; -} - -impl BivariateEval for BivariatePoly { - fn partial_x_evaluation(&self, alpha: Z) -> Poly { - let d = self.dimension(); - let coefs = self.coefficients(); + pub fn partial_x_evaluation(&self, alpha: Z) -> Poly { + let d = self.dim(); + let coefs = self.coeffs(); let mut res = coefs[(d - 1) * d..d * d].to_vec(); // @reviewers The old implementation checked each column length inside generic - // `matmul`. After the invariant check in `coefficients`, every row has exactly `d` - // coefficients, and `chunks_exact(d)` encodes the only shape this evaluator supports. + // `matmul`. `BivariatePoly` now has private fields and constructors that create or + // validate exactly `d * d` coefficients, so 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; @@ -119,13 +96,15 @@ impl BivariateEval for BivariatePoly { Poly::from_coefs(res) } - fn partial_y_evaluation(&self, alpha: Z) -> Poly { - let d = self.dimension(); - let coefs = self.coefficients(); + /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the Y variable: + /// G(X) = F(X, \alpha) with coefficients [sum_j a_0j \alpha^j, ..., sum_j a_dj \alpha^j]. + pub fn partial_y_evaluation(&self, alpha: Z) -> Poly { + let d = self.dim(); + let coefs = self.coeffs(); let mut res = Vec::with_capacity(d); // @reviewers The previous generic matrix path checked row/vector dimensions before - // multiplication. Here `coefficients` asserts the only valid bivariate shape once, - // and each Horner pass consumes exactly one row of length `d`. + // multiplication. The constructors establish the only valid bivariate shape once, and + // each Horner 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() { @@ -136,9 +115,13 @@ impl BivariateEval for BivariatePoly { Poly::from_coefs(res) } - fn partial_evaluations(&self, alpha: Z) -> (Poly, Poly) { - let d = self.dimension(); - let coefs = self.coefficients(); + /// Compute both partial evaluations at the same point. + /// + /// Returns `(F(\alpha, Y), F(X, \alpha))`, matching + /// [`Self::partial_x_evaluation`] and [`Self::partial_y_evaluation`]. + pub fn partial_evaluations(&self, alpha: Z) -> (Poly, Poly) { + let d = self.dim(); + let coefs = self.coeffs(); let last_row_start = (d - 1) * d; let last_row = &coefs[last_row_start..last_row_start + d]; @@ -151,9 +134,9 @@ impl BivariateEval for BivariatePoly { } partial_y[d - 1] = acc; - // @reviewers This fused path relies on the same `d * d` coefficient invariant as - // the separate evaluators. The only difference is that it computes both directions - // in one pass to avoid repeating the row scan. + // @reviewers This fused path relies on the same constructor-established `d * d` + // coefficient invariant as the separate evaluators. The only difference is that it + // computes both directions in one pass to avoid repeating the row scan. for (row_idx, row) in coefs[..last_row_start].chunks_exact(d).enumerate().rev() { let mut acc = row[d - 1]; for coef in row[..d - 1].iter().rev() { @@ -169,9 +152,11 @@ impl BivariateEval for BivariatePoly { (Poly::from_coefs(partial_x), Poly::from_coefs(partial_y)) } - fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Z { - let d = self.dimension(); - let coefs = self.coefficients(); + /// Given a degree T bivariate poly F(X,Y) and two points \alpha_x, \alpha_y, we compute + /// F(\alpha_x, \alpha_y) + pub fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Z { + let d = self.dim(); + let coefs = self.coeffs(); let mut acc = Z::ZERO; let mut alpha_x_power = Z::ONE; for row in coefs.chunks_exact(d) { @@ -184,7 +169,7 @@ impl BivariateEval for BivariatePoly { } // @reviewers This replaces the old vector-matrix-then-dot-product path. Both power // vectors had length `d` by construction, so the final rank-1 length check bought no - // extra validation beyond the coefficient invariant checked above. + // extra validation beyond the constructor-established coefficient invariant. acc } } @@ -212,10 +197,7 @@ mod tests { #[test] fn test_bivariate_supported_shapes() { let two = ResiduePolyF4Z128::ONE + ResiduePolyF4Z128::ONE; - let bpoly2 = BivariatePoly { - coefs: vec![ResiduePolyF4Z128::ONE; 4], - degree: 1, - }; + let bpoly2 = BivariatePoly::from_coeffs(vec![ResiduePolyF4Z128::ONE; 4], 1); let expected2 = Poly::from_coefs(vec![two; 2]); assert_eq!( bpoly2.partial_x_evaluation(ResiduePolyF4Z128::ONE), @@ -231,10 +213,7 @@ mod tests { ); let five = two + two + ResiduePolyF4Z128::ONE; - let bpoly5 = BivariatePoly { - coefs: vec![ResiduePolyF4Z128::ONE; 25], - degree: 4, - }; + let bpoly5 = BivariatePoly::from_coeffs(vec![ResiduePolyF4Z128::ONE; 25], 4); let expected5 = Poly::from_coefs(vec![five; 5]); assert_eq!( bpoly5.partial_x_evaluation(ResiduePolyF4Z128::ONE), @@ -250,17 +229,6 @@ mod tests { ); } - #[test] - #[should_panic(expected = "BivariatePoly expected 4 coefficients for degree 1")] - fn test_bivariate_invalid_coefficient_count() { - let bpoly = BivariatePoly { - coefs: vec![ResiduePolyF4Z128::ONE; 3], - degree: 1, - }; - - let _ = bpoly.partial_x_evaluation(ResiduePolyF4Z128::ONE); - } - #[rstest] fn test_bivariate_partial_evaluations(#[values(1, 4, 10)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); @@ -301,7 +269,10 @@ mod tests { let secret = ResiduePolyF4Z128::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); let ev_one = bpoly.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE); - let sum_coefs = bpoly.coefs.iter().fold(ResiduePoly::ZERO, |acc, x| acc + x); + let sum_coefs = bpoly + .coeffs() + .iter() + .fold(ResiduePoly::ZERO, |acc, x| acc + x); assert_eq!(ev_one, sum_coefs); } @@ -312,7 +283,10 @@ mod tests { let secret = ResiduePolyF4Z64::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); let ev_one = bpoly.full_evaluation(ResiduePolyF4Z64::ONE, ResiduePolyF4Z64::ONE); - let sum_coefs = bpoly.coefs.iter().fold(ResiduePoly::ZERO, |acc, x| acc + x); + let sum_coefs = bpoly + .coeffs() + .iter() + .fold(ResiduePoly::ZERO, |acc, x| acc + x); assert_eq!(ev_one, sum_coefs); } @@ -622,7 +596,7 @@ mod tests { }, //x4y4 ]; - let bpoly = BivariatePoly { coefs, degree: 4 }; + let bpoly = BivariatePoly::from_coeffs(coefs, 4); let point = ResiduePoly { coefs: [ diff --git a/core/threshold-execution/src/large_execution/vss.rs b/core/threshold-execution/src/large_execution/vss.rs index cb54be3c98..4fed939949 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}, }; @@ -1184,7 +1184,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}, From 8c6444d47e1eceafff841d7186677ea88f525313 Mon Sep 17 00:00:00 2001 From: David Palm Date: Fri, 8 May 2026 20:41:19 +0200 Subject: [PATCH 12/20] chore: tighter comments&docs chore: moved compute_powers* to matrix --- core/threshold-algebra/benches/bivariate.rs | 14 +++ core/threshold-algebra/src/bivariate.rs | 94 +++++++-------- .../src/galois_rings/common.rs | 2 +- .../src/galois_rings/degree_3.rs | 2 +- .../src/galois_rings/degree_4.rs | 2 +- .../src/galois_rings/degree_5.rs | 2 +- .../src/galois_rings/degree_6.rs | 2 +- .../src/galois_rings/degree_7.rs | 2 +- .../src/galois_rings/degree_8.rs | 2 +- core/threshold-algebra/src/matrix.rs | 111 +++++++++++++++++- core/threshold-algebra/src/syndrome.rs | 2 +- .../src/small_execution/prss.rs | 3 +- 12 files changed, 181 insertions(+), 57 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 1be5276d1d..cefaf754bd 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -48,6 +48,20 @@ fn bench_bivariate_evaluation(c: &mut Criterion) { group.bench_function(BenchmarkId::new("partial_y", degree), |b| { b.iter(|| black_box(poly.partial_y_evaluation(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_evaluation(p)), + black_box(poly.partial_y_evaluation(p)), + ) + }); + }); + // Fused: the path `DoublePoly::from_bivariate` actually uses. + group.bench_function(BenchmarkId::new("partial_evaluations", degree), |b| { + b.iter(|| black_box(poly.partial_evaluations(black_box(point)))); + }); group.bench_function(BenchmarkId::new("full", degree), |b| { b.iter(|| black_box(poly.full_evaluation(black_box(point), black_box(point)))); }); diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 2ecf88551f..bb6f68d807 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -1,18 +1,19 @@ use super::poly::Poly; -use super::structure_traits::One; use super::structure_traits::Ring; use super::structure_traits::Sample; use rand::{CryptoRng, Rng}; -use std::ops::Mul; /// Bivariate polynomial represented as a row-major square matrix of coefficients. /// /// The row view of the polynomials is the following: /// [[a_{00}, a_{01}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}]] +/// +/// Constructors establish and maintain the invariant `coefs.len() == (degree + 1)^2`, +/// so every row in `coefs.chunks_exact(degree + 1)` has exactly `degree + 1` entries. #[derive(Clone, Debug)] pub struct BivariatePoly { - /// Row-major coefficients. Contains `(degree + 1)^2` elements. + /// Row-major; `(degree + 1)^2` elements. coefs: Vec, degree: usize, } @@ -35,15 +36,7 @@ impl BivariatePoly { #[cfg(test)] fn from_coeffs(coefs: Vec, degree: usize) -> Self { let d = degree + 1; - let n = d * d; - assert_eq!( - coefs.len(), - n, - "BivariatePoly expected {n} coefficients for degree {}, got {}", - degree, - coefs.len() - ); - + assert_eq!(coefs.len(), d * d); BivariatePoly { coefs, degree } } @@ -56,28 +49,6 @@ impl BivariatePoly { } } -/// 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 -} - -/// 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); // start with - for i in 1..=degree { - powers_of_point.push(powers_of_point[i - 1] * point); - } - powers_of_point -} - impl BivariatePoly { /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the X variable: /// G(Y) = F(\alpha, Y) with coefficients [sum_i a_i0 \alpha^i, ..., sum_i a_id \alpha^i]. @@ -85,9 +56,7 @@ impl BivariatePoly { let d = self.dim(); let coefs = self.coeffs(); let mut res = coefs[(d - 1) * d..d * d].to_vec(); - // @reviewers The old implementation checked each column length inside generic - // `matmul`. `BivariatePoly` now has private fields and constructors that create or - // validate exactly `d * d` coefficients, so every row has exactly `d` coefficients. + // 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; @@ -102,9 +71,7 @@ impl BivariatePoly { let d = self.dim(); let coefs = self.coeffs(); let mut res = Vec::with_capacity(d); - // @reviewers The previous generic matrix path checked row/vector dimensions before - // multiplication. The constructors establish the only valid bivariate shape once, and - // each Horner pass consumes exactly one row of length `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() { @@ -134,9 +101,7 @@ impl BivariatePoly { } partial_y[d - 1] = acc; - // @reviewers This fused path relies on the same constructor-established `d * d` - // coefficient invariant as the separate evaluators. The only difference is that it - // computes both directions in one pass to avoid repeating the row scan. + // Computes both directions in one pass to avoid repeating the row scan. for (row_idx, row) in coefs[..last_row_start].chunks_exact(d).enumerate().rev() { let mut acc = row[d - 1]; for coef in row[..d - 1].iter().rev() { @@ -167,9 +132,6 @@ impl BivariatePoly { acc += alpha_x_power * row_acc; alpha_x_power *= alpha_x; } - // @reviewers This replaces the old vector-matrix-then-dot-product path. Both power - // vectors had length `d` by construction, so the final rank-1 length check bought no - // extra validation beyond the constructor-established coefficient invariant. acc } } @@ -229,6 +191,46 @@ mod tests { ); } + #[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 = one + one; + let three = two + one; + let four = two + two; + let five = four + one; + let seven = four + three; + 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 = five + five + five + one; + let twentytwo = sixteen + three + three; + assert_eq!( + bpoly.partial_x_evaluation(five), + Poly::from_coefs(vec![sixteen, twentytwo]) + ); + + // F(X, 5) = (1 + 2*5) + (3 + 4*5)*X = 11 + 23*X + let eleven = five + five + one; + let twentythree = eleven + four + four + four; + assert_eq!( + bpoly.partial_y_evaluation(five), + Poly::from_coefs(vec![eleven, twentythree]) + ); + + // Fused returns (F(alpha, Y), F(X, alpha)). + let (px, py) = bpoly.partial_evaluations(five); + assert_eq!(px, Poly::from_coefs(vec![sixteen, twentytwo])); + assert_eq!(py, Poly::from_coefs(vec![eleven, twentythree])); + + // F(5, 7) = 1 + 2*7 + 3*5 + 4*5*7 = 170 + let mut expected_full = ResiduePolyF4Z128::ZERO; + for _ in 0..170 { + expected_full += one; + } + assert_eq!(bpoly.full_evaluation(five, seven), expected_full); + } + #[rstest] fn test_bivariate_partial_evaluations(#[values(1, 4, 10)] degree: usize) { let mut rng = AesRng::seed_from_u64(0); 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 92900bee36..591305e358 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 e558e9696c..d628c0dcdf 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 81ea035785..e211ca71b4 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 df006e34ca..a1915ca818 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 cfd464561d..4925e44a00 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 8190867b33..d4a69eb6a0 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/matrix.rs b/core/threshold-algebra/src/matrix.rs index 47809c2e36..bef220e234 100644 --- a/core/threshold-algebra/src/matrix.rs +++ b/core/threshold-algebra/src/matrix.rs @@ -1,15 +1,38 @@ -use crate::structure_traits::Ring; +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; @@ -94,3 +117,89 @@ impl MatrixMul> for ArrayD { } } } + +#[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 784ba2e4f3..0499434443 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/small_execution/prss.rs b/core/threshold-execution/src/small_execution/prss.rs index 2d6e578c36..78c97e4920 100644 --- a/core/threshold-execution/src/small_execution/prss.rs +++ b/core/threshold-execution/src/small_execution/prss.rs @@ -21,8 +21,7 @@ use crate::{ }, }; use algebra::{ - bivariate::compute_powers_list, - matrix::MatrixMul, + matrix::{MatrixMul, compute_powers_list}, poly::Poly, structure_traits::{ErrorCorrect, Invert, Ring, RingWithExceptionalSequence}, }; From 6f6347963a4c6830b58f38d385bce048537e2e04 Mon Sep 17 00:00:00 2001 From: David Palm Date: Fri, 8 May 2026 21:23:06 +0200 Subject: [PATCH 13/20] chore: shorter names opt: do "outer Horner" for full_eval opt: see if avoiding indexing is faster (avoids bounds checks) --- core/threshold-algebra/benches/bivariate.rs | 12 +- core/threshold-algebra/src/bivariate.rs | 107 +++++++++--------- .../src/large_execution/vss.rs | 12 +- 3 files changed, 66 insertions(+), 65 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index cefaf754bd..08175e4e44 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -43,27 +43,27 @@ fn bench_bivariate_evaluation(c: &mut Criterion) { let (poly, point) = bivariate_setup(degree); group.bench_function(BenchmarkId::new("partial_x", degree), |b| { - b.iter(|| black_box(poly.partial_x_evaluation(black_box(point)))); + 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_evaluation(black_box(point)))); + 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_evaluation(p)), - black_box(poly.partial_y_evaluation(p)), + black_box(poly.partial_x_eval(p)), + black_box(poly.partial_y_eval(p)), ) }); }); // Fused: the path `DoublePoly::from_bivariate` actually uses. group.bench_function(BenchmarkId::new("partial_evaluations", degree), |b| { - b.iter(|| black_box(poly.partial_evaluations(black_box(point)))); + b.iter(|| black_box(poly.partial_evals(black_box(point)))); }); group.bench_function(BenchmarkId::new("full", degree), |b| { - b.iter(|| black_box(poly.full_evaluation(black_box(point), black_box(point)))); + b.iter(|| black_box(poly.full_eval(black_box(point), black_box(point)))); }); } diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index bb6f68d807..f871eb0848 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -52,7 +52,7 @@ impl BivariatePoly { impl BivariatePoly { /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the X variable: /// G(Y) = F(\alpha, Y) with coefficients [sum_i a_i0 \alpha^i, ..., sum_i a_id \alpha^i]. - pub fn partial_x_evaluation(&self, alpha: Z) -> Poly { + 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(); @@ -67,7 +67,7 @@ impl BivariatePoly { /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the Y variable: /// G(X) = F(X, \alpha) with coefficients [sum_j a_0j \alpha^j, ..., sum_j a_dj \alpha^j]. - pub fn partial_y_evaluation(&self, alpha: Z) -> Poly { + pub fn partial_y_eval(&self, alpha: Z) -> Poly { let d = self.dim(); let coefs = self.coeffs(); let mut res = Vec::with_capacity(d); @@ -85,33 +85,43 @@ impl BivariatePoly { /// Compute both partial evaluations at the same point. /// /// Returns `(F(\alpha, Y), F(X, \alpha))`, matching - /// [`Self::partial_x_evaluation`] and [`Self::partial_y_evaluation`]. - pub fn partial_evaluations(&self, alpha: Z) -> (Poly, Poly) { + /// [`Self::partial_x_eval`] and [`Self::partial_y_eval`]. + // TODO(dp): Not sure this is actually worth keeping, The perf gain from doing this vs partial_x + partial_y is ~3% which seems silly. + pub fn partial_evals(&self, alpha: Z) -> (Poly, Poly) { let d = self.dim(); let coefs = self.coeffs(); let last_row_start = (d - 1) * d; let last_row = &coefs[last_row_start..last_row_start + d]; + // partial_x is seeded with the last row (the high-i term in the outer Horner). let mut partial_x = last_row.to_vec(); let mut partial_y = vec![Z::ZERO; d]; + // Last row's partial_y is its own Horner; its partial_x slots are already correct. let mut acc = last_row[d - 1]; for coef in last_row[..d - 1].iter().rev() { acc = acc * alpha + *coef; } partial_y[d - 1] = acc; - // Computes both directions in one pass to avoid repeating the row scan. - for (row_idx, row) in coefs[..last_row_start].chunks_exact(d).enumerate().rev() { - let mut acc = row[d - 1]; - for coef in row[..d - 1].iter().rev() { - acc = acc * alpha + *coef; - } - partial_y[row_idx] = acc; - - for (res_coef, coef) in partial_x.iter_mut().zip(row) { - *res_coef = *res_coef * alpha + *coef; + // Walk each remaining row once. Each coefficient feeds two Horner chains: + // the row-wise chain producing the next `partial_y` coefficient, and the + // column-wise chain updating `partial_x`. This keeps the fused path equivalent + // to two separate partial evaluations while reusing each loaded coefficient. + for (partial_y_coef, row) in partial_y[..d - 1] + .iter_mut() + .rev() + .zip(coefs[..last_row_start].chunks_exact(d).rev()) + { + let last = row[d - 1]; + let mut acc = last; + partial_x[d - 1] = partial_x[d - 1] * alpha + last; + for (res_coef, coef) in partial_x[..d - 1].iter_mut().zip(&row[..d - 1]).rev() { + let c = *coef; + acc = acc * alpha + c; + *res_coef = *res_coef * alpha + c; } + *partial_y_coef = acc; } (Poly::from_coefs(partial_x), Poly::from_coefs(partial_y)) @@ -119,18 +129,21 @@ impl BivariatePoly { /// Given a degree T bivariate poly F(X,Y) and two points \alpha_x, \alpha_y, we compute /// F(\alpha_x, \alpha_y) - pub fn full_evaluation(&self, alpha_x: Z, alpha_y: Z) -> Z { + pub fn full_eval(&self, alpha_x: Z, alpha_y: Z) -> Z { let d = self.dim(); let coefs = self.coeffs(); - let mut acc = Z::ZERO; - let mut alpha_x_power = Z::ONE; - for row in coefs.chunks_exact(d) { + 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 += alpha_x_power * row_acc; - alpha_x_power *= alpha_x; + acc = acc * alpha_x + row_acc; } acc } @@ -161,32 +174,20 @@ mod tests { let two = ResiduePolyF4Z128::ONE + ResiduePolyF4Z128::ONE; 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.partial_x_evaluation(ResiduePolyF4Z128::ONE), - expected2 - ); - assert_eq!( - bpoly2.partial_y_evaluation(ResiduePolyF4Z128::ONE), - expected2 - ); - assert_eq!( - bpoly2.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), + bpoly2.full_eval(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), two + two ); let five = two + two + ResiduePolyF4Z128::ONE; 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.partial_x_evaluation(ResiduePolyF4Z128::ONE), - expected5 - ); - assert_eq!( - bpoly5.partial_y_evaluation(ResiduePolyF4Z128::ONE), - expected5 - ); - assert_eq!( - bpoly5.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), + bpoly5.full_eval(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE), five + five + five + five + five ); } @@ -206,7 +207,7 @@ mod tests { let sixteen = five + five + five + one; let twentytwo = sixteen + three + three; assert_eq!( - bpoly.partial_x_evaluation(five), + bpoly.partial_x_eval(five), Poly::from_coefs(vec![sixteen, twentytwo]) ); @@ -214,12 +215,12 @@ mod tests { let eleven = five + five + one; let twentythree = eleven + four + four + four; assert_eq!( - bpoly.partial_y_evaluation(five), + bpoly.partial_y_eval(five), Poly::from_coefs(vec![eleven, twentythree]) ); // Fused returns (F(alpha, Y), F(X, alpha)). - let (px, py) = bpoly.partial_evaluations(five); + let (px, py) = bpoly.partial_evals(five); assert_eq!(px, Poly::from_coefs(vec![sixteen, twentytwo])); assert_eq!(py, Poly::from_coefs(vec![eleven, twentythree])); @@ -228,7 +229,7 @@ mod tests { for _ in 0..170 { expected_full += one; } - assert_eq!(bpoly.full_evaluation(five, seven), expected_full); + assert_eq!(bpoly.full_eval(five, seven), expected_full); } #[rstest] @@ -238,10 +239,10 @@ mod tests { let point = ResiduePolyF4Z128::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - let (partial_x, partial_y) = bpoly.partial_evaluations(point); + let (partial_x, partial_y) = bpoly.partial_evals(point); - assert_eq!(partial_x, bpoly.partial_x_evaluation(point)); - assert_eq!(partial_y, bpoly.partial_y_evaluation(point)); + assert_eq!(partial_x, bpoly.partial_x_eval(point)); + assert_eq!(partial_y, bpoly.partial_y_eval(point)); } //Test that eval at 0 return the secret for ResiduePolyF4Z128 @@ -250,7 +251,7 @@ mod tests { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z128::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - let ev_zero = bpoly.full_evaluation(ResiduePolyF4Z128::ZERO, ResiduePolyF4Z128::ZERO); + let ev_zero = bpoly.full_eval(ResiduePolyF4Z128::ZERO, ResiduePolyF4Z128::ZERO); assert_eq!(ev_zero, secret); } @@ -260,7 +261,7 @@ mod tests { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z64::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - let ev_zero = bpoly.full_evaluation(ResiduePolyF4Z64::ZERO, ResiduePolyF4Z64::ZERO); + let ev_zero = bpoly.full_eval(ResiduePolyF4Z64::ZERO, ResiduePolyF4Z64::ZERO); assert_eq!(ev_zero, secret); } @@ -270,7 +271,7 @@ mod tests { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z128::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - let ev_one = bpoly.full_evaluation(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE); + let ev_one = bpoly.full_eval(ResiduePolyF4Z128::ONE, ResiduePolyF4Z128::ONE); let sum_coefs = bpoly .coeffs() .iter() @@ -284,7 +285,7 @@ mod tests { let mut rng = AesRng::seed_from_u64(0); let secret = ResiduePolyF4Z64::sample(&mut rng); let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - let ev_one = bpoly.full_evaluation(ResiduePolyF4Z64::ONE, ResiduePolyF4Z64::ONE); + let ev_one = bpoly.full_eval(ResiduePolyF4Z64::ONE, ResiduePolyF4Z64::ONE); let sum_coefs = bpoly .coeffs() .iter() @@ -621,7 +622,7 @@ mod tests { #[test] fn test_bivariate_partial_eval_x() { let (bpoly, point) = poly_setup(); - let res = bpoly.partial_x_evaluation(point); + let res = bpoly.partial_x_eval(point); let expected_result = Poly::::from_coefs(vec![ ResiduePoly { @@ -695,7 +696,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); + let res = bpoly.partial_y_eval(point); let expected_result = Poly::::from_coefs(vec![ ResiduePoly { @@ -768,9 +769,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); + let res = bpoly.full_eval(point_x, point_y); - let expected_res = bpoly.partial_x_evaluation(point_x).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-execution/src/large_execution/vss.rs b/core/threshold-execution/src/large_execution/vss.rs index 4fed939949..23b1eeaf60 100644 --- a/core/threshold-execution/src/large_execution/vss.rs +++ b/core/threshold-execution/src/large_execution/vss.rs @@ -110,7 +110,7 @@ impl DoublePoly { pub(crate) fn from_bivariate(poly: &BivariatePoly, point: Z) -> Self { // `partial_evaluations` returns `(F(alpha, Y), F(X, alpha))`, which are the // recipient's Y-share and X-share respectively. - let (share_in_y, share_in_x) = poly.partial_evaluations(point); + let (share_in_y, share_in_x) = poly.partial_evals(point); Self { share_in_x, // F(X, alpha) share_in_y, // F(alpha, Y) @@ -913,7 +913,7 @@ where (*own_role, *pi_role, *pj_role), vss.my_poly .iter() - .map(|poly| poly.full_evaluation(point_x, point_y)) + .map(|poly| poly.full_eval(point_x, point_y)) .collect(), ); } @@ -1037,7 +1037,7 @@ fn round_4_conflict_resolution( true => ValueOrPoly::Poly( vss.my_poly .iter() - .map(|poly| poly.partial_y_evaluation(point_pi)) + .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) @@ -1338,7 +1338,7 @@ pub(crate) mod tests { &result .my_poly .iter() - .map(|p| p.full_evaluation(x_0, y_0)) + .map(|p| p.full_eval(x_0, y_0)) .collect_vec(), expected_secret, ); @@ -1350,12 +1350,12 @@ pub(crate) mod tests { let expected_result_x = result .my_poly .iter() - .map(|p| p.partial_y_evaluation(embedded_pn)) + .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)) + .map(|p| p.partial_x_eval(embedded_pn)) .collect_vec(); assert_eq!( From 2f889a322b182e7bcca594d7fa93593fd72e636f Mon Sep 17 00:00:00 2001 From: David Palm Date: Sun, 10 May 2026 22:28:41 +0200 Subject: [PATCH 14/20] chore: tweak docs a wee bit --- core/threshold-algebra/src/bivariate.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index f871eb0848..3392e44722 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -4,13 +4,12 @@ use super::structure_traits::Sample; use rand::{CryptoRng, Rng}; -/// Bivariate polynomial represented as a row-major square matrix of coefficients. +/// Bivariate polynomial represented as a row-major square matrix of coefficients of ResiduePolynomials. /// /// The row view of the polynomials is the following: /// [[a_{00}, a_{01}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}]] /// -/// Constructors establish and maintain the invariant `coefs.len() == (degree + 1)^2`, -/// so every row in `coefs.chunks_exact(degree + 1)` has exactly `degree + 1` entries. +/// Invariant: `coefs.len() == (degree + 1)^2` #[derive(Clone, Debug)] pub struct BivariatePoly { /// Row-major; `(degree + 1)^2` elements. @@ -33,6 +32,7 @@ impl BivariatePoly { BivariatePoly { coefs, degree } } + // 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; From 3d87b2274a3aa6dcfff98e1fccf4312c94c1b30f Mon Sep 17 00:00:00 2001 From: David Palm Date: Sun, 10 May 2026 22:41:18 +0200 Subject: [PATCH 15/20] fix: did ToB rename the example we use? --- dylint.toml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/dylint.toml b/dylint.toml index 4fb9048b20..64d55a0eae 100644 --- a/dylint.toml +++ b/dylint.toml @@ -1,8 +1,7 @@ [workspace.metadata.dylint] libraries = [ - { git = "https://github.com/trailofbits/dylint", pattern = "examples/general/non_local_effect_before_error_return" }, + { git = "https://github.com/trailofbits/dylint", pattern = "examples/general/non_local_effect_before_unhandled_error" }, ] -[non_local_effect_before_error_return] +[non_local_effect_before_unhandled_error] work_limit = 100000000 - From 9839cf075b01e70dc1d1da6ac316b46ae6b59ccf Mon Sep 17 00:00:00 2001 From: David Palm Date: Sun, 10 May 2026 23:21:50 +0200 Subject: [PATCH 16/20] chore: fix dylint --- dylint.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dylint.toml b/dylint.toml index 64d55a0eae..81bbf2da9c 100644 --- a/dylint.toml +++ b/dylint.toml @@ -1,7 +1,7 @@ [workspace.metadata.dylint] libraries = [ - { git = "https://github.com/trailofbits/dylint", pattern = "examples/general/non_local_effect_before_unhandled_error" }, + { git = "https://github.com/trailofbits/dylint", rev = "2c20d164e06c446de2cecadff0c55e16a877c006", pattern = "examples/general/non_local_effect_before_error_return" }, ] -[non_local_effect_before_unhandled_error] +[non_local_effect_before_error_return] work_limit = 100000000 From e0e880558e1d8db685a6c321e830c360d3a641e8 Mon Sep 17 00:00:00 2001 From: David Date: Tue, 12 May 2026 17:42:21 +0200 Subject: [PATCH 17/20] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- core/threshold-execution/src/large_execution/vss.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/threshold-execution/src/large_execution/vss.rs b/core/threshold-execution/src/large_execution/vss.rs index 23b1eeaf60..32ca394a0f 100644 --- a/core/threshold-execution/src/large_execution/vss.rs +++ b/core/threshold-execution/src/large_execution/vss.rs @@ -108,7 +108,7 @@ where impl DoublePoly { pub(crate) fn from_bivariate(poly: &BivariatePoly, point: Z) -> Self { - // `partial_evaluations` returns `(F(alpha, Y), F(X, alpha))`, which are the + // `partial_evals` returns `(F(alpha, Y), F(X, alpha))`, which are the // recipient's Y-share and X-share respectively. let (share_in_y, share_in_x) = poly.partial_evals(point); Self { From 3f694d651af2343b99414d9e88ca4df48cbea29d Mon Sep 17 00:00:00 2001 From: David Palm Date: Fri, 22 May 2026 12:26:15 +0200 Subject: [PATCH 18/20] chore: address review feedback --- core/threshold-algebra/src/bivariate.rs | 46 ++++++++++++++----------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 0331debd2f..9c18a2e7e1 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -50,8 +50,12 @@ impl BivariatePoly { } impl BivariatePoly { - /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the X variable: - /// G(Y) = F(\alpha, Y) with coefficients [sum_i a_i0 \alpha^i, ..., sum_i a_id \alpha^i]. + /// 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(); @@ -65,8 +69,12 @@ impl BivariatePoly { Poly::from_coefs(res) } - /// Given a degree T bivariate poly F(X,Y) = sum a_ij X^i Y^j and a point \alpha, evaluate the Y variable: - /// G(X) = F(X, \alpha) with coefficients [sum_j a_0j \alpha^j, ..., sum_j a_dj \alpha^j]. + /// 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(); @@ -153,6 +161,7 @@ impl BivariatePoly { mod tests { use super::*; use crate::galois_rings::degree_8::ResiduePolyF8Z128; + use crate::structure_traits::{FromU128, ZConsts}; use crate::{ galois_rings::{ common::ResiduePoly, @@ -169,7 +178,7 @@ mod tests { //Checks the hot coefficient shapes used by bivariate evaluation. #[test] fn test_bivariate_supported_shapes() { - let two = ResiduePolyF4Z128::ONE + ResiduePolyF4Z128::ONE; + 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); @@ -179,14 +188,14 @@ mod tests { two + two ); - let five = two + two + ResiduePolyF4Z128::ONE; + 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), - five + five + five + five + five + ResiduePolyF4Z128::from_u128(25) ); } @@ -194,24 +203,24 @@ mod tests { 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 = one + one; - let three = two + one; - let four = two + two; - let five = four + one; - let seven = four + three; + 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 = five + five + five + one; - let twentytwo = sixteen + three + three; + 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 = five + five + one; - let twentythree = eleven + four + four + four; + 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]) @@ -223,10 +232,7 @@ mod tests { assert_eq!(py, Poly::from_coefs(vec![eleven, twentythree])); // F(5, 7) = 1 + 2*7 + 3*5 + 4*5*7 = 170 - let mut expected_full = ResiduePolyF4Z128::ZERO; - for _ in 0..170 { - expected_full += one; - } + let expected_full = ResiduePolyF4Z128::from_u128(170); assert_eq!(bpoly.full_eval(five, seven), expected_full); } From 671f134b828f47fc17527b8197f4e0daface0c71 Mon Sep 17 00:00:00 2001 From: David Palm Date: Fri, 22 May 2026 12:49:58 +0200 Subject: [PATCH 19/20] =?UTF-8?q?Remove=20partial=5Fevals=20=E2=80=93=20no?= =?UTF-8?q?t=20a=20good=20enough=20speedup?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core/threshold-algebra/benches/bivariate.rs | 4 -- core/threshold-algebra/src/bivariate.rs | 66 +------------------ .../src/large_execution/vss.rs | 9 ++- 3 files changed, 5 insertions(+), 74 deletions(-) diff --git a/core/threshold-algebra/benches/bivariate.rs b/core/threshold-algebra/benches/bivariate.rs index 08175e4e44..3c3a681b3b 100644 --- a/core/threshold-algebra/benches/bivariate.rs +++ b/core/threshold-algebra/benches/bivariate.rs @@ -58,10 +58,6 @@ fn bench_bivariate_evaluation(c: &mut Criterion) { ) }); }); - // Fused: the path `DoublePoly::from_bivariate` actually uses. - group.bench_function(BenchmarkId::new("partial_evaluations", degree), |b| { - b.iter(|| black_box(poly.partial_evals(black_box(point)))); - }); group.bench_function(BenchmarkId::new("full", degree), |b| { b.iter(|| black_box(poly.full_eval(black_box(point), black_box(point)))); }); diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 9c18a2e7e1..91443f5f88 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -90,53 +90,7 @@ impl BivariatePoly { Poly::from_coefs(res) } - /// Compute both partial evaluations at the same point. - /// - /// Returns `(F(\alpha, Y), F(X, \alpha))`, matching - /// [`Self::partial_x_eval`] and [`Self::partial_y_eval`]. - // TODO(dp): Not sure this is actually worth keeping, The perf gain from doing this vs partial_x + partial_y is ~3% which seems silly. - pub fn partial_evals(&self, alpha: Z) -> (Poly, Poly) { - let d = self.dim(); - let coefs = self.coeffs(); - let last_row_start = (d - 1) * d; - let last_row = &coefs[last_row_start..last_row_start + d]; - - // partial_x is seeded with the last row (the high-i term in the outer Horner). - let mut partial_x = last_row.to_vec(); - let mut partial_y = vec![Z::ZERO; d]; - - // Last row's partial_y is its own Horner; its partial_x slots are already correct. - let mut acc = last_row[d - 1]; - for coef in last_row[..d - 1].iter().rev() { - acc = acc * alpha + *coef; - } - partial_y[d - 1] = acc; - - // Walk each remaining row once. Each coefficient feeds two Horner chains: - // the row-wise chain producing the next `partial_y` coefficient, and the - // column-wise chain updating `partial_x`. This keeps the fused path equivalent - // to two separate partial evaluations while reusing each loaded coefficient. - for (partial_y_coef, row) in partial_y[..d - 1] - .iter_mut() - .rev() - .zip(coefs[..last_row_start].chunks_exact(d).rev()) - { - let last = row[d - 1]; - let mut acc = last; - partial_x[d - 1] = partial_x[d - 1] * alpha + last; - for (res_coef, coef) in partial_x[..d - 1].iter_mut().zip(&row[..d - 1]).rev() { - let c = *coef; - acc = acc * alpha + c; - *res_coef = *res_coef * alpha + c; - } - *partial_y_coef = acc; - } - - (Poly::from_coefs(partial_x), Poly::from_coefs(partial_y)) - } - - /// Given a degree T bivariate poly F(X,Y) and two points \alpha_x, \alpha_y, we compute - /// F(\alpha_x, \alpha_y) + /// 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(); @@ -226,29 +180,11 @@ mod tests { Poly::from_coefs(vec![eleven, twentythree]) ); - // Fused returns (F(alpha, Y), F(X, alpha)). - let (px, py) = bpoly.partial_evals(five); - assert_eq!(px, Poly::from_coefs(vec![sixteen, twentytwo])); - assert_eq!(py, 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); } - #[rstest] - fn test_bivariate_partial_evaluations(#[values(1, 4, 10)] degree: usize) { - let mut rng = AesRng::seed_from_u64(0); - let secret = ResiduePolyF4Z128::sample(&mut rng); - let point = ResiduePolyF4Z128::sample(&mut rng); - let bpoly = BivariatePoly::from_secret(&mut rng, secret, degree); - - let (partial_x, partial_y) = bpoly.partial_evals(point); - - assert_eq!(partial_x, bpoly.partial_x_eval(point)); - assert_eq!(partial_y, bpoly.partial_y_eval(point)); - } - //Test that eval at 0 return the secret for ResiduePolyF4Z128 #[rstest] fn test_bivariate_zero_128(#[values(4, 10, 20)] degree: usize) { diff --git a/core/threshold-execution/src/large_execution/vss.rs b/core/threshold-execution/src/large_execution/vss.rs index 626a0fa3e9..7b686ea3fd 100644 --- a/core/threshold-execution/src/large_execution/vss.rs +++ b/core/threshold-execution/src/large_execution/vss.rs @@ -102,18 +102,17 @@ 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 { - // `partial_evals` returns `(F(alpha, Y), F(X, alpha))`, which are the - // recipient's Y-share and X-share respectively. - let (share_in_y, share_in_x) = poly.partial_evals(point); Self { - share_in_x, // F(X, alpha) - share_in_y, // F(alpha, Y) + share_in_x: poly.partial_y_eval(point), // F(X, alpha) + share_in_y: poly.partial_x_eval(point), // F(alpha, Y) } } } From 0b911836143eaab8979b7a550163a497b2ac5415 Mon Sep 17 00:00:00 2001 From: David Palm Date: Fri, 22 May 2026 13:05:13 +0200 Subject: [PATCH 20/20] Update doc-comment --- core/threshold-algebra/src/bivariate.rs | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/core/threshold-algebra/src/bivariate.rs b/core/threshold-algebra/src/bivariate.rs index 91443f5f88..5471657f2e 100644 --- a/core/threshold-algebra/src/bivariate.rs +++ b/core/threshold-algebra/src/bivariate.rs @@ -4,12 +4,14 @@ use super::structure_traits::Sample; use rand::{CryptoRng, Rng}; -/// Bivariate polynomial represented as a row-major square matrix of coefficients of ResiduePolynomials. +/// Bivariate polynomial represented as a square, row-major coefficient matrix. /// -/// The row view of the polynomials is the following: -/// [[a_{00}, a_{01}, ..., a_{0d}], ..., [a_{d0}, ..., a_{dd}]] +/// The row view of the coefficients is: +/// [[a_{00}, a_{01}, ..., a_{0t}], ..., [a_{t0}, ..., a_{tt}]] /// -/// Invariant: `coefs.len() == (degree + 1)^2` +/// 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 { /// Row-major; `(degree + 1)^2` elements.