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