1use super::PolyMultiplier;
19use crate::fft::{EvaluationDomain, Evaluations, Polynomial};
20use snarkvm_fields::{Field, PrimeField};
21use snarkvm_utilities::{cfg_iter_mut, serialize::*};
22
23use anyhow::Result;
24use num_traits::CheckedDiv;
25use rand::Rng;
26use std::{
27 fmt,
28 ops::{Add, AddAssign, Deref, DerefMut, Div, Mul, MulAssign, Neg, Sub, SubAssign},
29};
30
31use itertools::Itertools;
32
33#[cfg(not(feature = "serial"))]
34use rayon::prelude::*;
35
36#[derive(Clone, PartialEq, Eq, Hash, Default, CanonicalSerialize, CanonicalDeserialize)]
38#[must_use]
39pub struct DensePolynomial<F: Field> {
40 pub coeffs: Vec<F>,
42}
43
44impl<F: Field> fmt::Debug for DensePolynomial<F> {
45 fn fmt(&self, f: &mut fmt::Formatter) -> Result<(), fmt::Error> {
46 for (i, coeff) in self.coeffs.iter().enumerate().filter(|(_, c)| !c.is_zero()) {
47 if i == 0 {
48 write!(f, "\n{coeff:?}",)?;
49 } else if i == 1 {
50 write!(f, " + \n{coeff:?} * x")?;
51 } else {
52 write!(f, " + \n{coeff:?} * x^{i}")?;
53 }
54 }
55 Ok(())
56 }
57}
58
59impl<F: Field> DensePolynomial<F> {
60 pub fn zero() -> Self {
62 Self { coeffs: Vec::new() }
63 }
64
65 pub fn is_zero(&self) -> bool {
67 self.coeffs.is_empty() || self.coeffs.iter().all(|coeff| coeff.is_zero())
68 }
69
70 pub fn from_coefficients_slice(coeffs: &[F]) -> Self {
72 Self::from_coefficients_vec(coeffs.to_vec())
73 }
74
75 pub fn from_coefficients_vec(mut coeffs: Vec<F>) -> Self {
77 while let Some(true) = coeffs.last().map(|c| c.is_zero()) {
79 coeffs.pop();
80 }
81 assert!(coeffs.last().map_or(true, |coeff| !coeff.is_zero()));
84
85 Self { coeffs }
86 }
87
88 pub fn degree(&self) -> usize {
90 if self.is_zero() {
91 0
92 } else {
93 assert!(self.coeffs.last().map_or(false, |coeff| !coeff.is_zero()));
94 self.coeffs.len() - 1
95 }
96 }
97
98 pub fn evaluate(&self, point: F) -> F {
100 if self.is_zero() {
101 return F::zero();
102 } else if point.is_zero() {
103 return self.coeffs[0];
104 }
105 let mut powers_of_point = Vec::with_capacity(1 + self.degree());
106 powers_of_point.push(F::one());
107 let mut cur = point;
108 for _ in 0..self.degree() {
109 powers_of_point.push(cur);
110 cur *= point;
111 }
112 let zero = F::zero();
113 let mapping = crate::cfg_into_iter!(powers_of_point).zip_eq(&self.coeffs).map(|(power, coeff)| power * coeff);
114 crate::cfg_reduce!(mapping, || zero, |a, b| a + b)
115 }
116
117 pub fn rand<R: Rng>(d: usize, rng: &mut R) -> Self {
122 let mut random_coeffs = (0..(d + 1)).map(|_| F::rand(rng)).collect_vec();
123 while random_coeffs[d].is_zero() {
124 random_coeffs[d] = F::rand(rng);
126 }
127 Self::from_coefficients_vec(random_coeffs)
128 }
129
130 pub fn coeffs(&self) -> &[F] {
132 &self.coeffs
133 }
134
135 #[cfg(test)]
137 fn naive_mul(&self, other: &Self) -> Self {
138 if self.is_zero() || other.is_zero() {
139 DensePolynomial::zero()
140 } else {
141 let mut result = vec![F::zero(); self.degree() + other.degree() + 1];
142 for (i, self_coeff) in self.coeffs.iter().enumerate() {
143 for (j, other_coeff) in other.coeffs.iter().enumerate() {
144 result[i + j] += *self_coeff * other_coeff;
145 }
146 }
147 DensePolynomial::from_coefficients_vec(result)
148 }
149 }
150}
151
152impl<F: PrimeField> DensePolynomial<F> {
153 pub fn mul_by_vanishing_poly(&self, domain: EvaluationDomain<F>) -> DensePolynomial<F> {
155 let mut shifted = vec![F::zero(); domain.size()];
156 shifted.extend_from_slice(&self.coeffs);
157 crate::cfg_iter_mut!(shifted[..self.coeffs.len()]).zip_eq(&self.coeffs).for_each(|(s, c)| *s -= c);
158 DensePolynomial::from_coefficients_vec(shifted)
159 }
160
161 pub fn divide_by_vanishing_poly(
164 &self,
165 domain: EvaluationDomain<F>,
166 ) -> Result<(DensePolynomial<F>, DensePolynomial<F>)> {
167 let self_poly = Polynomial::from(self);
168 let vanishing_poly = Polynomial::from(domain.vanishing_polynomial());
169 self_poly.divide_with_q_and_r(&vanishing_poly)
170 }
171
172 pub fn evaluate_over_domain_by_ref(&self, domain: EvaluationDomain<F>) -> Evaluations<F> {
174 let poly: Polynomial<'_, F> = self.into();
175 Polynomial::<F>::evaluate_over_domain(poly, domain)
176 }
177
178 pub fn evaluate_over_domain(self, domain: EvaluationDomain<F>) -> Evaluations<F> {
180 let poly: Polynomial<'_, F> = self.into();
181 Polynomial::<F>::evaluate_over_domain(poly, domain)
182 }
183}
184
185impl<F: Field> From<super::SparsePolynomial<F>> for DensePolynomial<F> {
186 fn from(other: super::SparsePolynomial<F>) -> Self {
187 let mut result = vec![F::zero(); other.degree() + 1];
188 for (i, coeff) in other.coeffs() {
189 result[*i] = *coeff;
190 }
191 DensePolynomial::from_coefficients_vec(result)
192 }
193}
194
195impl<'a, F: Field> Add<&'a DensePolynomial<F>> for &'_ DensePolynomial<F> {
196 type Output = DensePolynomial<F>;
197
198 fn add(self, other: &'a DensePolynomial<F>) -> DensePolynomial<F> {
199 let mut result = if self.is_zero() {
200 other.clone()
201 } else if other.is_zero() {
202 self.clone()
203 } else if self.degree() >= other.degree() {
204 let mut result = self.clone();
205 cfg_iter_mut!(result.coeffs).zip(&other.coeffs).for_each(|(a, b)| *a += b);
207 result
208 } else {
209 let mut result = other.clone();
210 cfg_iter_mut!(result.coeffs).zip(&self.coeffs).for_each(|(a, b)| *a += b);
212 result
213 };
214 while let Some(true) = result.coeffs.last().map(|c| c.is_zero()) {
216 result.coeffs.pop();
217 }
218 result
219 }
220}
221
222impl<'a, F: Field> AddAssign<&'a DensePolynomial<F>> for DensePolynomial<F> {
223 fn add_assign(&mut self, other: &'a DensePolynomial<F>) {
224 if self.is_zero() {
225 self.coeffs.clear();
226 self.coeffs.extend_from_slice(&other.coeffs);
227 } else if other.is_zero() {
228 } else if self.degree() >= other.degree() {
230 cfg_iter_mut!(self.coeffs).zip(&other.coeffs).for_each(|(a, b)| *a += b);
232 } else {
233 self.coeffs.resize(other.coeffs.len(), F::zero());
235 cfg_iter_mut!(self.coeffs).zip(&other.coeffs).for_each(|(a, b)| *a += b);
237 }
238 while let Some(true) = self.coeffs.last().map(|c| c.is_zero()) {
240 self.coeffs.pop();
241 }
242 }
243}
244
245impl<'a, F: Field> AddAssign<&'a Polynomial<'a, F>> for DensePolynomial<F> {
246 fn add_assign(&mut self, other: &'a Polynomial<F>) {
247 match other {
248 Polynomial::Sparse(p) => *self += &Self::from(p.clone().into_owned()),
249 Polynomial::Dense(p) => *self += p.as_ref(),
250 }
251 }
252}
253
254impl<'a, F: Field> AddAssign<(F, &'a Polynomial<'a, F>)> for DensePolynomial<F> {
255 fn add_assign(&mut self, (f, other): (F, &'a Polynomial<F>)) {
256 match other {
257 Polynomial::Sparse(p) => *self += (f, &Self::from(p.clone().into_owned())),
258 Polynomial::Dense(p) => *self += (f, p.as_ref()),
259 }
260 }
261}
262
263impl<'a, F: Field> AddAssign<(F, &'a DensePolynomial<F>)> for DensePolynomial<F> {
264 #[allow(clippy::suspicious_op_assign_impl)]
265 fn add_assign(&mut self, (f, other): (F, &'a DensePolynomial<F>)) {
266 if self.is_zero() {
267 self.coeffs.clear();
268 self.coeffs.extend_from_slice(&other.coeffs);
269 self.coeffs.iter_mut().for_each(|c| *c *= &f);
270 } else if other.is_zero() {
271 } else if self.degree() >= other.degree() {
273 cfg_iter_mut!(self.coeffs).zip(&other.coeffs).for_each(|(a, b)| {
275 *a += f * b;
276 });
277 } else {
278 self.coeffs.resize(other.coeffs.len(), F::zero());
280 cfg_iter_mut!(self.coeffs).zip(&other.coeffs).for_each(|(a, b)| {
282 *a += f * b;
283 });
284 }
285 while let Some(true) = self.coeffs.last().map(|c| c.is_zero()) {
287 self.coeffs.pop();
288 }
289 }
290}
291
292impl<F: Field> Neg for DensePolynomial<F> {
293 type Output = DensePolynomial<F>;
294
295 #[inline]
296 fn neg(mut self) -> DensePolynomial<F> {
297 for coeff in &mut self.coeffs {
298 *coeff = -*coeff;
299 }
300 self
301 }
302}
303
304impl<'a, F: Field> Sub<&'a DensePolynomial<F>> for &'_ DensePolynomial<F> {
305 type Output = DensePolynomial<F>;
306
307 #[inline]
308 fn sub(self, other: &'a DensePolynomial<F>) -> DensePolynomial<F> {
309 let mut result = if self.is_zero() {
310 let mut result = other.clone();
311 for coeff in &mut result.coeffs {
312 *coeff = -(*coeff);
313 }
314 result
315 } else if other.is_zero() {
316 self.clone()
317 } else if self.degree() >= other.degree() {
318 let mut result = self.clone();
319 cfg_iter_mut!(result.coeffs).zip(&other.coeffs).for_each(|(a, b)| *a -= b);
321 result
322 } else {
323 let mut result = self.clone();
324 result.coeffs.resize(other.coeffs.len(), F::zero());
325 cfg_iter_mut!(result.coeffs).zip(&other.coeffs).for_each(|(a, b)| {
327 *a -= b;
328 });
329 result
330 };
331 while let Some(true) = result.coeffs.last().map(|c| c.is_zero()) {
333 result.coeffs.pop();
334 }
335 result
336 }
337}
338
339impl<'a, F: Field> SubAssign<&'a DensePolynomial<F>> for DensePolynomial<F> {
340 #[inline]
341 fn sub_assign(&mut self, other: &'a DensePolynomial<F>) {
342 if self.is_zero() {
343 self.coeffs.resize(other.coeffs.len(), F::zero());
344 for (i, coeff) in other.coeffs.iter().enumerate() {
345 self.coeffs[i] -= coeff;
346 }
347 } else if other.is_zero() {
348 } else if self.degree() >= other.degree() {
350 cfg_iter_mut!(self.coeffs).zip(&other.coeffs).for_each(|(a, b)| *a -= b);
352 } else {
353 self.coeffs.resize(other.coeffs.len(), F::zero());
355 cfg_iter_mut!(self.coeffs).zip(&other.coeffs).for_each(|(a, b)| *a -= b);
357 }
358 while let Some(true) = self.coeffs.last().map(|c| c.is_zero()) {
360 self.coeffs.pop();
361 }
362 }
363}
364
365impl<'a, F: Field> AddAssign<&'a super::SparsePolynomial<F>> for DensePolynomial<F> {
366 #[inline]
367 fn add_assign(&mut self, other: &'a super::SparsePolynomial<F>) {
368 if self.degree() < other.degree() {
369 self.coeffs.resize(other.degree() + 1, F::zero());
370 }
371 for (i, b) in other.coeffs() {
372 self.coeffs[*i] += b;
373 }
374 while let Some(true) = self.coeffs.last().map(|c| c.is_zero()) {
376 self.coeffs.pop();
377 }
378 }
379}
380
381impl<'a, F: Field> Sub<&'a super::SparsePolynomial<F>> for DensePolynomial<F> {
382 type Output = Self;
383
384 #[inline]
385 fn sub(mut self, other: &'a super::SparsePolynomial<F>) -> Self::Output {
386 if self.degree() < other.degree() {
387 self.coeffs.resize(other.degree() + 1, F::zero());
388 }
389 for (i, b) in other.coeffs() {
390 self.coeffs[*i] -= b;
391 }
392 while let Some(true) = self.coeffs.last().map(|c| c.is_zero()) {
394 self.coeffs.pop();
395 }
396 self
397 }
398}
399
400impl<'a, F: Field> Div<&'a DensePolynomial<F>> for &'_ DensePolynomial<F> {
401 type Output = DensePolynomial<F>;
402
403 #[inline]
405 fn div(self, divisor: &'a DensePolynomial<F>) -> DensePolynomial<F> {
406 let a: Polynomial<_> = self.into();
407 let b: Polynomial<_> = divisor.into();
408 a.divide_with_q_and_r(&b).expect("division failed").0
409 }
410}
411
412impl<F: Field> Div<DensePolynomial<F>> for DensePolynomial<F> {
413 type Output = DensePolynomial<F>;
414
415 #[inline]
417 fn div(self, divisor: DensePolynomial<F>) -> DensePolynomial<F> {
418 let a: Polynomial<_> = self.into();
419 let b: Polynomial<_> = divisor.into();
420 a.divide_with_q_and_r(&b).expect("division failed").0
421 }
422}
423
424impl<F: Field> CheckedDiv for DensePolynomial<F> {
425 #[inline]
426 fn checked_div(&self, divisor: &DensePolynomial<F>) -> Option<DensePolynomial<F>> {
427 let a: Polynomial<_> = self.into();
428 let b: Polynomial<_> = divisor.into();
429 match a.divide_with_q_and_r(&b) {
430 Ok((divisor, remainder)) => {
431 if remainder.is_zero() {
432 Some(divisor)
433 } else {
434 None
435 }
436 }
437 Err(_) => None,
438 }
439 }
440}
441
442impl<'a, F: PrimeField> Mul<&'a DensePolynomial<F>> for &'_ DensePolynomial<F> {
444 type Output = DensePolynomial<F>;
445
446 #[inline]
447 #[allow(clippy::suspicious_arithmetic_impl)]
448 fn mul(self, other: &'a DensePolynomial<F>) -> DensePolynomial<F> {
449 if self.is_zero() || other.is_zero() {
450 DensePolynomial::zero()
451 } else {
452 let mut m = PolyMultiplier::new();
453 m.add_polynomial_ref(self, "");
454 m.add_polynomial_ref(other, "");
455 m.multiply().unwrap()
456 }
457 }
458}
459
460impl<F: Field> Mul<F> for DensePolynomial<F> {
462 type Output = Self;
463
464 #[inline]
465 fn mul(mut self, other: F) -> Self {
466 self.iter_mut().for_each(|c| *c *= other);
467 self
468 }
469}
470
471impl<F: Field> Mul<F> for &'_ DensePolynomial<F> {
473 type Output = DensePolynomial<F>;
474
475 #[inline]
476 fn mul(self, other: F) -> Self::Output {
477 let result = self.clone();
478 result * other
479 }
480}
481
482impl<F: Field> MulAssign<F> for DensePolynomial<F> {
484 #[allow(clippy::suspicious_arithmetic_impl)]
485 fn mul_assign(&mut self, other: F) {
486 cfg_iter_mut!(self).for_each(|c| *c *= other);
487 }
488}
489
490impl<F: Field> std::iter::Sum for DensePolynomial<F> {
492 fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
493 iter.fold(DensePolynomial::zero(), |a, b| &a + &b)
494 }
495}
496
497impl<F: Field> Deref for DensePolynomial<F> {
498 type Target = [F];
499
500 fn deref(&self) -> &[F] {
501 &self.coeffs
502 }
503}
504
505impl<F: Field> DerefMut for DensePolynomial<F> {
506 fn deref_mut(&mut self) -> &mut [F] {
507 &mut self.coeffs
508 }
509}
510
511#[cfg(test)]
512mod tests {
513 use crate::fft::polynomial::*;
514 use num_traits::CheckedDiv;
515 use snarkvm_curves::bls12_377::Fr;
516 use snarkvm_fields::{Field, One, Zero};
517 use snarkvm_utilities::rand::{TestRng, Uniform};
518
519 use rand::RngCore;
520
521 #[test]
522 fn double_polynomials_random() {
523 let rng = &mut TestRng::default();
524 for degree in 0..70 {
525 let p = DensePolynomial::<Fr>::rand(degree, rng);
526 let p_double = &p + &p;
527 let p_quad = &p_double + &p_double;
528 assert_eq!(&(&(&p + &p) + &p) + &p, p_quad);
529 }
530 }
531
532 #[test]
533 fn add_polynomials() {
534 let rng = &mut TestRng::default();
535 for a_degree in 0..70 {
536 for b_degree in 0..70 {
537 let p1 = DensePolynomial::<Fr>::rand(a_degree, rng);
538 let p2 = DensePolynomial::<Fr>::rand(b_degree, rng);
539 let res1 = &p1 + &p2;
540 let res2 = &p2 + &p1;
541 assert_eq!(res1, res2);
542 }
543 }
544 }
545
546 #[test]
547 fn add_polynomials_with_mul() {
548 let rng = &mut TestRng::default();
549 for a_degree in 0..70 {
550 for b_degree in 0..70 {
551 let mut p1 = DensePolynomial::rand(a_degree, rng);
552 let p2 = DensePolynomial::rand(b_degree, rng);
553 let f = Fr::rand(rng);
554 let f_p2 = DensePolynomial::from_coefficients_vec(p2.coeffs.iter().map(|c| f * c).collect());
555 let res2 = &f_p2 + &p1;
556 p1 += (f, &p2);
557 let res1 = p1;
558 assert_eq!(res1, res2);
559 }
560 }
561 }
562
563 #[test]
564 fn sub_polynomials() {
565 let rng = &mut TestRng::default();
566 let p1 = DensePolynomial::<Fr>::rand(5, rng);
567 let p2 = DensePolynomial::<Fr>::rand(3, rng);
568 let res1 = &p1 - &p2;
569 let res2 = &p2 - &p1;
570 assert_eq!(&res1 + &p2, p1, "Subtraction should be inverse of addition!");
571 assert_eq!(res1, -res2, "p2 - p1 = -(p1 - p2)");
572 }
573
574 #[test]
575 fn divide_polynomials_fixed() {
576 let dividend = DensePolynomial::from_coefficients_slice(&[
577 "4".parse().unwrap(),
578 "8".parse().unwrap(),
579 "5".parse().unwrap(),
580 "1".parse().unwrap(),
581 ]);
582 let divisor = DensePolynomial::from_coefficients_slice(&[Fr::one(), Fr::one()]); let result = dividend.checked_div(&divisor).unwrap();
584 let expected_result = DensePolynomial::from_coefficients_slice(&[
585 "4".parse().unwrap(),
586 "4".parse().unwrap(),
587 "1".parse().unwrap(),
588 ]);
589 assert_eq!(expected_result, result);
590 }
591
592 #[test]
593 #[allow(clippy::needless_borrow)]
594 fn divide_polynomials_random() {
595 let rng = &mut TestRng::default();
596
597 for a_degree in 0..70 {
598 for b_degree in 0..70 {
599 let dividend = DensePolynomial::<Fr>::rand(a_degree, rng);
600 let divisor = DensePolynomial::<Fr>::rand(b_degree, rng);
601 let (quotient, remainder) =
602 Polynomial::divide_with_q_and_r(&(÷nd).into(), &(&divisor).into()).unwrap();
603 assert_eq!(dividend, &(&divisor * "ient) + &remainder)
604 }
605 }
606 }
607
608 #[test]
609 fn evaluate_polynomials() {
610 let rng = &mut TestRng::default();
611 for a_degree in 0..70 {
612 let p = DensePolynomial::rand(a_degree, rng);
613 let point: Fr = Fr::from(10u64);
614 let mut total = Fr::zero();
615 for (i, coeff) in p.coeffs.iter().enumerate() {
616 total += point.pow([i as u64]) * coeff;
617 }
618 assert_eq!(p.evaluate(point), total);
619 }
620 }
621
622 #[test]
623 fn divide_poly_by_zero() {
624 let a = Polynomial::<Fr>::zero();
625 let b = Polynomial::<Fr>::zero();
626 assert!(a.divide_with_q_and_r(&b).is_err());
627 }
628
629 #[test]
630 fn mul_polynomials_random() {
631 let rng = &mut TestRng::default();
632 for a_degree in 0..70 {
633 for b_degree in 0..70 {
634 dbg!(a_degree);
635 dbg!(b_degree);
636 let a = DensePolynomial::<Fr>::rand(a_degree, rng);
637 let b = DensePolynomial::<Fr>::rand(b_degree, rng);
638 assert_eq!(&a * &b, a.naive_mul(&b))
639 }
640 }
641 }
642
643 #[test]
644 fn mul_polynomials_n_random() {
645 let rng = &mut TestRng::default();
646
647 let max_degree = 1 << 8;
648
649 for _ in 0..10 {
650 let mut multiplier = PolyMultiplier::new();
651 let a = DensePolynomial::<Fr>::rand(max_degree / 2, rng);
652 let mut mul_degree = a.degree() + 1;
653 multiplier.add_polynomial(a.clone(), "a");
654 let mut naive = a.clone();
655
656 let num_polys = (rng.next_u32() as usize) % 8;
658 let num_evals = (rng.next_u32() as usize) % 4;
659 println!("\nnum_polys {num_polys}, num_evals {num_evals}");
660
661 for _ in 1..num_polys {
662 let degree = (rng.next_u32() as usize) % max_degree;
663 mul_degree += degree + 1;
664 println!("poly degree {degree}");
665 let a = DensePolynomial::<Fr>::rand(degree, rng);
666 naive = naive.naive_mul(&a);
667 multiplier.add_polynomial(a.clone(), "a");
668 }
669
670 let mut eval_degree = mul_degree;
672 let domain = EvaluationDomain::new(mul_degree).unwrap();
673 println!("mul_degree {}, domain {}", mul_degree, domain.size());
674 for _ in 0..num_evals {
675 let a = DensePolynomial::<Fr>::rand(mul_degree / 8, rng);
676 eval_degree += a.len() + 1;
677 if eval_degree < domain.size() {
678 println!("eval degree {eval_degree}");
679 let mut a_evals = a.clone().coeffs;
680 domain.fft_in_place(&mut a_evals);
681 let a_evals = Evaluations::from_vec_and_domain(a_evals, domain);
682
683 naive = naive.naive_mul(&a);
684 multiplier.add_evaluation(a_evals, "a");
685 }
686 }
687
688 assert_eq!(multiplier.multiply().unwrap(), naive);
689 }
690 }
691
692 #[test]
693 fn mul_polynomials_corner_cases() {
694 let rng = &mut TestRng::default();
695
696 let a_degree = 70;
697
698 println!("Test single polynomial");
700 let a = DensePolynomial::<Fr>::rand(a_degree, rng);
701 let mut multiplier = PolyMultiplier::new();
702 multiplier.add_polynomial(a.clone(), "a");
703 assert_eq!(multiplier.multiply().unwrap(), a);
704
705 }
707
708 #[test]
709 fn mul_by_vanishing_poly() {
710 let rng = &mut TestRng::default();
711 for size in 1..10 {
712 let domain = EvaluationDomain::new(1 << size).unwrap();
713 for degree in 0..70 {
714 let p = DensePolynomial::<Fr>::rand(degree, rng);
715 let ans1 = p.mul_by_vanishing_poly(domain);
716 let ans2 = &p * &domain.vanishing_polynomial().into();
717 assert_eq!(ans1, ans2);
718 }
719 }
720 }
721}