1use crate::poly::domain::vanishing_poly::VanishingPolynomial;
  2use ark_ff::{batch_inversion_and_mul, PrimeField};
  3use ark_std::vec::Vec;
  4/// Struct describing Lagrange interpolation for a multiplicative coset I,
  5/// with |I| a power of 2.
  6/// TODO: Pull in lagrange poly explanation from libiop
  7#[derive(Clone)]
  8pub struct LagrangeInterpolator<F: PrimeField> {
  9    pub(crate) domain_order: usize,
 10    pub(crate) all_domain_elems: Vec<F>,
 11    pub(crate) v_inv_elems: Vec<F>,
 12    pub(crate) domain_vp: VanishingPolynomial<F>,
 13    poly_evaluations: Vec<F>,
 14}
 15
 16impl<F: PrimeField> LagrangeInterpolator<F> {
 17    /// Returns a lagrange interpolator, given the domain specification.
 18    pub fn new(
 19        domain_offset: F,
 20        domain_generator: F,
 21        domain_dim: u64,
 22        poly_evaluations: Vec<F>,
 23    ) -> Self {
 24        let domain_order = 1 << domain_dim;
 25        assert_eq!(poly_evaluations.len(), domain_order);
 26        let mut cur_elem = domain_offset;
 27        let mut all_domain_elems = vec![domain_offset];
 28        let mut v_inv_elems: Vec<F> = Vec::new();
 29        // Cache all elements in the domain
 30        for _ in 1..domain_order {
 31            cur_elem *= domain_generator;
 32            all_domain_elems.push(cur_elem);
 33        }
 34        // By computing the following elements as constants,
 35        // we can further reduce the interpolation costs.
 36        //
 37        // m = order of the interpolation domain
 38        // v_inv[i] = prod_{j != i} h(g^i - g^j)
 39        // We use the following facts to compute this:
 40        // v_inv[0] = m*h^{m-1}
 41        // v_inv[i] = g^{-1} * v_inv[i-1]
 42        // TODO: Include proof of the above two points
 43        let g_inv = domain_generator.inverse().unwrap();
 44        let m = F::from((1 << domain_dim) as u128);
 45        let mut v_inv_i = m * domain_offset.pow([(domain_order - 1) as u64]);
 46        for _ in 0..domain_order {
 47            v_inv_elems.push(v_inv_i);
 48            v_inv_i *= g_inv;
 49        }
 50
 51        // TODO: Cache the intermediate terms with Z_H(x) evaluations.
 52        let vp = VanishingPolynomial::new(domain_offset, domain_dim);
 53
 54        let lagrange_interpolation: LagrangeInterpolator<F> = LagrangeInterpolator {
 55            domain_order,
 56            all_domain_elems,
 57            v_inv_elems,
 58            domain_vp: vp,
 59            poly_evaluations,
 60        };
 61        lagrange_interpolation
 62    }
 63
 64    pub(crate) fn compute_lagrange_coefficients(&self, interpolation_point: F) -> Vec<F> {
 65        // Let t be the interpolation point, H be the multiplicative coset, with
 66        // elements of the form h*g^i. Compute each L_{i,H}(t) as Z_{H}(t) * v_i
 67        // / (t- h g^i) where:
 68        // - Z_{H}(t) = \prod_{j} (t-h*g^j) = (t^m-h^m), and
 69        // - v_{i} = 1 / \prod_{j \neq i} h(g^i-g^j).
 70        // Below we use the fact that v_{0} = 1/(m * h^(m-1)) and v_{i+1} = g * v_{i}.
 71        // We first compute the inverse of each coefficient, except for the Z_H(t) term.
 72        // We then batch invert the entire result, and multiply by Z_H(t).
 73        let mut inverted_lagrange_coeffs: Vec<F> = Vec::with_capacity(self.all_domain_elems.len());
 74        for i in 0..self.domain_order {
 75            let l = self.v_inv_elems[i];
 76            let r = self.all_domain_elems[i];
 77            inverted_lagrange_coeffs.push(l * (interpolation_point - r));
 78        }
 79        let vp_t = self.domain_vp.evaluate(&interpolation_point);
 80        let lagrange_coeffs = inverted_lagrange_coeffs.as_mut_slice();
 81        batch_inversion_and_mul::<F>(lagrange_coeffs, &vp_t);
 82        lagrange_coeffs.iter().cloned().collect()
 83    }
 84
 85    pub fn interpolate(&self, interpolation_point: F) -> F {
 86        let lagrange_coeffs = self.compute_lagrange_coefficients(interpolation_point);
 87        let mut interpolation = F::zero();
 88        for i in 0..self.domain_order {
 89            interpolation += lagrange_coeffs[i] * self.poly_evaluations[i];
 90        }
 91        interpolation
 92    }
 93}
 94
 95#[cfg(test)]
 96mod tests {
 97    use crate::{
 98        fields::{fp::FpVar, FieldVar},
 99        poly::{
100            domain::Radix2DomainVar,
101            evaluations::univariate::lagrange_interpolator::LagrangeInterpolator,
102        },
103        R1CSVar,
104    };
105    use ark_ff::{FftField, Field, One};
106    use ark_poly::{univariate::DensePolynomial, DenseUVPolynomial, Polynomial};
107    use ark_std::{test_rng, UniformRand};
108    use ark_test_curves::bls12_381::Fr;
109
110    #[test]
111    pub fn test_native_interpolate() {
112        let mut rng = test_rng();
113        let poly = DensePolynomial::rand(15, &mut rng);
114        let gen = Fr::get_root_of_unity(1 << 4).unwrap();
115        assert_eq!(gen.pow(&[1 << 4]), Fr::one());
116        let domain = Radix2DomainVar::new(
117            gen,
118            4, // 2^4 = 16
119            FpVar::constant(Fr::GENERATOR),
120        )
121        .unwrap();
122        // generate evaluations of `poly` on this domain
123        let mut coset_point = domain.offset().value().unwrap();
124        let mut oracle_evals = Vec::new();
125        for _ in 0..(1 << 4) {
126            oracle_evals.push(poly.evaluate(&coset_point));
127            coset_point *= gen;
128        }
129
130        let interpolator = LagrangeInterpolator::new(
131            domain.offset().value().unwrap(),
132            domain.gen,
133            domain.dim,
134            oracle_evals,
135        );
136
137        // the point to evaluate at
138        let interpolate_point = Fr::rand(&mut rng);
139
140        let expected = poly.evaluate(&interpolate_point);
141        let actual = interpolator.interpolate(interpolate_point);
142
143        assert_eq!(actual, expected)
144    }
145}