ark_poly/polynomial/univariate/
sparse.rs

1//! A sparse polynomial represented in coefficient form.
2use crate::{
3    polynomial::Polynomial,
4    univariate::{DenseOrSparsePolynomial, DensePolynomial},
5    DenseUVPolynomial, EvaluationDomain, Evaluations,
6};
7use ark_ff::{FftField, Field, Zero};
8use ark_serialize::{CanonicalDeserialize, CanonicalSerialize};
9use ark_std::{
10    cfg_iter_mut,
11    cmp::Ordering,
12    collections::BTreeMap,
13    fmt,
14    ops::{Add, AddAssign, Deref, DerefMut, Mul, Neg, SubAssign},
15    vec,
16    vec::*,
17};
18
19#[cfg(feature = "parallel")]
20use rayon::prelude::*;
21
22/// Stores a sparse polynomial in coefficient form.
23#[derive(Clone, PartialEq, Eq, Hash, Default, CanonicalSerialize, CanonicalDeserialize)]
24pub struct SparsePolynomial<F: Field> {
25    /// The coefficient a_i of `x^i` is stored as (i, a_i) in `self.coeffs`.
26    /// the entries in `self.coeffs` *must*  be sorted in increasing order of
27    /// `i`.
28    pub coeffs: Vec<(usize, F)>,
29}
30
31impl<F: Field> fmt::Debug for SparsePolynomial<F> {
32    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
33        for (i, coeff) in self.coeffs.iter().filter(|(_, c)| !c.is_zero()) {
34            if *i == 0 {
35                write!(f, "\n{coeff:?}")?;
36            } else if *i == 1 {
37                write!(f, " + \n{coeff:?} * x")?;
38            } else {
39                write!(f, " + \n{coeff:?} * x^{i}")?;
40            }
41        }
42        Ok(())
43    }
44}
45
46impl<F: Field> Deref for SparsePolynomial<F> {
47    type Target = [(usize, F)];
48
49    fn deref(&self) -> &[(usize, F)] {
50        &self.coeffs
51    }
52}
53
54impl<F: Field> DerefMut for SparsePolynomial<F> {
55    fn deref_mut(&mut self) -> &mut [(usize, F)] {
56        &mut self.coeffs
57    }
58}
59
60impl<F: Field> Polynomial<F> for SparsePolynomial<F> {
61    type Point = F;
62
63    /// Returns the degree of the polynomial.
64    fn degree(&self) -> usize {
65        if self.is_zero() {
66            0
67        } else {
68            assert!(self.coeffs.last().is_some_and(|(_, c)| !c.is_zero()));
69            self.coeffs.last().unwrap().0
70        }
71    }
72
73    /// Evaluates `self` at the given `point` in the field.
74    fn evaluate(&self, point: &F) -> F {
75        if self.is_zero() {
76            return F::zero();
77        }
78
79        // We need floor(log2(deg)) + 1 powers, starting from the 0th power p^2^0 = p
80        let num_powers = 0usize.leading_zeros() - self.degree().leading_zeros();
81        let mut powers_of_2 = Vec::with_capacity(num_powers as usize);
82
83        let mut p = *point;
84        powers_of_2.push(p);
85        for _ in 1..num_powers {
86            p.square_in_place();
87            powers_of_2.push(p);
88        }
89        // compute all coeff * point^{i} and then sum the results
90        let total = self
91            .coeffs
92            .iter()
93            .map(|(i, c)| {
94                debug_assert_eq!(
95                    F::pow_with_table(&powers_of_2[..], [*i as u64]).unwrap(),
96                    point.pow([*i as u64]),
97                    "pows not equal"
98                );
99                *c * F::pow_with_table(&powers_of_2[..], [*i as u64]).unwrap()
100            })
101            .sum();
102        total
103    }
104}
105
106impl<F: Field> Add for SparsePolynomial<F> {
107    type Output = Self;
108
109    fn add(self, other: Self) -> Self {
110        &self + &other
111    }
112}
113
114impl<'a, F: Field> Add<&'a SparsePolynomial<F>> for &SparsePolynomial<F> {
115    type Output = SparsePolynomial<F>;
116
117    fn add(self, other: &'a SparsePolynomial<F>) -> SparsePolynomial<F> {
118        if self.is_zero() {
119            return other.clone();
120        } else if other.is_zero() {
121            return self.clone();
122        }
123        // Single pass add algorithm (merging two sorted sets)
124        let mut result = SparsePolynomial::<F>::zero();
125        // our current index in each vector
126        let mut self_index = 0;
127        let mut other_index = 0;
128        loop {
129            // if we've reached the end of one vector, just append the other vector to our
130            // result.
131            if self_index == self.coeffs.len() && other_index == other.coeffs.len() {
132                return result;
133            } else if self_index == self.coeffs.len() {
134                result.append_coeffs(&other.coeffs[other_index..]);
135                return result;
136            } else if other_index == other.coeffs.len() {
137                result.append_coeffs(&self.coeffs[self_index..]);
138                return result;
139            }
140
141            // Get the current degree / coeff for each
142            let (self_term_degree, self_term_coeff) = self.coeffs[self_index];
143            let (other_term_degree, other_term_coeff) = other.coeffs[other_index];
144            // add the lower degree term to our sorted set.
145            match self_term_degree.cmp(&other_term_degree) {
146                Ordering::Less => {
147                    result.coeffs.push((self_term_degree, self_term_coeff));
148                    self_index += 1;
149                },
150                Ordering::Equal => {
151                    let term_sum = self_term_coeff + other_term_coeff;
152                    if !term_sum.is_zero() {
153                        result.coeffs.push((self_term_degree, term_sum));
154                    }
155                    self_index += 1;
156                    other_index += 1;
157                },
158                Ordering::Greater => {
159                    result.coeffs.push((other_term_degree, other_term_coeff));
160                    other_index += 1;
161                },
162            }
163        }
164    }
165}
166
167impl<'a, F: Field> AddAssign<&'a Self> for SparsePolynomial<F> {
168    // TODO: Reduce number of clones
169    fn add_assign(&mut self, other: &'a Self) {
170        self.coeffs = (self.clone() + other.clone()).coeffs;
171    }
172}
173
174impl<'a, F: Field> AddAssign<(F, &'a Self)> for SparsePolynomial<F> {
175    // TODO: Reduce number of clones
176    fn add_assign(&mut self, (f, other): (F, &'a Self)) {
177        self.coeffs = (self.clone() + other.clone()).coeffs;
178        for i in 0..self.coeffs.len() {
179            self.coeffs[i].1 *= f;
180        }
181    }
182}
183
184impl<F: Field> Neg for SparsePolynomial<F> {
185    type Output = Self;
186
187    #[inline]
188    fn neg(mut self) -> Self {
189        for (_, coeff) in &mut self.coeffs {
190            *coeff = -*coeff;
191        }
192        self
193    }
194}
195
196impl<'a, F: Field> SubAssign<&'a Self> for SparsePolynomial<F> {
197    // TODO: Reduce number of clones
198    #[inline]
199    fn sub_assign(&mut self, other: &'a Self) {
200        let self_copy = -self.clone();
201        self.coeffs = (self_copy + other.clone()).coeffs;
202    }
203}
204
205impl<F: Field> Mul<F> for &SparsePolynomial<F> {
206    type Output = SparsePolynomial<F>;
207
208    #[inline]
209    fn mul(self, elem: F) -> SparsePolynomial<F> {
210        if self.is_zero() || elem.is_zero() {
211            SparsePolynomial::zero()
212        } else {
213            let mut result = self.clone();
214            cfg_iter_mut!(result).for_each(|e| {
215                e.1 *= elem;
216            });
217            result
218        }
219    }
220}
221
222impl<F: Field> Zero for SparsePolynomial<F> {
223    /// Returns the zero polynomial.
224    fn zero() -> Self {
225        Self { coeffs: Vec::new() }
226    }
227
228    /// Checks if the given polynomial is zero.
229    fn is_zero(&self) -> bool {
230        self.coeffs.is_empty() || self.coeffs.iter().all(|(_, c)| c.is_zero())
231    }
232}
233
234impl<F: Field> SparsePolynomial<F> {
235    /// Constructs a new polynomial from a list of coefficients.
236    pub fn from_coefficients_slice(coeffs: &[(usize, F)]) -> Self {
237        Self::from_coefficients_vec(coeffs.to_vec())
238    }
239
240    /// Constructs a new polynomial from a list of coefficients.
241    /// The function does not combine like terms and so multiple monomials
242    /// of the same degree are ignored.
243    pub fn from_coefficients_vec(mut coeffs: Vec<(usize, F)>) -> Self {
244        // While there are zeros at the end of the coefficient vector, pop them off.
245        while coeffs.last().is_some_and(|(_, c)| c.is_zero()) {
246            coeffs.pop();
247        }
248        // Ensure that coeffs are in ascending order.
249        coeffs.sort_by(|(c1, _), (c2, _)| c1.cmp(c2));
250        // Check that either the coefficients vec is empty or that the last coeff is
251        // non-zero.
252        assert!(coeffs.last().map_or(true, |(_, c)| !c.is_zero()));
253
254        Self { coeffs }
255    }
256
257    /// Perform a naive n^2 multiplication of `self` by `other`.
258    pub fn mul(&self, other: &Self) -> Self {
259        if self.is_zero() || other.is_zero() {
260            Self::zero()
261        } else {
262            let mut result = BTreeMap::new();
263            for (i, self_coeff) in &self.coeffs {
264                for (j, other_coeff) in &other.coeffs {
265                    result
266                        .entry(i + j)
267                        .and_modify(|cur_coeff| *cur_coeff += *self_coeff * other_coeff)
268                        .or_insert_with(|| *self_coeff * other_coeff);
269                }
270            }
271            Self::from_coefficients_vec(result.into_iter().collect())
272        }
273    }
274
275    /// Returns the quotient of the division of `self` of degree n by `other` of k values using an algorithm in O(nk)
276    pub fn div(&self, other: &Self) -> DensePolynomial<F> {
277        let dividend: DenseOrSparsePolynomial<'_, F> = self.into();
278        let divisor: DenseOrSparsePolynomial<'_, F> = other.into();
279
280        dividend.naive_div(&divisor).expect("division failed").0
281    }
282
283    // append append_coeffs to self.
284    // Correctness relies on the lowest degree term in append_coeffs
285    // being higher than self.degree()
286    fn append_coeffs(&mut self, append_coeffs: &[(usize, F)]) {
287        assert!(append_coeffs.is_empty() || self.degree() < append_coeffs[0].0);
288        self.coeffs.extend_from_slice(append_coeffs);
289    }
290}
291
292impl<F: FftField> SparsePolynomial<F> {
293    /// Evaluate `self` over `domain`.
294    pub fn evaluate_over_domain_by_ref<D: EvaluationDomain<F>>(
295        &self,
296        domain: D,
297    ) -> Evaluations<F, D> {
298        let poly: DenseOrSparsePolynomial<'_, F> = self.into();
299        DenseOrSparsePolynomial::evaluate_over_domain(poly, domain)
300    }
301
302    /// Evaluate `self` over `domain`.
303    pub fn evaluate_over_domain<D: EvaluationDomain<F>>(self, domain: D) -> Evaluations<F, D> {
304        let poly: DenseOrSparsePolynomial<'_, F> = self.into();
305        DenseOrSparsePolynomial::evaluate_over_domain(poly, domain)
306    }
307}
308
309impl<F: Field> From<SparsePolynomial<F>> for DensePolynomial<F> {
310    fn from(other: SparsePolynomial<F>) -> Self {
311        let mut result = vec![F::zero(); other.degree() + 1];
312        for (i, coeff) in other.coeffs {
313            result[i] = coeff;
314        }
315        Self::from_coefficients_vec(result)
316    }
317}
318
319impl<F: Field> From<DensePolynomial<F>> for SparsePolynomial<F> {
320    fn from(dense_poly: DensePolynomial<F>) -> Self {
321        Self::from_coefficients_vec(
322            dense_poly
323                .coeffs()
324                .iter()
325                .enumerate()
326                .filter(|&(_, coeff)| !coeff.is_zero())
327                .map(|(i, coeff)| (i, *coeff))
328                .collect(),
329        )
330    }
331}
332
333#[cfg(test)]
334mod tests {
335    use crate::{
336        polynomial::Polynomial,
337        univariate::{DensePolynomial, SparsePolynomial},
338        EvaluationDomain, GeneralEvaluationDomain,
339    };
340    use ark_ff::{UniformRand, Zero};
341    use ark_std::{cmp::max, ops::Mul, rand::Rng, test_rng, vec};
342    use ark_test_curves::bls12_381::Fr;
343
344    // probability of rand sparse polynomial having a particular coefficient be 0
345    const ZERO_COEFF_PROBABILITY: f64 = 0.8f64;
346
347    fn rand_sparse_poly<R: Rng>(degree: usize, rng: &mut R) -> SparsePolynomial<Fr> {
348        // Initialize coeffs so that its guaranteed to have a x^{degree} term
349        let mut coeffs = vec![(degree, Fr::rand(rng))];
350        for i in 0..degree {
351            if !rng.gen_bool(ZERO_COEFF_PROBABILITY) {
352                coeffs.push((i, Fr::rand(rng)));
353            }
354        }
355        SparsePolynomial::from_coefficients_vec(coeffs)
356    }
357
358    #[test]
359    fn evaluate_at_point() {
360        let mut rng = test_rng();
361        // Test evaluation at point by comparing against DensePolynomial
362        for degree in 0..60 {
363            let sparse_poly = rand_sparse_poly(degree, &mut rng);
364            let dense_poly: DensePolynomial<Fr> = sparse_poly.clone().into();
365            let pt = Fr::rand(&mut rng);
366            assert_eq!(sparse_poly.evaluate(&pt), dense_poly.evaluate(&pt));
367        }
368    }
369
370    #[test]
371    fn add_polynomial() {
372        // Test adding polynomials by comparing against dense polynomial
373        let mut rng = test_rng();
374        for degree_a in 0..20 {
375            let sparse_poly_a = rand_sparse_poly(degree_a, &mut rng);
376            let dense_poly_a: DensePolynomial<Fr> = sparse_poly_a.clone().into();
377            for degree_b in 0..20 {
378                let sparse_poly_b = rand_sparse_poly(degree_b, &mut rng);
379                let dense_poly_b: DensePolynomial<Fr> = sparse_poly_b.clone().into();
380
381                // Test Add trait
382                let sparse_sum = sparse_poly_a.clone() + sparse_poly_b.clone();
383                assert_eq!(
384                    sparse_sum.degree(),
385                    max(degree_a, degree_b),
386                    "degree_a = {}, degree_b = {}",
387                    degree_a,
388                    degree_b
389                );
390                let actual_dense_sum: DensePolynomial<Fr> = sparse_sum.into();
391                let expected_dense_sum = dense_poly_a.clone() + dense_poly_b;
392                assert_eq!(
393                    actual_dense_sum, expected_dense_sum,
394                    "degree_a = {}, degree_b = {}",
395                    degree_a, degree_b
396                );
397                // Test AddAssign Trait
398                let mut sparse_add_assign_sum = sparse_poly_a.clone();
399                sparse_add_assign_sum += &sparse_poly_b;
400                let actual_add_assign_dense_sum: DensePolynomial<Fr> = sparse_add_assign_sum.into();
401                assert_eq!(
402                    actual_add_assign_dense_sum, expected_dense_sum,
403                    "degree_a = {}, degree_b = {}",
404                    degree_a, degree_b
405                );
406            }
407        }
408    }
409
410    #[test]
411    fn polynomial_additive_identity() {
412        // Test adding polynomials with its negative equals 0
413        let mut rng = test_rng();
414        for degree in 0..70 {
415            // Test with Neg trait
416            let sparse_poly = rand_sparse_poly(degree, &mut rng);
417            let neg = -sparse_poly.clone();
418            assert!((sparse_poly + neg).is_zero());
419
420            // Test with SubAssign trait
421            let sparse_poly = rand_sparse_poly(degree, &mut rng);
422            let mut result = sparse_poly.clone();
423            result -= &sparse_poly;
424            assert!(result.is_zero());
425        }
426    }
427
428    #[test]
429    fn mul_random_element() {
430        let rng = &mut test_rng();
431        for degree in 0..20 {
432            let a = rand_sparse_poly(degree, rng);
433            let e = Fr::rand(rng);
434            assert_eq!(
435                &a * e,
436                a.mul(&SparsePolynomial::from_coefficients_slice(&[(0, e)]))
437            )
438        }
439    }
440
441    #[test]
442    fn mul_polynomial() {
443        // Test multiplying polynomials over their domains, and over the native
444        // representation. The expected result is obtained by comparing against
445        // dense polynomial
446        let mut rng = test_rng();
447        for degree_a in 0..20 {
448            let sparse_poly_a = rand_sparse_poly(degree_a, &mut rng);
449            let dense_poly_a: DensePolynomial<Fr> = sparse_poly_a.clone().into();
450            for degree_b in 0..20 {
451                let sparse_poly_b = rand_sparse_poly(degree_b, &mut rng);
452                let dense_poly_b: DensePolynomial<Fr> = sparse_poly_b.clone().into();
453
454                // Test multiplying the polynomials over their native representation
455                let sparse_prod = sparse_poly_a.mul(&sparse_poly_b);
456                assert_eq!(
457                    sparse_prod.degree(),
458                    degree_a + degree_b,
459                    "degree_a = {}, degree_b = {}",
460                    degree_a,
461                    degree_b
462                );
463                let dense_prod = dense_poly_a.naive_mul(&dense_poly_b);
464                assert_eq!(sparse_prod.degree(), dense_prod.degree());
465                assert_eq!(
466                    sparse_prod,
467                    SparsePolynomial::<Fr>::from(dense_prod),
468                    "degree_a = {}, degree_b = {}",
469                    degree_a,
470                    degree_b
471                );
472
473                // Test multiplying the polynomials over their evaluations and interpolating
474                let domain = GeneralEvaluationDomain::new(sparse_prod.degree() + 1).unwrap();
475                let poly_a_evals = sparse_poly_a.evaluate_over_domain_by_ref(domain);
476                let poly_b_evals = sparse_poly_b.evaluate_over_domain_by_ref(domain);
477                let poly_prod_evals = sparse_prod.evaluate_over_domain_by_ref(domain);
478                assert_eq!(poly_a_evals.mul(&poly_b_evals), poly_prod_evals);
479            }
480        }
481    }
482
483    #[test]
484    fn div_polynomial() {
485        let mut rng = test_rng();
486        for degree_a in 0..10 {
487            let sparse_poly_a = rand_sparse_poly(degree_a, &mut rng);
488            let dense_poly_a: DensePolynomial<Fr> = sparse_poly_a.clone().into();
489            for degree_b in 0..degree_a {
490                let sparse_poly_b = rand_sparse_poly(degree_b, &mut rng);
491                let dense_poly_b: DensePolynomial<Fr> = sparse_poly_b.clone().into();
492
493                // Test dividing the polynomials over their native representation
494                let sparse_quotient = sparse_poly_a.div(&sparse_poly_b);
495                assert_eq!(
496                    sparse_quotient.degree(),
497                    degree_a - degree_b,
498                    "degree_a  = {}, degree_b = {}, degree_q = {}",
499                    degree_a - degree_b,
500                    degree_b,
501                    sparse_quotient.degree(),
502                );
503
504                // Test that both dense and sparse division will return the same value
505                // Since the FFT division is tested in the dense.rs, we then know this division is valid too
506                // (we are over FFT fields so dense div will be FFT)
507                let dense_quotient = &dense_poly_a / &dense_poly_b;
508                assert_eq!(sparse_quotient.degree(), dense_quotient.degree());
509                assert_eq!(
510                    sparse_quotient,
511                    dense_quotient,
512                    "sparse quotient different from dense quotient for the same division! sparse {:?}, dense {:?}",
513                    sparse_quotient,
514                    dense_quotient,
515                );
516            }
517        }
518    }
519
520    #[test]
521    fn evaluate_over_domain() {
522        // Test that polynomial evaluation over a domain, and interpolation returns the
523        // same poly.
524        let mut rng = test_rng();
525        for poly_degree_dim in 0..5 {
526            let poly_degree = (1 << poly_degree_dim) - 1;
527            let sparse_poly = rand_sparse_poly(poly_degree, &mut rng);
528
529            for domain_dim in poly_degree_dim..(poly_degree_dim + 2) {
530                let domain_size = 1 << domain_dim;
531                let domain = GeneralEvaluationDomain::new(domain_size).unwrap();
532
533                let sparse_evals = sparse_poly.evaluate_over_domain_by_ref(domain);
534
535                // Test interpolation works, by checking against DensePolynomial
536                let dense_poly: DensePolynomial<Fr> = sparse_poly.clone().into();
537                let dense_evals = dense_poly.clone().evaluate_over_domain(domain);
538                assert_eq!(
539                    sparse_evals.clone().interpolate(),
540                    dense_evals.clone().interpolate(),
541                    "poly_degree_dim = {}, domain_dim = {}",
542                    poly_degree_dim,
543                    domain_dim
544                );
545                assert_eq!(
546                    sparse_evals.interpolate(),
547                    dense_poly,
548                    "poly_degree_dim = {}, domain_dim = {}",
549                    poly_degree_dim,
550                    domain_dim
551                );
552                // Consistency check that the dense polynomials interpolation is correct.
553                assert_eq!(
554                    dense_evals.interpolate(),
555                    dense_poly,
556                    "poly_degree_dim = {}, domain_dim = {}",
557                    poly_degree_dim,
558                    domain_dim
559                );
560            }
561        }
562    }
563
564    #[test]
565    fn evaluate_over_small_domain() {
566        // Test that polynomial evaluation over a domain, and interpolation returns the
567        // same poly.
568        let mut rng = test_rng();
569        for poly_degree_dim in 1..5 {
570            let poly_degree = (1 << poly_degree_dim) - 1;
571            let sparse_poly = rand_sparse_poly(poly_degree, &mut rng);
572
573            for domain_dim in 0..poly_degree_dim {
574                let domain_size = 1 << domain_dim;
575                let domain = GeneralEvaluationDomain::new(domain_size).unwrap();
576
577                let sparse_evals = sparse_poly.evaluate_over_domain_by_ref(domain);
578
579                // Test that sparse evaluation and dense evaluation agree
580                let dense_poly: DensePolynomial<Fr> = sparse_poly.clone().into();
581                let dense_evals = dense_poly.clone().evaluate_over_domain(domain);
582                assert_eq!(
583                    sparse_evals, dense_evals,
584                    "poly_degree_dim = {}, domain_dim = {}",
585                    poly_degree_dim, domain_dim
586                );
587
588                // Test interpolation works, by checking that interpolated polynomial agrees with the original on the domain
589                let (_q, r) = (dense_poly.clone() + -sparse_evals.interpolate())
590                    .divide_by_vanishing_poly(domain);
591                assert_eq!(
592                    r,
593                    DensePolynomial::<Fr>::zero(),
594                    "poly_degree_dim = {}, domain_dim = {}",
595                    poly_degree_dim,
596                    domain_dim
597                );
598
599                // Consistency check that the dense polynomials interpolation is correct.
600                let (_q, r) = (dense_poly.clone() + -dense_evals.interpolate())
601                    .divide_by_vanishing_poly(domain);
602                assert_eq!(
603                    r,
604                    DensePolynomial::<Fr>::zero(),
605                    "poly_degree_dim = {}, domain_dim = {}",
606                    poly_degree_dim,
607                    domain_dim
608                );
609            }
610        }
611    }
612}
ark_poly/polynomial/univariate/sparse.rs

ark_poly/polynomial/univariate/
sparse.rs
