1use ark_ff::PrimeField;
  2use ark_poly::{univariate::DensePolynomial, DenseUVPolynomial};
  3
  4/// Computes the lagrange interpolated polynomial from the given points `p_i`
  5pub fn compute_lagrange_interpolated_poly<F: PrimeField>(p_i: &[F]) -> DensePolynomial<F> {
  6    // domain is 0..p_i.len(), to fit `interpolate_uni_poly` from hyperplonk
  7    let domain: Vec<usize> = (0..p_i.len()).collect();
  8
  9    // compute l(x), common to every basis polynomial
 10    let mut l_x = DensePolynomial::from_coefficients_vec(vec![F::ONE]);
 11    for x_m in domain.clone() {
 12        let prod_m = DensePolynomial::from_coefficients_vec(vec![-F::from(x_m as u64), F::ONE]);
 13        l_x = &l_x * &prod_m;
 14    }
 15
 16    // compute each w_j - barycentric weights
 17    let mut w_j_vector: Vec<F> = vec![];
 18    for x_j in domain.clone() {
 19        let mut w_j = F::ONE;
 20        for x_m in domain.clone() {
 21            if x_m != x_j {
 22                let prod = (F::from(x_j as u64) - F::from(x_m as u64))
 23                    .inverse()
 24                    .unwrap(); // an inverse always exists since x_j != x_m (!=0)
 25                               // hence, we call unwrap() here without checking the Option's content
 26                w_j *= prod;
 27            }
 28        }
 29        w_j_vector.push(w_j);
 30    }
 31
 32    // compute each polynomial within the sum L(x)
 33    let mut lagrange_poly = DensePolynomial::from_coefficients_vec(vec![F::ZERO]);
 34    for (j, w_j) in w_j_vector.iter().enumerate() {
 35        let x_j = domain[j];
 36        let y_j = p_i[j];
 37        // we multiply by l(x) here, otherwise the below division will not work - deg(0)/deg(d)
 38        let poly_numerator = &(&l_x * (*w_j)) * (y_j);
 39        let poly_denominator =
 40            DensePolynomial::from_coefficients_vec(vec![-F::from(x_j as u64), F::ONE]);
 41        let poly = &poly_numerator / &poly_denominator;
 42        lagrange_poly = &lagrange_poly + &poly;
 43    }
 44
 45    lagrange_poly
 46}
 47
 48#[cfg(test)]
 49mod tests {
 50
 51    use crate::utils::espresso::sum_check::verifier::interpolate_uni_poly;
 52    use crate::utils::lagrange_poly::compute_lagrange_interpolated_poly;
 53    use ark_pallas::Fr;
 54    use ark_poly::{univariate::DensePolynomial, DenseUVPolynomial, Polynomial};
 55    use ark_std::{vec::Vec, UniformRand};
 56    use espresso_subroutines::poly_iop::prelude::PolyIOPErrors;
 57
 58    #[test]
 59    fn test_compute_lagrange_interpolated_poly() {
 60        let mut prng = ark_std::test_rng();
 61        for degree in 1..30 {
 62            let poly = DensePolynomial::<Fr>::rand(degree, &mut prng);
 63            // range (which is exclusive) is from 0 to degree + 1, since we need degree + 1 evaluations
 64            let evals = (0..(degree + 1))
 65                .map(|i| poly.evaluate(&Fr::from(i as u64)))
 66                .collect::<Vec<Fr>>();
 67            let lagrange_poly = compute_lagrange_interpolated_poly(&evals);
 68            for _ in 0..10 {
 69                let query = Fr::rand(&mut prng);
 70                let lagrange_eval = lagrange_poly.evaluate(&query);
 71                let eval = poly.evaluate(&query);
 72                assert_eq!(eval, lagrange_eval);
 73                assert_eq!(lagrange_poly.degree(), poly.degree());
 74            }
 75        }
 76    }
 77
 78    #[test]
 79    fn test_interpolation() -> Result<(), PolyIOPErrors> {
 80        let mut prng = ark_std::test_rng();
 81
 82        // test a polynomial with 20 known points, i.e., with degree 19
 83        let poly = DensePolynomial::<Fr>::rand(20 - 1, &mut prng);
 84        let evals = (0..20)
 85            .map(|i| poly.evaluate(&Fr::from(i)))
 86            .collect::<Vec<Fr>>();
 87        let query = Fr::rand(&mut prng);
 88
 89        assert_eq!(poly.evaluate(&query), interpolate_uni_poly(&evals, query)?);
 90        assert_eq!(
 91            compute_lagrange_interpolated_poly(&evals).evaluate(&query),
 92            interpolate_uni_poly(&evals, query)?
 93        );
 94
 95        // test a polynomial with 33 known points, i.e., with degree 32
 96        let poly = DensePolynomial::<Fr>::rand(33 - 1, &mut prng);
 97        let evals = (0..33)
 98            .map(|i| poly.evaluate(&Fr::from(i)))
 99            .collect::<Vec<Fr>>();
100        let query = Fr::rand(&mut prng);
101
102        assert_eq!(poly.evaluate(&query), interpolate_uni_poly(&evals, query)?);
103        assert_eq!(
104            compute_lagrange_interpolated_poly(&evals).evaluate(&query),
105            interpolate_uni_poly(&evals, query)?
106        );
107
108        // test a polynomial with 64 known points, i.e., with degree 63
109        let poly = DensePolynomial::<Fr>::rand(64 - 1, &mut prng);
110        let evals = (0..64)
111            .map(|i| poly.evaluate(&Fr::from(i)))
112            .collect::<Vec<Fr>>();
113        let query = Fr::rand(&mut prng);
114
115        assert_eq!(poly.evaluate(&query), interpolate_uni_poly(&evals, query)?);
116        assert_eq!(
117            compute_lagrange_interpolated_poly(&evals).evaluate(&query),
118            interpolate_uni_poly(&evals, query)?
119        );
120
121        Ok(())
122    }
123}