binius_math/
univariate.rs

1// Copyright 2024-2025 Irreducible Inc.
2
3use binius_field::{BinaryField, Field};
4
5use super::BinarySubspace;
6
7/// Evaluate a univariate polynomial specified by its monomial coefficients.
8///
9/// # Arguments
10/// * `coeffs` - Slice of coefficients ordered from low-degree terms to high-degree terms
11/// * `x` - Point at which to evaluate the polynomial
12pub fn evaluate_univariate<F: Field>(coeffs: &[F], x: F) -> F {
13	let Some((&highest_degree, rest)) = coeffs.split_last() else {
14		return F::ZERO;
15	};
16
17	// Evaluate using Horner's method
18	rest.iter()
19		.rev()
20		.fold(highest_degree, |acc, &coeff| acc * x + coeff)
21}
22
23/// Optimized Lagrange evaluation for power-of-2 domains in binary fields.
24///
25/// Computes the Lagrange polynomial evaluations L̃(z, i) for a power-of-2 domain at point `z`.
26/// Uses the provided binary subspace as the evaluation domain.
27///
28/// # Key Optimization
29/// For power-of-2 domains, all barycentric weights are identical due to the additive group
30/// structure. For each i ∈ {0, ..., 2^k - 1}, the set {i ⊕ j | j ≠ i} = {1, ..., 2^k - 1}.
31/// This allows us to:
32/// 1. Compute a single barycentric weight w = 1 / ∏_{j=1}^{n-1} j
33/// 2. Use prefix/suffix products to avoid redundant computation
34/// 3. Replace inversions with multiplications for better performance
35///
36/// # Complexity
37/// - Time: O(n) where n = subspace size, using 4n - 2 multiplications and 1 inversion
38/// - Space: O(n) for prefix/suffix arrays
39///
40/// # Parameters
41/// - `subspace`: The binary subspace defining the evaluation domain
42/// - `z`: The evaluation point
43///
44/// # Returns
45/// A vector of Lagrange polynomial evaluations, one for each domain element
46pub fn lagrange_evals<F: BinaryField>(subspace: &BinarySubspace<F>, z: F) -> Vec<F> {
47	let domain: Vec<F> = subspace.iter().collect();
48	let n = domain.len();
49
50	// Compute single barycentric weight for the additive subgroup
51	// All points have the same weight due to subgroup structure
52	let w = domain[1..]
53		.iter()
54		.fold(F::ONE, |acc, &d| acc * d)
55		.invert()
56		.unwrap_or(F::ONE);
57
58	// Compute prefix products: prefix[i] = ∏_{j=0}^{i-1} (z - domain[j])
59	let mut prefixes = vec![F::ONE; n];
60	for i in 1..n {
61		prefixes[i] = prefixes[i - 1] * (z - domain[i - 1]);
62	}
63
64	// Compute suffix products: suffix[i] = ∏_{j=i+1}^{n-1} (z - domain[j])
65	let mut suffixes = vec![F::ONE; n];
66	for i in (0..n - 1).rev() {
67		suffixes[i] = suffixes[i + 1] * (z - domain[i + 1]);
68	}
69
70	// Combine prefix, suffix, and weight: L_i(z) = prefix[i] * suffix[i] * w
71	let mut result = vec![F::ZERO; n];
72	for i in 0..n {
73		result[i] = prefixes[i] * suffixes[i] * w;
74	}
75
76	result
77}
78
79#[cfg(test)]
80mod tests {
81	use binius_field::{BinaryField128bGhash, Field, Random, util::powers};
82	use rand::prelude::*;
83
84	use super::*;
85	use crate::{BinarySubspace, inner_product::inner_product, test_utils::random_scalars};
86
87	fn evaluate_univariate_with_powers<F: Field>(coeffs: &[F], x: F) -> F {
88		inner_product(coeffs.iter().copied(), powers(x).take(coeffs.len()))
89	}
90
91	type F = BinaryField128bGhash;
92
93	#[test]
94	fn test_evaluate_univariate_against_reference() {
95		let mut rng = StdRng::seed_from_u64(0);
96
97		for n_coeffs in [0, 1, 2, 5, 10] {
98			let coeffs = random_scalars(&mut rng, n_coeffs);
99			let x = F::random(&mut rng);
100			assert_eq!(
101				evaluate_univariate(&coeffs, x),
102				evaluate_univariate_with_powers(&coeffs, x)
103			);
104		}
105	}
106
107	#[test]
108	fn test_lagrange_evals() {
109		let mut rng = StdRng::seed_from_u64(0);
110
111		// Test mathematical properties across different domain sizes
112		for log_domain_size in [3, 4, 5, 6] {
113			// Create subspace for this test
114			let subspace = BinarySubspace::<F>::with_dim(log_domain_size).unwrap();
115			let domain: Vec<F> = subspace.iter().collect();
116
117			// Test 1: Partition of Unity - Lagrange polynomials sum to 1
118			let eval_point = F::random(&mut rng);
119			let lagrange_coeffs = lagrange_evals(&subspace, eval_point);
120			let sum: F = lagrange_coeffs.iter().copied().sum();
121			assert_eq!(
122				sum,
123				F::ONE,
124				"Partition of unity failed for domain size {}",
125				1 << log_domain_size
126			);
127
128			// Test 2: Interpolation Property - L_i(x_j) = δ_ij
129			for (j, &domain_point) in domain.iter().enumerate() {
130				let lagrange_at_domain = lagrange_evals(&subspace, domain_point);
131				for (i, &coeff) in lagrange_at_domain.iter().enumerate() {
132					let expected = if i == j { F::ONE } else { F::ZERO };
133					assert_eq!(
134						coeff, expected,
135						"Interpolation property failed: L_{i}({j}) ≠ {expected}"
136					);
137				}
138			}
139		}
140
141		// Test 3: Polynomial Interpolation Accuracy
142		let log_domain_size = 6;
143		let subspace = BinarySubspace::<F>::with_dim(log_domain_size).unwrap();
144		let domain: Vec<F> = subspace.iter().collect();
145		let coeffs = random_scalars(&mut rng, 10);
146
147		// Evaluate polynomial at domain points
148		let domain_evals: Vec<F> = domain
149			.iter()
150			.map(|&point| evaluate_univariate(&coeffs, point))
151			.collect();
152
153		// Test interpolation at random point
154		let test_point = F::random(&mut rng);
155		let lagrange_coeffs = lagrange_evals(&subspace, test_point);
156		let interpolated =
157			inner_product(domain_evals.iter().copied(), lagrange_coeffs.iter().copied());
158		let direct = evaluate_univariate(&coeffs, test_point);
159
160		assert_eq!(interpolated, direct, "Polynomial interpolation accuracy failed");
161	}
162}
binius_math/univariate.rs

binius_math/
univariate.rs
