From eea6702ac41701f33c7835a1e05f9c7a31115726 Mon Sep 17 00:00:00 2001 From: Colin Roberts Date: Thu, 9 May 2024 10:49:03 -0600 Subject: [PATCH] feat: basic polynomial arithmetic (#48) * feat: basic polynomial arithmetic * fix: remove polynomial const generics --- math/field.sage | 22 ++- math/polynomial.sage | 36 ++++ src/field/mod.rs | 10 +- src/lib.rs | 7 + src/polynomial.rs | 217 ---------------------- src/polynomial/arithmetic.rs | 253 ++++++++++++++++++++++++++ src/polynomial/mod.rs | 340 +++++++++++++++++++++++++++++++++++ 7 files changed, 663 insertions(+), 222 deletions(-) create mode 100644 math/polynomial.sage delete mode 100644 src/polynomial.rs create mode 100644 src/polynomial/arithmetic.rs create mode 100644 src/polynomial/mod.rs diff --git a/math/field.sage b/math/field.sage index ad696947..afea94da 100644 --- a/math/field.sage +++ b/math/field.sage @@ -1,4 +1,4 @@ -# Define a finite field, our `PlutoField` which is +# Define a finite field, our `GF101` which is # a field of 101 elements (101 being prime). F = GF(101) print(F) @@ -39,6 +39,26 @@ omega_n = primitive_element ^ quotient print("The", n, "th root of unity is: ", omega_n) # Check that this is actually a root of unity: +assert omega_n^n == 1 +print(omega_n, "^", n, " = ", omega_n^n) +###################################################################### + +###################################################################### +# Let's find a mth root of unity (for l = 2) +# First, check that m divides 101 - 1 = 100 +l = 2 +assert (101 - 1) % l == 0 +quotient = (101 - 1) // l +print("The quotient is: ", quotient) + +# Find a primitive root of unity using the formula: +# omega = primitive_element^quotient +omega_l = primitive_element^quotient +print("The ", l, "th root of unity is: ", omega_l) + +# Check that this is actually a root of unity: +assert omega_l^l == 1 +print(omega_l, "^", l, " = ", omega_l^l) assert omega_n ^ n == 1 print(omega_n, "^", n, " = ", omega_n ^ n) ###################################################################### diff --git a/math/polynomial.sage b/math/polynomial.sage new file mode 100644 index 00000000..02d38ea6 --- /dev/null +++ b/math/polynomial.sage @@ -0,0 +1,36 @@ +# Define a finite field, our `GF101` which is +# a field of 101 elements (101 being prime). +F = GF(101) + +# Define the polynomial ring over GF(101) +R. = PolynomialRing(F) + +# Create the polynomials +a = 4*x^3 + 3*x^2 + 2*x^1 + 1 +b = 9*x^4 + 8*x^3 + 7*x^2 + 6*x^1 + 5 + +# Perform polynomial division +q_ab, r_ab = a.quo_rem(b) + +# Print the results +print("For a / b :") +print("Quotient:", q_ab) +print("Remainder:", r_ab) + +# Perform polynomial division +q_ba, r_ba = b.quo_rem(a) + +# Print the results +print("\nFor b / a :") +print("Quotient:", q_ba) +print("Remainder:", r_ba) + +# Polynomials easy to do by hand +p = x^2 + 2*x + 1 +q = x + 1 + +# Perform polynomial division +quot, rem = p.quo_rem(q) +print("\nFor p / q :") +print("Quotient:", quot) +print("Remainder:", rem) \ No newline at end of file diff --git a/src/field/mod.rs b/src/field/mod.rs index 30e7c5fc..9eb15f03 100644 --- a/src/field/mod.rs +++ b/src/field/mod.rs @@ -4,6 +4,8 @@ use std::{ ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, Sub, SubAssign}, }; +use super::*; + pub mod gf_101; pub mod gf_101_2; @@ -74,10 +76,10 @@ pub trait FiniteField: // - There exists a multiplicative subgroup generated by a primitive element 'a'. // // According to the Sylow theorems (https://en.wikipedia.org/wiki/Sylow_theorems): - // A non-trivial multiplicative subgroup of prime order 'q' exists if and only if - // 'p - 1' is divisible by 'q'. - // The primitive q-th root of unity 'w' is defined as: w = a^((p - 1) / q), - // and the roots of unity are generated by 'w', such that {w^i | i in [0, q - 1]}. + // A non-trivial multiplicative subgroup of prime order 'n' exists if and only if + // 'p - 1' is divisible by 'n'. + // The primitive n-th root of unity 'w' is defined as: w = a^((p - 1) / n), + // and the roots of unity are generated by 'w', such that {w^i | i in [0, n - 1]}. fn primitive_root_of_unity(n: Self::Storage) -> Self { let p_minus_one = Self::ORDER - Self::Storage::from(1); assert!(p_minus_one % n == 0.into(), "n must divide p - 1"); diff --git a/src/lib.rs b/src/lib.rs index 341159f6..e8f83366 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -12,3 +12,10 @@ pub mod curves; pub mod field; pub mod kzg; pub mod polynomial; + +use core::{ + fmt, + hash::Hash, + iter::{Product, Sum}, + ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, Sub, SubAssign}, +}; diff --git a/src/polynomial.rs b/src/polynomial.rs deleted file mode 100644 index bb796d74..00000000 --- a/src/polynomial.rs +++ /dev/null @@ -1,217 +0,0 @@ -use std::{collections::HashSet, ops::Div}; - -use crate::field::FiniteField; - -// https://people.inf.ethz.ch/gander/papers/changing.pdf - -/// A polynomial of degree N-1 -#[derive(Debug, Clone, PartialEq, Eq, Hash)] -pub struct Polynomial { - pub coefficients: [F; N], - pub basis: B, -} - -// Type state patern for polynomial bases: -pub trait Basis { - type Data; -} - -pub struct Monomial; -impl Basis for Monomial { - type Data = (); -} - -// https://en.wikipedia.org/wiki/Lagrange_polynomial -pub struct Lagrange { - pub nodes: [F; N], - pub weights: [F; N], - pub L: Box F>, -} - -impl Basis for Lagrange { - type Data = Self; -} - -#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] -pub struct Nodes([F; N]); - -/// A polynomial in monomial basis -impl Polynomial { - // TODO check that the polynomial degree divides the field order - pub fn new(coefficients: [F; N]) -> Self { Self { coefficients, basis: Monomial } } - - /// Convert a polynomial from monomial to Lagrange basis using the Barycentric form - pub fn to_lagrange(&self, nodes: Nodes) -> Polynomial, F> { - // check that nodes are distinct - let mut set = HashSet::new(); - for &node in nodes.0.iter() { - if !set.insert(node) { - panic!("Nodes must be distinct"); - } - } - let mut coeffs = [F::ONE; N]; - for j in 0..N { - coeffs[j] = self.evaluate(nodes.0[j]); - } - Polynomial::, F>::new(coeffs, nodes.0) - } - - /// Evaluate the polynomial at field element x - pub fn evaluate(&self, x: F) -> F { - let mut result = F::ONE; - for (i, c) in self.coefficients.into_iter().enumerate() { - result += c * x.pow(F::Storage::from(i as u32)); - } - result - } -} - -/// A polynomial in Lagrange basis -impl Polynomial, F> { - // TODO check that the polynomial degree divides the field order - pub fn new(coefficients: [F; N], nodes: [F; N]) -> Self { - let mut weights = [F::ONE; N]; - for j in 0..N { - for m in 0..N { - if j != m { - weights[j] *= F::ONE.div(nodes[j] - nodes[m]); - } - } - } - // l(x) = \Sigmm_{m}(x-x_m) - let l = move |x: F| { - let mut val = F::ONE; - for i in 0..N { - val *= x - nodes[i]; - } - val - }; - - // L(x) = l(x) * \Sigma_{j=0}^{k} (w_j / (x - x_j)) y_j - let L = move |x: F| { - let mut val = F::ONE; - for j in 0..N { - // check if we are dividing by zero - if nodes[j] == x { - return coefficients[j]; - } - val += coefficients[j] * weights[j] / (x - nodes[j]); - } - l(x) * val - }; - - let basis = Lagrange { nodes, weights, L: Box::new(L) }; - Self { coefficients, basis } - } - - pub fn to_monomial(&self) -> Polynomial { - // This is the inverse of the conversion from monomial to Lagrange basis - // This uses something called the Vandermonde matrix which is defined as: - // - // / 1 | x_0 | x_0^2 | x_0^3 | ... | x_0^(N-1) \ - // | 1 | x_1 | x_1^2 | x_1^3 | ... | x_1^(N-1) | - // | 1 | x_2 | x_2^2 | x_2^3 | ... | x_2^(N-1) | - // v = | . | . | . | . | ... | . | - // | . | . | . | . | ... | . | - // | . | . | . | . | ... | . | - // \ 1 | x_N | x_N^2 | x_N^3 | ... | x_N^(N-1) / - // - // where x_i are the nodes of the Lagrange basis - // - // Then the monomial basis m is given V^T * l = m, where l is the Lagrange basis - // because we know the monomial basis we need to compute to monomial coefficients a_m = V^{-1} * - // a_l where a_l are the coefficients of the Lagrange basis - - // It also is the case that the the columns of the inverse matrix are the coefficients of the - // Lagrange polynomial basis TODO Finish this. - // let nodes = self.basis.nodes; - // let mut evaluations = [F::ONE; N]; - - // // Evaluate the polynomial at N distinct points - // for i in 0..N { - // let x = F::primitive_root().exp_u64(i as u64); - // evaluations[i] = self.evaluate(x); - // } - - // Polynomial::::new(evaluations) - todo!("Finish this after we get the roots of unity from other PRs") - } - - /// Evaluate the polynomial at field element x - pub fn evaluate(&self, x: F) -> F { - let L = &self.basis.L; - L(x) - } -} -mod test { - use super::*; - use crate::field::gf_101::GF101; - - #[test] - fn evaluation() { - // Coefficients of the polynomial 1 + 2x + 3x^2 - let a = GF101::new(1); - let b = GF101::new(2); - let c = GF101::new(3); - let polynomial = Polynomial::<3, Monomial, GF101>::new([a, b, c]); - let y = polynomial.evaluate(GF101::new(2)); - let r = GF101::new(17); - assert_eq!(y, r); // 1 + 2*(2) + 3(2)^2 = 17 - } - - #[test] - fn evaluation_with_zero() { - // Coefficients of the polynomial 1 + 3x^2 - let a = GF101::new(1); - let b = GF101::new(0); - let c = GF101::new(3); - let polynomial = Polynomial::<3, Monomial, GF101>::new([a, b, c]); - let y = polynomial.evaluate(GF101::new(0)); - let r = GF101::new(1); - assert_eq!(y, r); // 1 + 3(0)^2 = 1 - } - - #[test] - fn lagrange_evaluation() { - // Coefficients of the polynomial 1 + 2x + 3x^2 - let a = GF101::new(1); - let b = GF101::new(2); - let c = GF101::new(3); - let polynomial = Polynomial::<3, Monomial, GF101>::new([a, b, c]); - - // Nodes for the Lagrange basis - let nodes = Nodes::<3, GF101>([GF101::new(1), GF101::new(2), GF101::new(3)]); - let lagrange = polynomial.to_lagrange(nodes); - let r = lagrange.evaluate(GF101::new(2)); - assert_eq!(r, GF101::new(17)); - } - - #[test] - #[should_panic] - fn non_unique_nodes() { - // Coefficients of the polynomial 1 + 2x - let a = GF101::new(1); - let b = GF101::new(2); - - let polynomial = Polynomial::<2, Monomial, GF101>::new([a, b]); - // This should panic because the nodes are not distinct - let nodes = Nodes::<2, GF101>([GF101::new(1), GF101::new(1)]); - let lagrange = polynomial.to_lagrange(nodes); - } - - #[test] - fn test_by_hand() { - // Coefficients of the polynomial 1 + x + 2x^2 - let a = GF101::new(1); - let b = GF101::new(1); - let c = GF101::new(2); - let polynomial = Polynomial::<3, Monomial, GF101>::new([a, b, c]); - println!("monomial coefficients {:?}", polynomial.coefficients); - - // Nodes for the Lagrange basis - let nodes = Nodes::<3, GF101>([GF101::new(1), GF101::new(3), GF101::new(2)]); - - let lagrange = polynomial.to_lagrange(nodes); - println!("lagrange coefficients {:?}", lagrange.coefficients); - } -} diff --git a/src/polynomial/arithmetic.rs b/src/polynomial/arithmetic.rs new file mode 100644 index 00000000..39ab27f5 --- /dev/null +++ b/src/polynomial/arithmetic.rs @@ -0,0 +1,253 @@ +use std::array; + +use super::*; + +// TODO: Implement `Mul`, `MulAssign`, `Product`. +// * `Mul` and other products will require FFT. + +impl Add for Polynomial { + type Output = Self; + + fn add(self, rhs: Self) -> Self { + let d = self.degree().max(rhs.degree()); + let mut coefficients = vec![F::ZERO; d + 1]; + for i in 0..d + 1 { + coefficients[i] = *self.coefficients.get(i).unwrap_or(&F::ZERO) + + *rhs.coefficients.get(i).unwrap_or(&F::ZERO); + } + Self { coefficients, basis: self.basis } + } +} + +impl AddAssign for Polynomial { + fn add_assign(&mut self, mut rhs: Self) { + let d = self.degree().max(rhs.degree()); + let mut coefficients = vec![F::ZERO; d + 1]; + if self.degree() < d { + self.coefficients.resize(d + 1, F::ZERO); + } else { + rhs.coefficients.resize(d + 1, F::ZERO); + } + for i in 0..d + 1 { + self.coefficients[i] += rhs.coefficients[i]; + } + } +} + +impl Sum for Polynomial { + fn sum>(iter: I) -> Self { iter.reduce(|x, y| x + y).unwrap() } +} + +impl Sub for Polynomial { + type Output = Self; + + fn sub(self, rhs: Self) -> Self { + let d = self.degree().max(rhs.degree()); + let mut coefficients = vec![F::ZERO; d + 1]; + for i in 0..d + 1 { + coefficients[i] = *self.coefficients.get(i).unwrap_or(&F::ZERO) + - *rhs.coefficients.get(i).unwrap_or(&F::ZERO); + } + Self { coefficients, basis: self.basis } + } +} + +impl SubAssign for Polynomial { + fn sub_assign(&mut self, mut rhs: Self) { + let d = self.degree().max(rhs.degree()); + let mut coefficients = vec![F::ZERO; d + 1]; + if self.degree() < d { + self.coefficients.resize(d + 1, F::ZERO); + } else { + rhs.coefficients.resize(d + 1, F::ZERO); + } + for i in 0..d + 1 { + self.coefficients[i] -= rhs.coefficients[i]; + } + } +} + +impl Neg for Polynomial { + type Output = Self; + + fn neg(self) -> Self { + Self { + coefficients: self.coefficients.into_iter().map(|c| -c).collect(), + basis: self.basis, + } + } +} + +impl Div for Polynomial { + type Output = Self; + + fn div(self, rhs: Self) -> Self::Output { + // // Euclidean division + // // Initial quotient value + // let mut q = F::zero(); + + // // Initial remainder value is our numerator polynomial + // let mut p = self; + + // // Leading coefficient of the denominator + // let c = rhs.leading_coefficient(); + + // // Create quotient poly + // let mut diff = p.degree() as isize - rhs.degree() as isize; + // if diff < 0 { + // return Polynomial::::new(vec![F::zero(); 1]); + // } + // let mut q_coeffs = vec![F::zero(); diff as usize + 1]; + + // while diff >= 0 { + // let s = p.leading_coefficient() * c.inverse().unwrap(); + // q_coeffs[diff as usize] = s; + // p -= rhs.pow_mult(s, diff as usize); + // p.trim_zeros(); + // diff -= 1; + // } + // Polynomial::::new(q_coeffs) + self.quotient_and_remainder(rhs).0 + } +} + +impl Rem for Polynomial { + type Output = Self; + + fn rem(self, rhs: Self) -> Self { self.quotient_and_remainder(rhs).1 } +} + +#[cfg(test)] +mod tests { + use field::gf_101::GF101; + + use super::*; + + fn poly_a() -> Polynomial { + Polynomial::::new(vec![ + GF101::new(1), + GF101::new(2), + GF101::new(3), + GF101::new(4), + ]) + } + + fn poly_b() -> Polynomial { + Polynomial::::new(vec![ + GF101::new(5), + GF101::new(6), + GF101::new(7), + GF101::new(8), + GF101::new(9), + ]) + } + + #[test] + fn add() { + let c = poly_a() + poly_b(); + println!("{:?}", c.coefficients); + assert_eq!(c.coefficients, [ + GF101::new(6), + GF101::new(8), + GF101::new(10), + GF101::new(12), + GF101::new(9) + ]); + } + + #[test] + fn add_assign() { + let mut a = poly_a(); + a += poly_b(); + assert_eq!(a.coefficients, [ + GF101::new(6), + GF101::new(8), + GF101::new(10), + GF101::new(12), + GF101::new(9) + ]); + } + + #[test] + fn sum() { + let sum = [poly_a(), poly_b()].into_iter().sum::>(); + assert_eq!(sum.coefficients, [ + GF101::new(6), + GF101::new(8), + GF101::new(10), + GF101::new(12), + GF101::new(9) + ]); + } + + #[test] + fn sub() { + let c = poly_a() - poly_b(); + assert_eq!(c.coefficients, [ + GF101::new(97), + GF101::new(97), + GF101::new(97), + GF101::new(97), + GF101::new(92) + ]); + } + + #[test] + fn sub_assign() { + let mut a = poly_a(); + a -= poly_b(); + assert_eq!(a.coefficients, [ + GF101::new(97), + GF101::new(97), + GF101::new(97), + GF101::new(97), + GF101::new(92) + ]); + } + + #[test] + fn neg() { + let a = -poly_a(); + assert_eq!(a.coefficients, [GF101::new(100), GF101::new(99), GF101::new(98), GF101::new(97)]); + } + + #[test] + fn div() { + println!("poly_a: {}", poly_a()); + println!("poly_b: {}", poly_b()); + let q_ab = poly_a() / poly_b(); + assert_eq!(q_ab.coefficients, [GF101::new(0)]); + + let q_ba = poly_b() / poly_a(); + assert_eq!(q_ba.coefficients, [GF101::new(95), GF101::new(78)]); + println!("q_ba coefficients: {:?}", q_ba.coefficients); + + let p = Polynomial::::new(vec![GF101::new(1), GF101::new(2), GF101::new(1)]); + println!("p: {}", p); + let q = Polynomial::::new(vec![GF101::new(1), GF101::new(1)]); + println!("q: {}", q); + let r = p / q; + println!("r: {}", r); + assert_eq!(r.coefficients, [GF101::new(1), GF101::new(1)]); + } + + #[test] + fn rem() { + println!("poly_a: {}", poly_a()); + println!("poly_b: {}", poly_b()); + let q_ab = poly_a() % poly_b(); + assert_eq!(q_ab.coefficients, poly_a().coefficients); + + let q_ba = poly_b() % poly_a(); + assert_eq!(q_ba.coefficients, [GF101::new(11), GF101::new(41), GF101::new(71)]); + println!("q_ba coefficients: {:?}", q_ba.coefficients); + + let p = Polynomial::::new(vec![GF101::new(1), GF101::new(2), GF101::new(1)]); + println!("p: {}", p); + let q = Polynomial::::new(vec![GF101::new(1), GF101::new(1)]); + println!("q: {}", q); + let r = p % q; + println!("r: {}", r); + assert_eq!(r.coefficients, [GF101::new(0)]); + } +} diff --git a/src/polynomial/mod.rs b/src/polynomial/mod.rs new file mode 100644 index 00000000..71aaabef --- /dev/null +++ b/src/polynomial/mod.rs @@ -0,0 +1,340 @@ +use std::{ + collections::HashSet, + fmt::{Display, Formatter}, +}; + +use self::field::gf_101::GF101; +use super::*; +use crate::field::FiniteField; + +pub mod arithmetic; + +// https://people.inf.ethz.ch/gander/papers/changing.pdf + +/// A polynomial of degree N-1 which is generic over the basis and the field. +#[derive(Debug, Clone, PartialEq, Eq, Hash)] +pub struct Polynomial { + pub coefficients: Vec, + pub basis: B, +} + +// Type state patern for polynomial bases: +pub trait Basis { + type Data; +} + +#[derive(Debug, Clone, PartialEq, Eq, Hash)] +pub struct Monomial; +impl Basis for Monomial { + type Data = (); +} + +// https://en.wikipedia.org/wiki/Lagrange_polynomial +pub struct Lagrange { + pub nodes: Vec, +} + +impl Basis for Lagrange { + type Data = Self; +} + +impl Polynomial { + pub fn degree(&self) -> usize { self.coefficients.len() - 1 } +} + +/// A polynomial in monomial basis +impl Polynomial { + pub fn new(mut coefficients: Vec) -> Self { + // Simplify the polynomial to have leading coefficient non-zero + let mut poly = Self { coefficients, basis: Monomial }; + poly.trim_zeros(); + poly + } + + fn trim_zeros(&mut self) { + let last_nonzero_index = self.coefficients.iter().rposition(|&c| c != F::ZERO); + match last_nonzero_index { + Some(index) => self.coefficients.truncate(index + 1), + None => self.coefficients.truncate(1), + } + } + + /// Convert a polynomial from monomial to Lagrange basis using the Barycentric form + pub fn to_lagrange(&self) -> Polynomial, F> { + // Get the `N-1`th roots of unity for the field and this degree of polynomial. + let n = self.coefficients.len(); + let primitive_root = F::primitive_root_of_unity(F::Storage::from(n as u32)); + + // Evaluate the polynomial at the roots of unity to get the coefficients of the Lagrange basis + let mut coeffs = vec![F::ZERO; n]; + for j in 0..self.coefficients.len() { + coeffs[j] = self.evaluate(primitive_root.pow(F::Storage::from(j as u32))); + } + Polynomial::, F>::new(coeffs) + } + + /// Evaluate the polynomial at field element x + pub fn evaluate(&self, x: F) -> F { + let mut result = F::ZERO; + for (i, c) in self.coefficients.iter().enumerate() { + result += *c * x.pow(F::Storage::from(i as u32)); + } + result + } + + pub fn leading_coefficient(&self) -> F { *self.coefficients.last().unwrap() } + + pub fn pow_mult(&self, coeff: F, pow: usize) -> Polynomial { + let mut coefficients = vec![F::ZERO; self.coefficients.len() + pow]; + self.coefficients.iter().enumerate().for_each(|(i, c)| { + coefficients[i + pow] = *c * coeff; + }); + Polynomial::::new(coefficients) + } + + fn quotient_and_remainder(self, rhs: Self) -> (Self, Self) { + // Euclidean division + // Initial quotient value + let mut q = Self::new(vec![]); + + // Initial remainder value is our numerator polynomial + let mut p = self.clone(); + + // Leading coefficient of the denominator + let c = rhs.leading_coefficient(); + + // Create quotient poly + let mut diff = p.degree() as isize - rhs.degree() as isize; + if diff < 0 { + return (Self::new(vec![F::ZERO]), p); + } + let mut q_coeffs = vec![F::ZERO; diff as usize + 1]; + + while diff >= 0 { + let s = p.leading_coefficient() * c.inverse().unwrap(); + q_coeffs[diff as usize] = s; + p -= rhs.pow_mult(s, diff as usize); + p.trim_zeros(); + diff = p.degree() as isize - rhs.degree() as isize; + } + q.coefficients = q_coeffs; + (q, p) + } +} + +impl Display for Polynomial { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + let mut first = true; + for (i, c) in self.coefficients.iter().enumerate() { + if !first { + write!(f, " + ")?; + } + first = false; + if i == 0 { + write!(f, "{}", c)?; + } else { + write!(f, "{}x^{}", c, i)?; + } + } + Ok(()) + } +} + +/// A polynomial in Lagrange basis +impl Polynomial, F> { + pub fn new(coefficients: Vec) -> Self { + // Check that the polynomial degree divides the field order so that there are roots of unity. + let n = coefficients.len(); + assert_eq!( + (F::ORDER - F::Storage::from(1_u32)) % F::Storage::from(n as u32), + F::Storage::from(0) + ); + let primitive_root = F::primitive_root_of_unity(F::Storage::from(n as u32)); + let nodes: Vec = (0..n).map(|i| primitive_root.pow(F::Storage::from(i as u32))).collect(); + + Self { coefficients, basis: Lagrange { nodes } } + } + + pub fn to_monomial(&self) -> Polynomial { + // This is the inverse of the conversion from monomial to Lagrange basis + // This uses something called the Vandermonde matrix which is defined as: + // + // / 1 | x_0 | x_0^2 | x_0^3 | ... | x_0^(N-1) \ + // | 1 | x_1 | x_1^2 | x_1^3 | ... | x_1^(N-1) | + // | 1 | x_2 | x_2^2 | x_2^3 | ... | x_2^(N-1) | + // v = | . | . | . | . | ... | . | + // | . | . | . | . | ... | . | + // | . | . | . | . | ... | . | + // \ 1 | x_N | x_N^2 | x_N^3 | ... | x_N^(N-1) / + // + // where x_i are the nodes of the Lagrange basis + // + // Then the monomial basis m is given V^T * l = m, where l is the Lagrange basis + // because we know the monomial basis we need to compute to monomial coefficients a_m = V^{-1} * + // a_l where a_l are the coefficients of the Lagrange basis + + // It also is the case that the the columns of the inverse matrix are the coefficients of the + // Lagrange polynomial basis TODO Finish this. + // let nodes = self.basis.nodes; + // let mut evaluations = [F::ZERO; N]; + + // // Evaluate the polynomial at N distinct points + // for i in 0..N { + // let x = F::primitive_root().exp_u64(i as u64); + // evaluations[i] = self.evaluate(x); + // } + + // Polynomial::::new(evaluations) + todo!("Finish this after we get the roots of unity from other PRs") + } + + /// Evaluate the polynomial at field element x + pub fn evaluate(&self, x: F) -> F { + let n = self.coefficients.len(); + let mut weights = vec![F::ZERO; n]; + for j in 0..n { + for m in 0..n { + if j != m { + weights[j] *= F::ONE.div(self.basis.nodes[j] - self.basis.nodes[m]); + } + } + } + // l(x) = \Sigmm_{m}(x-x_m) + let l = move |x: F| { + let mut val = F::ONE; + for i in 0..n { + val *= x - self.basis.nodes[i]; + } + val + }; + + // L(x) = l(x) * \Sigma_{j=0}^{k} (w_j / (x - x_j)) y_j + let mut val = F::ZERO; + for j in 0..n { + // check if we are dividing by zero + if self.basis.nodes[j] == x { + return self.coefficients[j]; + } + val += self.coefficients[j] * weights[j] / (x - self.basis.nodes[j]); + } + l(x) * val + } +} +mod test { + use super::*; + use crate::field::gf_101::GF101; + + fn deg_three_poly() -> Polynomial { + // Coefficients of the polynomial 1 + 2x + 3x^2 + 4x^3 + let a = GF101::new(1); + let b = GF101::new(2); + let c = GF101::new(3); + let d = GF101::new(4); + Polynomial::::new(vec![a, b, c, d]) + } + + #[test] + fn evaluation() { + // Should get: 1 + 2*(2) + 3*(2)^2 + 4*(2)^3 = 49 + let y = deg_three_poly().evaluate(GF101::new(2)); + let r = GF101::new(49); + assert_eq!(y, r); + } + + #[test] + fn evaluation_with_zero() { + // Coefficients of the polynomial 1 + 3x^2 + let a = GF101::new(1); + let b = GF101::new(0); + let c = GF101::new(3); + let polynomial = Polynomial::::new(vec![a, b, c]); + let y = polynomial.evaluate(GF101::new(0)); + + // Should get: 1 + 3(0)^2 = 1 + let r = GF101::new(1); + assert_eq!(y, r); + } + + #[test] + fn lagrange_evaluation() { + // Convert to Lagrange basis using roots of unity + let lagrange = deg_three_poly().to_lagrange(); + + // Should get: 1 + 2*(2) + 3*(2)^2 + 4*(2)^2= 49 + let r = lagrange.evaluate(GF101::new(2)); + assert_eq!(r, GF101::new(49)); + } + + #[test] + #[should_panic] + fn no_roots_of_unity() { + // Coefficients of the polynomial 1 + 2x + let a = GF101::new(1); + let b = GF101::new(2); + let c = GF101::new(3); + let polynomial = Polynomial::::new(vec![a, b, c]); + polynomial.to_lagrange(); + } + + #[test] + fn check_coefficients() { + assert_eq!(deg_three_poly().coefficients, [ + GF101::new(1), + GF101::new(2), + GF101::new(3), + GF101::new(4) + ]); + + assert_eq!(deg_three_poly().to_lagrange().coefficients, [ + GF101::new(10), + GF101::new(79), + GF101::new(99), + GF101::new(18) + ]); + } + + #[test] + fn degree() { + assert_eq!(deg_three_poly().degree(), 3); + } + + #[test] + fn leading_coefficient() { + assert_eq!(deg_three_poly().leading_coefficient(), GF101::new(4)); + } + + #[test] + fn pow_mult() { + let poly = deg_three_poly(); + let pow_mult = poly.pow_mult(GF101::new(5), 2); + assert_eq!(pow_mult.coefficients, [ + GF101::new(0), + GF101::new(0), + GF101::new(5), + GF101::new(10), + GF101::new(15), + GF101::new(20) + ]); + } + + #[test] + fn trim_zeros() { + let mut poly = deg_three_poly(); + poly.coefficients.push(GF101::ZERO); + assert_eq!(poly.coefficients, [ + GF101::new(1), + GF101::new(2), + GF101::new(3), + GF101::new(4), + GF101::ZERO + ]); + poly.trim_zeros(); + assert_eq!(poly.coefficients, [GF101::new(1), GF101::new(2), GF101::new(3), GF101::new(4)]); + } + + #[test] + fn trim_to_zero() { + let mut poly = Polynomial::::new(vec![GF101::ZERO, GF101::ZERO]); + assert_eq!(poly.coefficients, [GF101::ZERO]); + } +}