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