hyperplonk/pcs/src/util.rs
pcs/src/util.rsRaw
1//! Useful utilities for KZG PCS
2
3use crate::PCSErrors;
4use ark_ff::PrimeField;
5use ark_poly::{
6 univariate::DensePolynomial, DenseMultilinearExtension, EvaluationDomain, Evaluations,
7 MultilinearExtension, Polynomial, Radix2EvaluationDomain,
8};
9use ark_std::{end_timer, log2, start_timer};
10use poly_iop::bit_decompose;
11
12/// Compute W \circ l.
13///
14/// Given an MLE W, and a list of univariate polynomials l, generate the
15/// univariate polynomial that composes W with l.
16///
17/// Returns an error if l's length does not matches number of variables in W.
18pub(crate) fn compute_w_circ_l<F: PrimeField>(
19 w: &DenseMultilinearExtension<F>,
20 l: &[DensePolynomial<F>],
21) -> Result<DensePolynomial<F>, PCSErrors> {
22 let timer = start_timer!(|| "compute W \\circ l");
23
24 if w.num_vars != l.len() {
25 return Err(PCSErrors::InvalidParameters(format!(
26 "l's length ({}) does not match num_variables ({})",
27 l.len(),
28 w.num_vars(),
29 )));
30 }
31
32 let mut res_eval: Vec<F> = vec![];
33
34 // TODO: consider to pass this in from caller
35 // uni_degree is (product of each prefix's) + (2 * MLEs)
36 // = (l.len() - (num_vars - log(l.len())) + 2) * l[0].degree
37 let uni_degree = (l.len() - w.num_vars + log2(l.len()) as usize + 2) * l[0].degree();
38
39 let domain = match Radix2EvaluationDomain::<F>::new(uni_degree) {
40 Some(p) => p,
41 None => {
42 return Err(PCSErrors::InvalidParameters(
43 "failed to build radix 2 domain".to_string(),
44 ))
45 },
46 };
47 for point in domain.elements() {
48 // we reverse the order here because the coefficient vec are stored in
49 // bit-reversed order
50 let l_eval: Vec<F> = l.iter().rev().map(|x| x.evaluate(&point)).collect();
51 res_eval.push(w.evaluate(l_eval.as_ref()).unwrap())
52 }
53 let evaluation = Evaluations::from_vec_and_domain(res_eval, domain);
54 let res = evaluation.interpolate();
55
56 end_timer!(timer);
57 Ok(res)
58}
59
60/// Return the number of variables that one need for an MLE to
61/// batch the list of MLEs
62#[inline]
63pub(crate) fn get_batched_nv(num_var: usize, polynomials_len: usize) -> usize {
64 num_var + log2(polynomials_len) as usize
65}
66
67/// merge a set of polynomials. Returns an error if the
68/// polynomials do not share a same number of nvs.
69pub(crate) fn merge_polynomials<F: PrimeField>(
70 polynomials: &[impl MultilinearExtension<F>],
71) -> Result<DenseMultilinearExtension<F>, PCSErrors> {
72 let nv = polynomials[0].num_vars();
73 for poly in polynomials.iter() {
74 if nv != poly.num_vars() {
75 return Err(PCSErrors::InvalidParameters(
76 "num_vars do not match for polynomials".to_string(),
77 ));
78 }
79 }
80
81 let merged_nv = get_batched_nv(nv, polynomials.len());
82 let mut scalars = vec![];
83 for poly in polynomials.iter() {
84 scalars.extend_from_slice(poly.to_evaluations().as_slice());
85 }
86 scalars.extend_from_slice(vec![F::zero(); (1 << merged_nv) - scalars.len()].as_ref());
87 Ok(DenseMultilinearExtension::from_evaluations_vec(
88 merged_nv, scalars,
89 ))
90}
91
92/// Given a list of points, build `l(points)` which is a list of univariate
93/// polynomials that goes through the points
94pub(crate) fn build_l<F: PrimeField>(
95 num_var: usize,
96 points: &[&[F]],
97) -> Result<Vec<DensePolynomial<F>>, PCSErrors> {
98 let prefix_len = log2(points.len()) as usize;
99
100 let uni_degree = points.len();
101 let small_domain = match Radix2EvaluationDomain::<F>::new(uni_degree) {
102 Some(p) => p,
103 None => {
104 return Err(PCSErrors::InvalidParameters(
105 "failed to build radix 2 domain".to_string(),
106 ))
107 },
108 };
109
110 // The following code print out the roots for testing
111 // println!("domain root0: {}", small_domain.element(0));
112 // println!("domain root1: {}", small_domain.element(1));
113 // println!("domain root2: {}", small_domain.element(2));
114 // println!("domain root3: {}", small_domain.element(3));
115
116 let mut uni_polys = Vec::new();
117
118 // 1.1 build the indexes and the univariate polys that go through the indexes
119 let indexes: Vec<Vec<bool>> = (0..points.len())
120 .map(|x| bit_decompose(x as u64, prefix_len))
121 .collect();
122 for i in 0..prefix_len {
123 let eval: Vec<F> = indexes
124 .iter()
125 .map(|x| F::from(x[prefix_len - i - 1]))
126 .collect();
127
128 uni_polys.push(Evaluations::from_vec_and_domain(eval, small_domain).interpolate());
129 }
130
131 // 1.2 build the actual univariate polys that go through the points
132 for i in 0..num_var {
133 let mut eval: Vec<F> = points.iter().map(|x| x[i]).collect();
134 eval.extend_from_slice(vec![F::zero(); small_domain.size as usize - eval.len()].as_slice());
135 uni_polys.push(Evaluations::from_vec_and_domain(eval, small_domain).interpolate())
136 }
137
138 Ok(uni_polys)
139}
140
141#[cfg(test)]
142mod test {
143 use super::*;
144 use ark_bls12_381::Fr;
145 use ark_ff::field_new;
146 use ark_poly::UVPolynomial;
147 use ark_std::{One, Zero};
148
149 #[test]
150 fn test_w_circ_l() -> Result<(), PCSErrors> {
151 test_w_circ_l_helper::<Fr>()
152 }
153
154 fn test_w_circ_l_helper<F: PrimeField>() -> Result<(), PCSErrors> {
155 {
156 // Example from page 53:
157 // W = 3x1x2 + 2x2 whose evaluations are
158 // 0, 0 |-> 0
159 // 0, 1 |-> 2
160 // 1, 0 |-> 0
161 // 1, 1 |-> 5
162 let w_eval = vec![F::zero(), F::from(2u64), F::zero(), F::from(5u64)];
163 let w = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
164
165 // l0 = t + 2
166 // l1 = -2t + 4
167 let l0 = DensePolynomial::from_coefficients_vec(vec![F::from(2u64), F::one()]);
168 let l1 = DensePolynomial::from_coefficients_vec(vec![F::from(4u64), -F::from(2u64)]);
169
170 // res = -6t^2 - 4t + 32
171 let res = compute_w_circ_l(&w, [l0, l1].as_ref())?;
172 let res_rec = DensePolynomial::from_coefficients_vec(vec![
173 F::from(32u64),
174 -F::from(4u64),
175 -F::from(6u64),
176 ]);
177 assert_eq!(res, res_rec);
178 }
179 {
180 // A random example
181 // W = x1x2x3 - 2x1x2 + 3x2x3 - 4x1x3 + 5x1 - 6x2 + 7x3
182 // 0, 0, 0 |-> 0
183 // 0, 0, 1 |-> 7
184 // 0, 1, 0 |-> -6
185 // 0, 1, 1 |-> 4
186 // 1, 0, 0 |-> 5
187 // 1, 0, 1 |-> 8
188 // 1, 1, 0 |-> -3
189 // 1, 1, 1 |-> 4
190 let w_eval = vec![
191 F::zero(),
192 F::from(7u64),
193 -F::from(6u64),
194 F::from(4u64),
195 F::from(5u64),
196 F::from(8u64),
197 -F::from(3u64),
198 F::from(4u64),
199 ];
200 let w = DenseMultilinearExtension::from_evaluations_vec(3, w_eval);
201
202 // l0 = t + 2
203 // l1 = 3t - 4
204 // l2 = -5t + 6
205 let l0 = DensePolynomial::from_coefficients_vec(vec![F::from(2u64), F::one()]);
206 let l1 = DensePolynomial::from_coefficients_vec(vec![-F::from(4u64), F::from(3u64)]);
207 let l2 = DensePolynomial::from_coefficients_vec(vec![F::from(6u64), -F::from(5u64)]);
208 let res = compute_w_circ_l(&w, [l0, l1, l2].as_ref())?;
209
210 // res = -15t^3 - 23t^2 + 130t - 76
211 let res_rec = DensePolynomial::from_coefficients_vec(vec![
212 -F::from(76u64),
213 F::from(130u64),
214 -F::from(23u64),
215 -F::from(15u64),
216 ]);
217
218 assert_eq!(res, res_rec);
219 }
220 Ok(())
221 }
222
223 #[test]
224 fn test_merge_poly() -> Result<(), PCSErrors> {
225 test_merge_poly_helper::<Fr>()
226 }
227 fn test_merge_poly_helper<F: PrimeField>() -> Result<(), PCSErrors> {
228 // Example from page 53:
229 // W1 = 3x1x2 + 2x2 whose evaluations are
230 // 0, 0 |-> 0
231 // 0, 1 |-> 2
232 // 1, 0 |-> 0
233 // 1, 1 |-> 5
234 let w_eval = vec![F::zero(), F::from(2u64), F::zero(), F::from(5u64)];
235 let w1 = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
236
237 // W2 = x1x2 + x1 whose evaluations are
238 // 0, 0 |-> 0
239 // 0, 1 |-> 0
240 // 1, 0 |-> 1
241 // 1, 1 |-> 2
242 let w_eval = vec![F::zero(), F::zero(), F::from(1u64), F::from(2u64)];
243 let w2 = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
244
245 // W3 = x1 + x2 whose evaluations are
246 // 0, 0 |-> 0
247 // 0, 1 |-> 1
248 // 1, 0 |-> 1
249 // 1, 1 |-> 2
250 let w_eval = vec![F::zero(), F::one(), F::from(1u64), F::from(2u64)];
251 let w3 = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
252
253 {
254 // W = (3x1x2 + 2x2)(1-x0) + (x1x2 + x1)x0
255 // = -2x0x1x2 + x0x1 - 2x0x2 + 3x1x2 + 2x2
256 // with evaluation map
257 //
258 // x0 x1 x2
259 // 0, 0, 0 |-> 0
260 // 0, 0, 1 |-> 2
261 // 0, 1, 0 |-> 0
262 // 0, 1, 1 |-> 5
263 // 1, 0, 0 |-> 0
264 // 1, 0, 1 |-> 0
265 // 1, 1, 0 |-> 1
266 // 1, 1, 1 |-> 2
267 //
268 let w = merge_polynomials(&[w1.clone(), w2.clone()])?;
269 // w is [0,2,0,5,0,0,1,2]
270 let w_eval = vec![
271 F::zero(),
272 F::from(2u64),
273 F::zero(),
274 F::from(5u64),
275 F::zero(),
276 F::zero(),
277 F::from(1u64),
278 F::from(2u64),
279 ];
280 let w_rec = DenseMultilinearExtension::from_evaluations_vec(3, w_eval);
281
282 assert_eq!(w, w_rec);
283 }
284
285 {
286 // W = (3x1x2 + 2x2) * (1-y1) * (1-y2)
287 // + (x1x2 + x1) * (1-y1) * y2
288 // + (x1 + x2) * y1 * (1-y2)
289 //
290 // with evaluation map
291 //
292 // y1 y2 x1 x2
293 // 0, 0, 0, 0 |-> 0
294 // 0, 0, 0, 1 |-> 2
295 // 0, 0, 1, 0 |-> 0
296 // 0, 0, 1, 1 |-> 5
297 // 0, 1, 0, 0 |-> 0
298 // 0, 1, 0, 1 |-> 0
299 // 0, 1, 1, 0 |-> 1
300 // 0, 1, 1, 1 |-> 2
301 // 1, 0, 0, 0 |-> 0
302 // 1, 0, 0, 1 |-> 1
303 // 1, 0, 1, 0 |-> 1
304 // 1, 0, 1, 1 |-> 2
305 // 1, 1, 0, 0 |-> 0
306 // 1, 1, 0, 1 |-> 0
307 // 1, 1, 1, 0 |-> 0
308 // 1, 1, 1, 1 |-> 0
309 //
310 let w = merge_polynomials(&[w1, w2, w3])?;
311 // w is [0,2,0,5,0,0,1,2, 0,1,1,2]
312 let w_eval = vec![
313 F::zero(),
314 F::from(2u64),
315 F::zero(),
316 F::from(5u64),
317 F::zero(),
318 F::zero(),
319 F::from(1u64),
320 F::from(2u64),
321 F::zero(),
322 F::one(),
323 F::from(1u64),
324 F::from(2u64),
325 F::zero(),
326 F::zero(),
327 F::zero(),
328 F::zero(),
329 ];
330 let w_rec = DenseMultilinearExtension::from_evaluations_vec(4, w_eval);
331
332 assert_eq!(w, w_rec);
333 }
334 Ok(())
335 }
336
337 #[test]
338 fn test_build_l() -> Result<(), PCSErrors> {
339 test_build_l_helper::<Fr>()
340 }
341
342 fn test_build_l_helper<F: PrimeField>() -> Result<(), PCSErrors> {
343 // point 1 is [1, 2]
344 let point1 = [Fr::from(1u64), Fr::from(2u64)];
345
346 // point 2 is [3, 4]
347 let point2 = [Fr::from(3u64), Fr::from(4u64)];
348
349 // point 3 is [5, 6]
350 let point3 = [Fr::from(5u64), Fr::from(6u64)];
351
352 {
353 let l = build_l(2, &[&point1, &point2])?;
354
355 // roots: [1, -1]
356 // l0 = -1/2 * x + 1/2
357 // l1 = -x + 2
358 // l2 = -x + 3
359 let l0 = DensePolynomial::from_coefficients_vec(vec![
360 Fr::one() / Fr::from(2u64),
361 -Fr::one() / Fr::from(2u64),
362 ]);
363 let l1 = DensePolynomial::from_coefficients_vec(vec![Fr::from(2u64), -Fr::one()]);
364 let l2 = DensePolynomial::from_coefficients_vec(vec![Fr::from(3u64), -Fr::one()]);
365
366 assert_eq!(l0, l[0], "l0 not equal");
367 assert_eq!(l1, l[1], "l1 not equal");
368 assert_eq!(l2, l[2], "l2 not equal");
369 }
370
371 {
372 let l = build_l(2, &[&point1, &point2, &point3])?;
373
374 // sage: q = 52435875175126190479447740508185965837690552500527637822603658699938581184513
375 // sage: P.<x> = PolynomialRing(Zmod(q))
376 // sage: root1 = 1
377 // sage: root2 = 0x8D51CCCE760304D0EC030002760300000001000000000000
378 // sage: root3 = -1
379 // sage: root4 = -root2
380 // Arkwork's code is a bit wired: it also interpolate (root4, 0)
381 // which returns a degree 3 polynomial, instead of degree 2
382
383 // ========================
384 // l0: [0, 0, 1]
385 // ========================
386 // sage: points = [(root1, 0), (root2, 0), (root3, 1), (root4, 0)]
387 // sage: P.lagrange_polynomial(points)
388 // 13108968793781547619861935127046491459422638125131909455650914674984645296128*x^3 +
389 // 39326906381344642859585805381139474378267914375395728366952744024953935888385*x^2 +
390 // 13108968793781547619861935127046491459422638125131909455650914674984645296128*x +
391 // 39326906381344642859585805381139474378267914375395728366952744024953935888385
392 let l0 = DensePolynomial::from_coefficients_vec(vec![
393 field_new!(
394 Fr,
395 "39326906381344642859585805381139474378267914375395728366952744024953935888385"
396 ),
397 field_new!(
398 Fr,
399 "13108968793781547619861935127046491459422638125131909455650914674984645296128"
400 ),
401 field_new!(
402 Fr,
403 "39326906381344642859585805381139474378267914375395728366952744024953935888385"
404 ),
405 field_new!(
406 Fr,
407 "13108968793781547619861935127046491459422638125131909455650914674984645296128"
408 ),
409 ]);
410
411 // ========================
412 // l1: [0, 1, 0]
413 // ========================
414 // sage: points = [(root1, 0), (root2, 1), (root3, 0), (root4, 0)]
415 // sage: P.lagrange_polynomial(points)
416 // 866286206518413079694067382671935694567563117191340490752*x^3 +
417 // 13108968793781547619861935127046491459422638125131909455650914674984645296128*x^2 +
418 // 52435875175126190478581454301667552757996485117855702128036095582747240693761*x +
419 // 39326906381344642859585805381139474378267914375395728366952744024953935888385
420 let l1 = DensePolynomial::from_coefficients_vec(vec![
421 field_new!(
422 Fr,
423 "39326906381344642859585805381139474378267914375395728366952744024953935888385"
424 ),
425 field_new!(
426 Fr,
427 "52435875175126190478581454301667552757996485117855702128036095582747240693761"
428 ),
429 field_new!(
430 Fr,
431 "13108968793781547619861935127046491459422638125131909455650914674984645296128"
432 ),
433 field_new!(
434 Fr,
435 "866286206518413079694067382671935694567563117191340490752"
436 ),
437 ]);
438
439 // ========================
440 // l2: [1, 3, 5]
441 // ========================
442 // sage: points = [(root1, 1), (root2, 3), (root3, 5), (root4, 0)]
443 // sage: P.lagrange_polynomial(points)
444 // 2598858619555239239082202148015807083702689351574021472255*x^3 +
445 // 13108968793781547619861935127046491459422638125131909455650914674984645296129*x^2 +
446 // 52435875175126190476848881888630726598608350352511830738900969348364559712256*x +
447 // 39326906381344642859585805381139474378267914375395728366952744024953935888387
448 let l2 = DensePolynomial::from_coefficients_vec(vec![
449 field_new!(
450 Fr,
451 "39326906381344642859585805381139474378267914375395728366952744024953935888387"
452 ),
453 field_new!(
454 Fr,
455 "52435875175126190476848881888630726598608350352511830738900969348364559712256"
456 ),
457 field_new!(
458 Fr,
459 "13108968793781547619861935127046491459422638125131909455650914674984645296129"
460 ),
461 field_new!(
462 Fr,
463 "2598858619555239239082202148015807083702689351574021472255"
464 ),
465 ]);
466
467 // ========================
468 // l3: [2, 4, 6]
469 // ========================
470 // sage: points = [(root1, 2), (root2, 4), (root3, 6), (root4, 0)]
471 // sage: P.lagrange_polynomial(points)
472 // 3465144826073652318776269530687742778270252468765361963007*x^3 +
473 // x^2 +
474 // 52435875175126190475982595682112313518914282969839895044333406231173219221504*x +
475 // 3
476 let l3 = DensePolynomial::from_coefficients_vec(vec![
477 Fr::from(3u64),
478 field_new!(
479 Fr,
480 "52435875175126190475982595682112313518914282969839895044333406231173219221504"
481 ),
482 Fr::one(),
483 field_new!(
484 Fr,
485 "3465144826073652318776269530687742778270252468765361963007"
486 ),
487 ]);
488
489 assert_eq!(l0, l[0], "l0 not equal");
490 assert_eq!(l1, l[1], "l1 not equal");
491 assert_eq!(l2, l[2], "l2 not equal");
492 assert_eq!(l3, l[3], "l3 not equal");
493 }
494 Ok(())
495 }
496
497 #[test]
498 fn test_qx() -> Result<(), PCSErrors> {
499 // Example from page 53:
500 // W1 = 3x1x2 + 2x2
501 let w_eval = vec![Fr::zero(), Fr::from(2u64), Fr::zero(), Fr::from(5u64)];
502 let w1 = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
503
504 // W2 = x1x2 + x1
505 let w_eval = vec![Fr::zero(), Fr::zero(), Fr::from(1u64), Fr::from(2u64)];
506 let w2 = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
507
508 // W3 = x1 + x2
509 let w_eval = vec![Fr::zero(), Fr::one(), Fr::from(1u64), Fr::from(2u64)];
510 let w3 = DenseMultilinearExtension::from_evaluations_vec(2, w_eval);
511
512 let r = Fr::from(42u64);
513
514 // point 1 is [1, 2]
515 let point1 = [Fr::from(1u64), Fr::from(2u64)];
516
517 // point 2 is [3, 4]
518 let point2 = [Fr::from(3u64), Fr::from(4u64)];
519
520 // point 3 is [5, 6]
521 let point3 = [Fr::from(5u64), Fr::from(6u64)];
522
523 {
524 // w = (3x1x2 + 2x2)(1-x0) + (x1x2 + x1)x0
525 // with evaluations: [0,2,0,5,0,0,1,2]
526 let w = merge_polynomials(&[w1.clone(), w2.clone()])?;
527
528 let l = build_l(2, &[&point1, &point2])?;
529
530 // sage: P.<x> = PolynomialRing(ZZ)
531 // sage: l0 = -1/2 * x + 1/2
532 // sage: l1 = -x + 2
533 // sage: l2 = -x + 3
534 // sage: w = (3 * l1 * l2 + 2 * l2) * (1-l0) + (l1 * l2 + l1) * l0
535 // sage: w
536 // x^3 - 7/2*x^2 - 7/2*x + 16
537 //
538 // q(x) = x^3 - 7/2*x^2 - 7/2*x + 16
539 let q_x = compute_w_circ_l(&w, &l)?;
540
541 let point: Vec<Fr> = l.iter().rev().map(|poly| poly.evaluate(&r)).collect();
542
543 assert_eq!(
544 q_x.evaluate(&r),
545 w.evaluate(&point).unwrap(),
546 "q(r) != w(l(r))"
547 );
548 }
549
550 {
551 // W = (3x1x2 + 2x2) * (1-y1) * (1-y2)
552 // + (x1x2 + x1) * (1-y1) * y2
553 // + (x1 + x2) * y1 * (1-y2)
554 let w = merge_polynomials(&[w1, w2, w3])?;
555
556 let l = build_l(2, &[&point1, &point2, &point3])?;
557
558 // l0 =
559 // 13108968793781547619861935127046491459422638125131909455650914674984645296128*x^3 +
560 // 39326906381344642859585805381139474378267914375395728366952744024953935888385*x^2 +
561 // 13108968793781547619861935127046491459422638125131909455650914674984645296128*x +
562 // 39326906381344642859585805381139474378267914375395728366952744024953935888385
563 //
564 // l1 =
565 // 866286206518413079694067382671935694567563117191340490752*x^3 +
566 // 13108968793781547619861935127046491459422638125131909455650914674984645296128*x^2 +
567 // 52435875175126190478581454301667552757996485117855702128036095582747240693761*x +
568 // 39326906381344642859585805381139474378267914375395728366952744024953935888385
569 //
570 // l2 =
571 // 2598858619555239239082202148015807083702689351574021472255*x^3 +
572 // 13108968793781547619861935127046491459422638125131909455650914674984645296129*x^2 +
573 // 52435875175126190476848881888630726598608350352511830738900969348364559712256*x +
574 // 39326906381344642859585805381139474378267914375395728366952744024953935888387
575 //
576 // l3 =
577 // 3465144826073652318776269530687742778270252468765361963007*x^3 +
578 // x^2 +
579 // 52435875175126190475982595682112313518914282969839895044333406231173219221504*x +
580 // 3
581 //
582 // q_x = (3*l2*l3 + 2*l3) * (1-l0) *(1-l1)
583 // + (l2*l3+l2)*(1-l0)*l1
584 // + (l2+l3)*l0*(1-l1)
585 // q_x(42) = 42675783400755005965526147011103024780845819057955866345013183657072368533932
586 let q_x = compute_w_circ_l(&w, &l)?;
587
588 let point: Vec<Fr> = vec![
589 l[3].evaluate(&r),
590 l[2].evaluate(&r),
591 l[1].evaluate(&r),
592 l[0].evaluate(&r),
593 ];
594
595 assert_eq!(
596 q_x.evaluate(&r),
597 field_new!(
598 Fr,
599 "42675783400755005965526147011103024780845819057955866345013183657072368533932"
600 ),
601 );
602 assert_eq!(
603 q_x.evaluate(&r),
604 w.evaluate(&point).unwrap(),
605 "q(r) != w(l(r))"
606 );
607 }
608 Ok(())
609 }
610}