From 035d15bc3d3598b07cba92ec78d1b5be1ee734e6 Mon Sep 17 00:00:00 2001 From: Daniel Lubarov Date: Fri, 23 Apr 2021 10:26:57 -0700 Subject: [PATCH] Interpolants of arbitrary (point, value) lists Closes #10. This combines Lagrange interpolation with FFTs as mentioned there. I was previously thinking that all our polynomial encodings might as well just use power-of-two length vectors, so they'll be "FFT-ready", with no need to trim/pad. This sort of breaks that assumption though, as e.g. I think we'll want to compute interpolants with three coefficients in the batch opening argument. I think we can still skip trimming/padding in most cases, since it the majority of our polynomials will have power-of-two-minus-1 degrees with high probability. But we'll now have one or two uses where that's not the case. --- src/field/crandall_field.rs | 13 ++++ src/field/fft.rs | 2 +- src/field/field.rs | 21 +++++-- src/field/lagrange.rs | 111 +++++++++++++++++++++++++++++++++++ src/field/mod.rs | 1 + src/polynomial/polynomial.rs | 23 ++++++-- 6 files changed, 158 insertions(+), 13 deletions(-) create mode 100644 src/field/lagrange.rs diff --git a/src/field/crandall_field.rs b/src/field/crandall_field.rs index 0ec184e8..657193d9 100644 --- a/src/field/crandall_field.rs +++ b/src/field/crandall_field.rs @@ -6,6 +6,7 @@ use num::Integer; use crate::field::field::Field; use std::hash::{Hash, Hasher}; +use std::iter::{Product, Sum}; /// EPSILON = 9 * 2**28 - 1 const EPSILON: u64 = 2415919103; @@ -257,6 +258,12 @@ impl AddAssign for CrandallField { } } +impl Sum for CrandallField { + fn sum>(iter: I) -> Self { + iter.fold(Self::ZERO, |acc, x| acc + x) + } +} + impl Sub for CrandallField { type Output = Self; @@ -291,6 +298,12 @@ impl MulAssign for CrandallField { } } +impl Product for CrandallField { + fn product>(iter: I) -> Self { + iter.fold(Self::ONE, |acc, x| acc * x) + } +} + impl Div for CrandallField { type Output = Self; diff --git a/src/field/fft.rs b/src/field/fft.rs index f38ea204..8bcde967 100644 --- a/src/field/fft.rs +++ b/src/field/fft.rs @@ -169,7 +169,7 @@ mod tests { for i in 0..degree { coefficients.push(F::from_canonical_usize(i * 1337 % 100)); } - let coefficients = PolynomialCoeffs::pad(coefficients); + let coefficients = PolynomialCoeffs::new_padded(coefficients); let points = fft(coefficients.clone()); assert_eq!(points, evaluate_naive(&coefficients)); diff --git a/src/field/field.rs b/src/field/field.rs index 6c12bf7d..15a1a9f4 100644 --- a/src/field/field.rs +++ b/src/field/field.rs @@ -1,10 +1,13 @@ -use crate::util::bits_u64; -use rand::rngs::OsRng; -use rand::Rng; use std::fmt::{Debug, Display}; use std::hash::Hash; +use std::iter::{Product, Sum}; use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign}; +use rand::rngs::OsRng; +use rand::Rng; + +use crate::util::bits_u64; + /// A finite field with prime order less than 2^64. pub trait Field: 'static @@ -14,10 +17,12 @@ pub trait Field: + Neg + Add + AddAssign + + Sum + Sub + SubAssign + Mul + MulAssign + + Product + Div + DivAssign + Debug @@ -42,6 +47,10 @@ pub trait Field: *self == Self::ZERO } + fn is_nonzero(&self) -> bool { + *self != Self::ZERO + } + fn is_one(&self) -> bool { *self == Self::ONE } @@ -90,12 +99,12 @@ pub trait Field: x_inv } - fn primitive_root_of_unity(n_power: usize) -> Self { - assert!(n_power <= Self::TWO_ADICITY); + fn primitive_root_of_unity(n_log: usize) -> Self { + assert!(n_log <= Self::TWO_ADICITY); let base = Self::POWER_OF_TWO_GENERATOR; // TODO: Just repeated squaring should be a bit faster, to avoid conditionals. base.exp(Self::from_canonical_u64( - 1u64 << (Self::TWO_ADICITY - n_power), + 1u64 << (Self::TWO_ADICITY - n_log), )) } diff --git a/src/field/lagrange.rs b/src/field/lagrange.rs new file mode 100644 index 00000000..36241da7 --- /dev/null +++ b/src/field/lagrange.rs @@ -0,0 +1,111 @@ +use crate::field::fft::ifft; +use crate::field::field::Field; +use crate::polynomial::polynomial::{PolynomialCoeffs, PolynomialValues}; +use crate::util::log2_ceil; + +/// Computes the interpolant of an arbitrary list of (point, value) pairs. +/// +/// Note that the implementation assumes that `F` is two-adic, in particular that +/// `2^{F::TWO_ADICITY} >= points.len()`. This leads to a simple FFT-based implementation. +pub(crate) fn interpolant(points: &[(F, F)]) -> PolynomialCoeffs { + let n = points.len(); + let n_log = log2_ceil(n); + let n_padded = 1 << n_log; + + let g = F::primitive_root_of_unity(n_log); + let subgroup = F::cyclic_subgroup_known_order(g, n_padded); + let subgroup_evals = subgroup + .into_iter() + .map(|x| interpolate(points, x)) + .collect(); + + let mut coeffs = ifft(PolynomialValues { + values: subgroup_evals, + }); + coeffs.trim(); + coeffs +} + +/// Interpolate the polynomial defined by an arbitrary set of (point, value) pairs at the given +/// point `x`. +fn interpolate(points: &[(F, F)], x: F) -> F { + (0..points.len()) + .map(|i| { + let y_i = points[i].1; + let l_i_x = eval_basis(points, i, x); + y_i * l_i_x + }) + .sum() +} + +/// Evaluate the `i`th Lagrange basis, i.e. the one that vanishes except on the `i`th point. +fn eval_basis(points: &[(F, F)], i: usize, x: F) -> F { + let n = points.len(); + let x_i = points[i].0; + let mut numerator = F::ONE; + let mut denominator_parts = Vec::with_capacity(n - 1); + + for j in 0..n { + if i != j { + let x_j = points[j].0; + numerator *= x - x_j; + denominator_parts.push(x_i - x_j); + } + } + + let denominator_inv = F::batch_multiplicative_inverse(&denominator_parts) + .into_iter() + .product(); + numerator * denominator_inv +} + +#[cfg(test)] +mod tests { + use crate::field::crandall_field::CrandallField; + use crate::field::field::Field; + use crate::field::lagrange::interpolant; + use crate::polynomial::polynomial::PolynomialCoeffs; + + #[test] + fn interpolant_random() { + type F = CrandallField; + + for deg in 0..10 { + let domain = (0..deg).map(|_| F::rand()).collect::>(); + let coeffs = (0..deg).map(|_| F::rand()).collect(); + let coeffs = PolynomialCoeffs { coeffs }; + + let points = eval_naive(&coeffs, &domain); + assert_eq!(interpolant(&points), coeffs); + } + } + + #[test] + fn interpolant_random_overspecified() { + type F = CrandallField; + + for deg in 0..10 { + let points = deg + 5; + let domain = (0..points).map(|_| F::rand()).collect::>(); + let coeffs = (0..deg).map(|_| F::rand()).collect(); + let coeffs = PolynomialCoeffs { coeffs }; + + let points = eval_naive(&coeffs, &domain); + assert_eq!(interpolant(&points), coeffs); + } + } + + fn eval_naive(coeffs: &PolynomialCoeffs, domain: &[F]) -> Vec<(F, F)> { + domain + .iter() + .map(|&x| { + let eval = x + .powers() + .zip(&coeffs.coeffs) + .map(|(x_power, &coeff)| coeff * x_power) + .sum(); + (x, eval) + }) + .collect() + } +} diff --git a/src/field/mod.rs b/src/field/mod.rs index 9f58ef08..6c2828a6 100644 --- a/src/field/mod.rs +++ b/src/field/mod.rs @@ -2,6 +2,7 @@ pub(crate) mod cosets; pub mod crandall_field; pub mod fft; pub mod field; +pub(crate) mod lagrange; #[cfg(test)] mod field_testing; diff --git a/src/polynomial/polynomial.rs b/src/polynomial/polynomial.rs index 0444b1e8..7034a8b9 100644 --- a/src/polynomial/polynomial.rs +++ b/src/polynomial/polynomial.rs @@ -2,7 +2,7 @@ use crate::field::fft::{fft, ifft}; use crate::field::field::Field; use crate::util::log2_strict; -/// A polynomial in point-value form. The number of values must be a power of two. +/// A polynomial in point-value form. /// /// The points are implicitly `g^i`, where `g` generates the subgroup whose size equals the number /// of points. @@ -13,7 +13,6 @@ pub struct PolynomialValues { impl PolynomialValues { pub fn new(values: Vec) -> Self { - assert!(values.len().is_power_of_two()); PolynomialValues { values } } @@ -36,7 +35,7 @@ impl PolynomialValues { } } -/// A polynomial in coefficient form. The number of coefficients must be a power of two. +/// A polynomial in coefficient form. #[derive(Clone, Debug, Eq, PartialEq)] pub struct PolynomialCoeffs { pub(crate) coeffs: Vec, @@ -44,11 +43,11 @@ pub struct PolynomialCoeffs { impl PolynomialCoeffs { pub fn new(coeffs: Vec) -> Self { - assert!(coeffs.len().is_power_of_two()); PolynomialCoeffs { coeffs } } - pub(crate) fn pad(mut coeffs: Vec) -> Self { + /// Create a new polynomial with its coefficient list padded to the next power of two. + pub(crate) fn new_padded(mut coeffs: Vec) -> Self { while !coeffs.len().is_power_of_two() { coeffs.push(F::ZERO); } @@ -70,7 +69,6 @@ impl PolynomialCoeffs { } pub(crate) fn chunks(&self, chunk_size: usize) -> Vec { - assert!(chunk_size.is_power_of_two()); self.coeffs .chunks(chunk_size) .map(|chunk| PolynomialCoeffs::new(chunk.to_vec())) @@ -97,4 +95,17 @@ impl PolynomialCoeffs { } Self { coeffs } } + + /// Removes leading zero coefficients. + pub fn trim(&mut self) { + self.coeffs.drain(self.degree_plus_one()..); + } + + /// Degree of the polynomial + 1. + fn degree_plus_one(&self) -> usize { + (0usize..self.len()) + .rev() + .find(|&i| self.coeffs[i].is_nonzero()) + .map_or(0, |i| i + 1) + } }