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

ark_poly/polynomial/univariate/
sparse.rs
