ark-r1cs-std/src/poly/evaluations/univariate/mod.rs
src/poly/evaluations/univariate/mod.rsRaw
1pub mod lagrange_interpolator;
2
3use crate::alloc::AllocVar;
4use crate::eq::EqGadget;
5use crate::fields::fp::FpVar;
6use crate::fields::FieldVar;
7use crate::poly::domain::Radix2DomainVar;
8use crate::poly::evaluations::univariate::lagrange_interpolator::LagrangeInterpolator;
9use crate::R1CSVar;
10use ark_ff::{batch_inversion, PrimeField};
11use ark_relations::r1cs::SynthesisError;
12use ark_std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Sub, SubAssign};
13use ark_std::vec::Vec;
14
15#[derive(Clone)]
16/// Stores a UV polynomial in evaluation form.
17pub struct EvaluationsVar<F: PrimeField> {
18 /// Evaluations of univariate polynomial over domain
19 pub evals: Vec<FpVar<F>>,
20 /// Optional Lagrange Interpolator. Useful for lagrange interpolation.
21 pub lagrange_interpolator: Option<LagrangeInterpolator<F>>,
22 domain: Radix2DomainVar<F>,
23 /// Contains all domain elements of `domain.base_domain`.
24 ///
25 /// This is a cache for lagrange interpolation when offset is non-constant. Will be `None` if offset is constant
26 /// or `interpolate` is set to `false`.
27 subgroup_points: Option<Vec<F>>,
28}
29
30impl<F: PrimeField> EvaluationsVar<F> {
31 /// Construct `Self` from evaluations and a domain.
32 /// `interpolate` indicates if user wants to interpolate this polynomial
33 /// using lagrange interpolation.
34 pub fn from_vec_and_domain(
35 evaluations: Vec<FpVar<F>>,
36 domain: Radix2DomainVar<F>,
37 interpolate: bool,
38 ) -> Self {
39 assert_eq!(
40 evaluations.len(),
41 1 << domain.dim,
42 "evaluations and domain has different dimensions"
43 );
44
45 let mut ev = Self {
46 evals: evaluations,
47 lagrange_interpolator: None,
48 domain,
49 subgroup_points: None,
50 };
51 if interpolate {
52 ev.generate_interpolation_cache();
53 }
54 ev
55 }
56
57 /// Precompute necessary calculation for lagrange interpolation and mark it ready to interpolate
58 pub fn generate_interpolation_cache(&mut self) {
59 if self.domain.offset().is_constant() {
60 let poly_evaluations_val: Vec<_> =
61 self.evals.iter().map(|v| v.value().unwrap()).collect();
62 let domain = &self.domain;
63 let lagrange_interpolator = if let &FpVar::Constant(x) = domain.offset() {
64 LagrangeInterpolator::new(x, domain.gen, domain.dim, poly_evaluations_val)
65 } else {
66 panic!("Domain offset needs to be constant.")
67 };
68 self.lagrange_interpolator = Some(lagrange_interpolator)
69 } else {
70 // calculate all elements of base subgroup so that in later part we don't need to calculate the exponents again
71 let mut subgroup_points = Vec::with_capacity(self.domain.size() as usize);
72 subgroup_points.push(F::one());
73 for i in 1..self.domain.size() as usize {
74 subgroup_points.push(subgroup_points[i - 1] * self.domain.gen)
75 }
76 self.subgroup_points = Some(subgroup_points)
77 }
78 }
79
80 /// Compute lagrange coefficients for each evaluation, given `interpolation_point`.
81 /// Only valid if the domain offset is constant.
82 fn compute_lagrange_coefficients(
83 &self,
84 interpolation_point: &FpVar<F>,
85 ) -> Result<Vec<FpVar<F>>, SynthesisError> {
86 // ref: https://github.com/alexchmit/perfect-constraints/blob/79692f2652a95a57f2c7187f5b5276345e680230/fractal/src/algebra/lagrange_interpolation.rs#L159
87 let cs = interpolation_point.cs();
88 let t = interpolation_point;
89 let lagrange_interpolator = self
90 .lagrange_interpolator
91 .as_ref()
92 .expect("lagrange interpolator has not been initialized. \
93 Call `self.generate_interpolation_cache` first or set `interpolate` to true in constructor. ");
94 let lagrange_coeffs =
95 lagrange_interpolator.compute_lagrange_coefficients(t.value().unwrap());
96 let mut lagrange_coeffs_fg = Vec::new();
97 // Now we convert these lagrange coefficients to gadgets, and then constrain them.
98 // The i-th lagrange coefficients constraint is:
99 // (v_inv[i] * t - v_inv[i] * domain_elem[i]) * (coeff) = 1/Z_I(t)
100 let vp_t = lagrange_interpolator.domain_vp.evaluate_constraints(t)?;
101 // let inv_vp_t = vp_t.inverse()?;
102 for i in 0..lagrange_interpolator.domain_order {
103 let constant: F =
104 (-lagrange_interpolator.all_domain_elems[i]) * lagrange_interpolator.v_inv_elems[i];
105 let mut a_element: FpVar<F> =
106 t * &FpVar::constant(lagrange_interpolator.v_inv_elems[i]);
107 a_element += FpVar::constant(constant);
108
109 let lag_coeff: FpVar<F> =
110 FpVar::new_witness(ns!(cs, "generate lagrange coefficient"), || {
111 Ok(lagrange_coeffs[i])
112 })?;
113 // Enforce the actual constraint (A_element) * (lagrange_coeff) = 1/Z_I(t)
114 assert_eq!(
115 (lagrange_interpolator.v_inv_elems[i] * t.value().unwrap()
116 - lagrange_interpolator.v_inv_elems[i]
117 * lagrange_interpolator.all_domain_elems[i])
118 * lagrange_coeffs[i],
119 vp_t.value().unwrap()
120 );
121 a_element.mul_equals(&lag_coeff, &vp_t)?;
122 lagrange_coeffs_fg.push(lag_coeff);
123 }
124 Ok(lagrange_coeffs_fg)
125 }
126
127 /// Returns constraints for Interpolating and evaluating at `interpolation_point`
128 pub fn interpolate_and_evaluate(
129 &self,
130 interpolation_point: &FpVar<F>,
131 ) -> Result<FpVar<F>, SynthesisError> {
132 // specialize: if domain offset is constant, we can optimize to have fewer constraints
133 if self.domain.offset().is_constant() {
134 self.lagrange_interpolate_with_constant_offset(interpolation_point)
135 } else {
136 // if domain offset is not constant, then we use standard lagrange interpolation code
137 self.lagrange_interpolate_with_non_constant_offset(interpolation_point)
138 }
139 }
140
141 fn lagrange_interpolate_with_constant_offset(
142 &self,
143 interpolation_point: &FpVar<F>,
144 ) -> Result<FpVar<F>, SynthesisError> {
145 let lagrange_interpolator = self
146 .lagrange_interpolator
147 .as_ref()
148 .expect("lagrange interpolator has not been initialized. ");
149 let lagrange_coeffs = self.compute_lagrange_coefficients(interpolation_point)?;
150 let mut interpolation: FpVar<F> = FpVar::zero();
151 for i in 0..lagrange_interpolator.domain_order {
152 let intermediate = &lagrange_coeffs[i] * &self.evals[i];
153 interpolation += &intermediate
154 }
155
156 Ok(interpolation)
157 }
158
159 /// Generate interpolation constraints. We assume at compile time we know the base coset (i.e. `gen`) but not know `offset`.
160 fn lagrange_interpolate_with_non_constant_offset(
161 &self,
162 interpolation_point: &FpVar<F>, // = alpha in the following code
163 ) -> Result<FpVar<F>, SynthesisError> {
164 // first, make sure `subgroup_points` is made
165 let subgroup_points = self.subgroup_points.as_ref()
166 .expect("lagrange interpolator has not been initialized. \
167 Call `self.generate_interpolation_cache` first or set `interpolate` to true in constructor. ");
168 // Let denote interpolation_point as alpha.
169 // Lagrange polynomial for coset element `a` is
170 // \frac{1}{size * offset ^ size} * \frac{alpha^size - offset^size}{alpha * a^{-1} - 1}
171 // Notice that a = (offset * a') where a' is the corresponding element of base coset
172
173 // let `lhs` be \frac{alpha^size - offset^size}{size * offset ^ size}. This part is shared by all lagrange polynomials
174 let coset_offset_to_size = self
175 .domain
176 .offset()
177 .pow_by_constant(&[self.domain.size()])?; // offset^size
178 let alpha_to_size = interpolation_point.pow_by_constant(&[self.domain.size()])?;
179 let lhs_numerator = &alpha_to_size - &coset_offset_to_size;
180 // This enforces that `alpha` is not in the coset.
181 // This also means that the denominator is
182 lhs_numerator.enforce_not_equal(&FpVar::zero())?;
183
184 // `domain.offset()` is non-zero by construction, so `coset_offset_to_size` is also non-zero, which means `lhs_denominator` is non-zero
185 let lhs_denominator = &coset_offset_to_size * FpVar::constant(F::from(self.domain.size()));
186
187 // unchecked is okay because the denominator is non-zero.
188 let lhs = lhs_numerator.mul_by_inverse_unchecked(&lhs_denominator)?;
189
190 // `rhs` for coset element `a` is \frac{1}{alpha * a^{-1} - 1} = \frac{1}{alpha * offset^{-1} * a'^{-1} - 1}
191 // domain.offset() is non-zero by construction.
192 let alpha_coset_offset_inv =
193 interpolation_point.mul_by_inverse_unchecked(&self.domain.offset())?;
194
195 // `res` stores the sum of all lagrange polynomials evaluated at alpha
196 let mut res = FpVar::<F>::zero();
197
198 let domain_size = self.domain.size() as usize;
199 for i in 0..domain_size {
200 // a'^{-1} where a is the base coset element
201 let subgroup_point_inv = subgroup_points[(domain_size - i) % domain_size];
202 debug_assert_eq!(subgroup_points[i] * subgroup_point_inv, F::one());
203 // alpha * offset^{-1} * a'^{-1} - 1
204 let lag_denom = &alpha_coset_offset_inv * subgroup_point_inv - F::one();
205 // lag_denom cannot be zero, so we use `unchecked`.
206 //
207 // Proof: lag_denom is zero if and only if alpha * (coset_offset * subgroup_point)^{-1} == 1.
208 // This can happen only if `alpha` is itself in the coset.
209 //
210 // Earlier we asserted that `lhs_numerator` is not zero.
211 // Since `lhs_numerator` is just the vanishing polynomial for the coset evaluated at `alpha`,
212 // and since this is non-zero, `alpha` is not in the coset.
213 let lag_coeff = lhs.mul_by_inverse_unchecked(&lag_denom)?;
214
215 let lag_interpoland = &self.evals[i] * lag_coeff;
216 res += lag_interpoland
217 }
218
219 Ok(res)
220 }
221}
222
223impl<'a, 'b, F: PrimeField> Add<&'a EvaluationsVar<F>> for &'b EvaluationsVar<F> {
224 type Output = EvaluationsVar<F>;
225
226 fn add(self, rhs: &'a EvaluationsVar<F>) -> Self::Output {
227 let mut result = self.clone();
228 result += rhs;
229 result
230 }
231}
232
233impl<'a, F: PrimeField> AddAssign<&'a EvaluationsVar<F>> for EvaluationsVar<F> {
234 /// Performs the `+=` operations, assuming `domain.offset` is equal.
235 fn add_assign(&mut self, other: &'a EvaluationsVar<F>) {
236 // offset might be unknown at compile time, so we assume offset is equal
237 assert!(
238 self.domain.gen == other.domain.gen && self.domain.dim == other.domain.dim,
239 "domains are unequal"
240 );
241
242 self.lagrange_interpolator = None;
243 self.evals
244 .iter_mut()
245 .zip(&other.evals)
246 .for_each(|(a, b)| *a = &*a + b)
247 }
248}
249
250impl<'a, 'b, F: PrimeField> Sub<&'a EvaluationsVar<F>> for &'b EvaluationsVar<F> {
251 type Output = EvaluationsVar<F>;
252
253 fn sub(self, rhs: &'a EvaluationsVar<F>) -> Self::Output {
254 let mut result = self.clone();
255 result -= rhs;
256 result
257 }
258}
259
260impl<'a, F: PrimeField> SubAssign<&'a EvaluationsVar<F>> for EvaluationsVar<F> {
261 /// Performs the `-=` operations, assuming `domain.offset` is equal.
262 fn sub_assign(&mut self, other: &'a EvaluationsVar<F>) {
263 // offset might be unknown at compile time, so we assume offset is equal
264 assert!(
265 self.domain.gen == other.domain.gen && self.domain.dim == other.domain.dim,
266 "domains are unequal"
267 );
268
269 self.lagrange_interpolator = None;
270 self.evals
271 .iter_mut()
272 .zip(&other.evals)
273 .for_each(|(a, b)| *a = &*a - b)
274 }
275}
276
277impl<'a, 'b, F: PrimeField> Mul<&'a EvaluationsVar<F>> for &'b EvaluationsVar<F> {
278 type Output = EvaluationsVar<F>;
279
280 /// Performs the `*` operations, assuming `domain.offset` is equal.
281 fn mul(self, rhs: &'a EvaluationsVar<F>) -> Self::Output {
282 let mut result = self.clone();
283 result *= rhs;
284 result
285 }
286}
287
288impl<'a, F: PrimeField> MulAssign<&'a EvaluationsVar<F>> for EvaluationsVar<F> {
289 /// Performs the `*=` operations, assuming `domain.offset` is equal.
290 fn mul_assign(&mut self, other: &'a EvaluationsVar<F>) {
291 // offset might be unknown at compile time, so we assume offset is equal
292 assert!(
293 self.domain.gen == other.domain.gen && self.domain.dim == other.domain.dim,
294 "domains are unequal"
295 );
296
297 self.lagrange_interpolator = None;
298 self.evals
299 .iter_mut()
300 .zip(&other.evals)
301 .for_each(|(a, b)| *a = &*a * b)
302 }
303}
304
305impl<'a, 'b, F: PrimeField> Div<&'a EvaluationsVar<F>> for &'b EvaluationsVar<F> {
306 type Output = EvaluationsVar<F>;
307
308 fn div(self, rhs: &'a EvaluationsVar<F>) -> Self::Output {
309 let mut result = self.clone();
310 result /= rhs;
311 result
312 }
313}
314
315impl<'a, F: PrimeField> DivAssign<&'a EvaluationsVar<F>> for EvaluationsVar<F> {
316 /// Performs the `/=` operations, assuming `domain.offset` is equal.
317 fn div_assign(&mut self, other: &'a EvaluationsVar<F>) {
318 // offset might be unknown at compile time, so we assume offset is equal
319 assert!(
320 self.domain.gen == other.domain.gen && self.domain.dim == other.domain.dim,
321 "domains are unequal"
322 );
323 let cs = self.evals[0].cs();
324 // the prover can generate result = (1 / other) * self offline
325 let mut result_val: Vec<_> = other.evals.iter().map(|x| x.value().unwrap()).collect();
326 batch_inversion(&mut result_val);
327 result_val
328 .iter_mut()
329 .zip(&self.evals)
330 .for_each(|(a, self_var)| *a *= self_var.value().unwrap());
331 let result_var: Vec<_> = result_val
332 .iter()
333 .map(|x| FpVar::new_witness(ns!(cs, "div result"), || Ok(*x)).unwrap())
334 .collect();
335 // enforce constraint
336 for i in 0..result_var.len() {
337 result_var[i]
338 .mul_equals(&other.evals[i], &self.evals[i])
339 .unwrap();
340 }
341
342 self.lagrange_interpolator = None;
343 self.evals = result_var
344 }
345}
346
347#[cfg(test)]
348mod tests {
349 use crate::alloc::AllocVar;
350 use crate::fields::fp::FpVar;
351 use crate::fields::FieldVar;
352 use crate::poly::domain::Radix2DomainVar;
353 use crate::poly::evaluations::univariate::EvaluationsVar;
354 use crate::R1CSVar;
355 use ark_ff::{FftField, Field, One, UniformRand};
356 use ark_poly::polynomial::univariate::DensePolynomial;
357 use ark_poly::{Polynomial, UVPolynomial};
358 use ark_relations::r1cs::ConstraintSystem;
359 use ark_std::test_rng;
360 use ark_test_curves::bls12_381::Fr;
361
362 #[test]
363 fn test_interpolate_constant_offset() {
364 let mut rng = test_rng();
365 let poly = DensePolynomial::rand(15, &mut rng);
366 let gen = Fr::get_root_of_unity(1 << 4).unwrap();
367 assert_eq!(gen.pow(&[1 << 4]), Fr::one());
368 let domain = Radix2DomainVar::new(
369 gen,
370 4, // 2^4 = 16
371 FpVar::constant(Fr::rand(&mut rng)),
372 )
373 .unwrap();
374 let mut coset_point = domain.offset().value().unwrap();
375 let mut oracle_evals = Vec::new();
376 for _ in 0..(1 << 4) {
377 oracle_evals.push(poly.evaluate(&coset_point));
378 coset_point *= gen;
379 }
380 let cs = ConstraintSystem::new_ref();
381 let evaluations_fp: Vec<_> = oracle_evals
382 .iter()
383 .map(|x| FpVar::new_input(ns!(cs, "evaluations"), || Ok(x)).unwrap())
384 .collect();
385 let evaluations_var = EvaluationsVar::from_vec_and_domain(evaluations_fp, domain, true);
386
387 let interpolate_point = Fr::rand(&mut rng);
388 let interpolate_point_fp =
389 FpVar::new_input(ns!(cs, "interpolate point"), || Ok(interpolate_point)).unwrap();
390
391 let expected = poly.evaluate(&interpolate_point);
392
393 let actual = evaluations_var
394 .interpolate_and_evaluate(&interpolate_point_fp)
395 .unwrap()
396 .value()
397 .unwrap();
398
399 assert_eq!(actual, expected);
400 assert!(cs.is_satisfied().unwrap());
401 println!("number of constraints: {}", cs.num_constraints())
402 }
403
404 #[test]
405 fn test_interpolate_non_constant_offset() {
406 let mut rng = test_rng();
407 let poly = DensePolynomial::rand(15, &mut rng);
408 let gen = Fr::get_root_of_unity(1 << 4).unwrap();
409 assert_eq!(gen.pow(&[1 << 4]), Fr::one());
410 let cs = ConstraintSystem::new_ref();
411 let domain = Radix2DomainVar::new(
412 gen,
413 4, // 2^4 = 16
414 FpVar::new_witness(ns!(cs, "offset"), || Ok(Fr::rand(&mut rng))).unwrap(),
415 )
416 .unwrap();
417 let mut coset_point = domain.offset().value().unwrap();
418 let mut oracle_evals = Vec::new();
419 for _ in 0..(1 << 4) {
420 oracle_evals.push(poly.evaluate(&coset_point));
421 coset_point *= gen;
422 }
423
424 let evaluations_fp: Vec<_> = oracle_evals
425 .iter()
426 .map(|x| FpVar::new_input(ns!(cs, "evaluations"), || Ok(x)).unwrap())
427 .collect();
428 let evaluations_var = EvaluationsVar::from_vec_and_domain(evaluations_fp, domain, true);
429
430 let interpolate_point = Fr::rand(&mut rng);
431 let interpolate_point_fp =
432 FpVar::new_input(ns!(cs, "interpolate point"), || Ok(interpolate_point)).unwrap();
433
434 let expected = poly.evaluate(&interpolate_point);
435
436 let actual = evaluations_var
437 .interpolate_and_evaluate(&interpolate_point_fp)
438 .unwrap()
439 .value()
440 .unwrap();
441
442 assert_eq!(actual, expected);
443 assert!(cs.is_satisfied().unwrap());
444 println!("number of constraints: {}", cs.num_constraints())
445 }
446
447 #[test]
448 fn test_division() {
449 let mut rng = test_rng();
450 let gen = Fr::get_root_of_unity(1 << 4).unwrap();
451 assert_eq!(gen.pow(&[1 << 4]), Fr::one());
452 let domain = Radix2DomainVar::new(
453 gen,
454 4, // 2^4 = 16
455 FpVar::constant(Fr::GENERATOR),
456 )
457 .unwrap();
458 let cs = ConstraintSystem::new_ref();
459
460 let ev_a = EvaluationsVar::from_vec_and_domain(
461 (0..16)
462 .map(|_| FpVar::new_input(ns!(cs, "poly_a"), || Ok(Fr::rand(&mut rng))).unwrap())
463 .collect(),
464 domain.clone(),
465 false,
466 );
467 let ev_b = EvaluationsVar::from_vec_and_domain(
468 (0..16)
469 .map(|_| FpVar::new_input(ns!(cs, "poly_a"), || Ok(Fr::rand(&mut rng))).unwrap())
470 .collect(),
471 domain.clone(),
472 false,
473 );
474
475 let a_div_b = (&ev_a) / (&ev_b);
476 assert!(cs.is_satisfied().unwrap());
477 let b_div_a = (&ev_b) / (&ev_a);
478
479 let one = &a_div_b * &b_div_a;
480 for ev in one.evals.iter() {
481 assert!(Fr::is_one(&ev.value().unwrap()))
482 }
483
484 assert!(cs.is_satisfied().unwrap());
485 }
486}