ark_r1cs_std/poly/evaluations/univariate/
lagrange_interpolator.rs

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}
ark_r1cs_std/poly/evaluations/univariate/lagrange_interpolator.rs

ark_r1cs_std/poly/evaluations/univariate/
lagrange_interpolator.rs
