ark_poly/polynomial/univariate/
mod.rs

1//! Work with sparse and dense polynomials.
2
3use core::cmp::min;
4
5use crate::{DenseUVPolynomial, EvaluationDomain, Evaluations, Polynomial};
6use ark_ff::{FftField, Field, Zero};
7use ark_std::{borrow::Cow, cfg_iter_mut, vec, vec::*};
8use DenseOrSparsePolynomial::{DPolynomial, SPolynomial};
9
10mod dense;
11mod sparse;
12
13pub use dense::DensePolynomial;
14pub use sparse::SparsePolynomial;
15
16#[cfg(feature = "parallel")]
17use rayon::prelude::*;
18
19/// Represents either a sparse polynomial or a dense one.
20#[derive(Clone)]
21pub enum DenseOrSparsePolynomial<'a, F: Field> {
22    /// Represents the case where `self` is a sparse polynomial
23    SPolynomial(Cow<'a, SparsePolynomial<F>>),
24    /// Represents the case where `self` is a dense polynomial
25    DPolynomial(Cow<'a, DensePolynomial<F>>),
26}
27
28impl<'a, F: 'a + Field> From<DensePolynomial<F>> for DenseOrSparsePolynomial<'a, F> {
29    fn from(other: DensePolynomial<F>) -> Self {
30        DPolynomial(Cow::Owned(other))
31    }
32}
33
34impl<'a, F: 'a + Field> From<&'a DensePolynomial<F>> for DenseOrSparsePolynomial<'a, F> {
35    fn from(other: &'a DensePolynomial<F>) -> Self {
36        DPolynomial(Cow::Borrowed(other))
37    }
38}
39
40impl<'a, F: 'a + Field> From<SparsePolynomial<F>> for DenseOrSparsePolynomial<'a, F> {
41    fn from(other: SparsePolynomial<F>) -> Self {
42        SPolynomial(Cow::Owned(other))
43    }
44}
45
46impl<'a, F: Field> From<&'a SparsePolynomial<F>> for DenseOrSparsePolynomial<'a, F> {
47    fn from(other: &'a SparsePolynomial<F>) -> Self {
48        SPolynomial(Cow::Borrowed(other))
49    }
50}
51
52impl<'a, F: Field> From<DenseOrSparsePolynomial<'a, F>> for DensePolynomial<F> {
53    fn from(other: DenseOrSparsePolynomial<'a, F>) -> Self {
54        match other {
55            DPolynomial(p) => p.into_owned(),
56            SPolynomial(p) => p.into_owned().into(),
57        }
58    }
59}
60
61impl<'a, F: 'a + Field> TryInto<SparsePolynomial<F>> for DenseOrSparsePolynomial<'a, F> {
62    type Error = ();
63
64    fn try_into(self) -> Result<SparsePolynomial<F>, ()> {
65        match self {
66            SPolynomial(p) => Ok(p.into_owned()),
67            _ => Err(()),
68        }
69    }
70}
71
72impl<F: Field> DenseOrSparsePolynomial<'_, F> {
73    /// Checks if the given polynomial is zero.
74    pub fn is_zero(&self) -> bool {
75        match self {
76            SPolynomial(s) => s.is_zero(),
77            DPolynomial(d) => d.is_zero(),
78        }
79    }
80
81    /// Return the degree of `self`.
82    pub fn degree(&self) -> usize {
83        match self {
84            SPolynomial(s) => s.degree(),
85            DPolynomial(d) => d.degree(),
86        }
87    }
88
89    #[inline]
90    fn leading_coefficient(&self) -> Option<&F> {
91        match self {
92            SPolynomial(p) => p.last().map(|(_, c)| c),
93            DPolynomial(p) => p.last(),
94        }
95    }
96
97    #[inline]
98    fn iter_with_index(&self) -> Vec<(usize, F)> {
99        match self {
100            SPolynomial(p) => p.to_vec(),
101            DPolynomial(p) => p.iter().copied().enumerate().collect(),
102        }
103    }
104
105    /// Naive division. Division by a SparsePoly (of k elements) takes O(kn) in
106    /// the worst case, division by a DensePoly takes O(n^2) in the worst case.
107    fn naive_div(&self, divisor: &Self) -> Option<(DensePolynomial<F>, DensePolynomial<F>)> {
108        if self.is_zero() {
109            Some((DensePolynomial::zero(), DensePolynomial::zero()))
110        } else if divisor.is_zero() {
111            panic!("Dividing by zero polynomial")
112        } else if self.degree() < divisor.degree() {
113            Some((DensePolynomial::zero(), self.clone().into()))
114        } else {
115            // Now we know that self.degree() >= divisor.degree();
116            let mut quotient = vec![F::zero(); self.degree() - divisor.degree() + 1];
117            let mut remainder: DensePolynomial<F> = self.clone().into();
118            // Can unwrap here because we know self is not zero.
119            let divisor_leading_inv = divisor.leading_coefficient().unwrap().inverse().unwrap();
120            while !remainder.is_zero() && remainder.degree() >= divisor.degree() {
121                let cur_q_coeff = *remainder.coeffs.last().unwrap() * divisor_leading_inv;
122                let cur_q_degree = remainder.degree() - divisor.degree();
123                quotient[cur_q_degree] = cur_q_coeff;
124
125                for (i, div_coeff) in divisor.iter_with_index() {
126                    remainder[cur_q_degree + i] -= &(cur_q_coeff * div_coeff);
127                }
128                while remainder.coeffs.last().map(|c| c.is_zero()) == Some(true) {
129                    remainder.coeffs.pop();
130                }
131            }
132            Some((DensePolynomial::from_coefficients_vec(quotient), remainder))
133        }
134    }
135}
136
137impl<'a, F: FftField> DenseOrSparsePolynomial<'a, F> {
138    const SWITCH_TO_HENSEL_DIV: usize = 1 << 8;
139    /// Divide self by another (sparse or dense) polynomial, and returns the
140    /// quotient and remainder.
141    /// If the divisor is a SparsePolynomial, a naive division algorithm is used.
142    /// If the divisor is a DensePolynomial of degree >= 2^8, a O(M(n)) algorithm
143    /// is used, with M(n) the multiplication complexity.
144    /// Since the field is FFT friendly, then M(n) = O(n log n).
145    pub fn divide_with_q_and_r(
146        &self,
147        divisor: &Self,
148    ) -> Option<(DensePolynomial<F>, DensePolynomial<F>)> {
149        match divisor {
150            // Hensel div is faster than naive for big divisors (degree of around 200)
151            DPolynomial(p) if p.degree() >= Self::SWITCH_TO_HENSEL_DIV => {
152                let q = self.hensel_div(divisor)?;
153                let r = &self.clone().into() - &(&q * &divisor.clone().into());
154
155                Some((q, r))
156            },
157            _ => self.naive_div(divisor),
158        }
159    }
160
161    pub fn divide(&self, divisor: &Self) -> Option<DensePolynomial<F>> {
162        match divisor {
163            // Hensel div is faster than naive for big divisors (degree of around 200)
164            DPolynomial(p) if p.degree() >= Self::SWITCH_TO_HENSEL_DIV => self.hensel_div(divisor),
165            _ => self.naive_div(divisor).map(|(q, _)| q),
166        }
167    }
168
169    // Fast poly division, following "Modern Computer Algebra", 3rd edition, section 9.1, algorithm 9.5
170    fn hensel_div(&self, divisor: &Self) -> Option<DensePolynomial<F>> {
171        if self.is_zero() {
172            Some(DensePolynomial::zero())
173        } else if divisor.is_zero() {
174            panic!("Dividing by zero polynomial")
175        } else if self.degree() < divisor.degree() {
176            Some(DensePolynomial::zero())
177        } else {
178            // Convert to DensePoly to get FFT multiplication
179            let dividend_poly: DensePolynomial<F> = self.clone().into();
180            let divisor_poly: DensePolynomial<F> = divisor.clone().into();
181
182            let reverted_dividend = Self::reverse_coeffs(&dividend_poly, self.degree() + 1);
183            let reverted_divisor = Self::reverse_coeffs(&divisor_poly, divisor.degree() + 1);
184
185            // compute rev(divisor)^-1 mod X^(deg q + 1), with deg q = deg dividend - deg divisor
186            let inversion_degree = self.degree() - divisor.degree() + 1;
187            let reverted_divisor_inverse = Self::inverse_mod(&reverted_divisor, inversion_degree);
188
189            // rev(q) = rev(divisor)^-1 * rev(dividend) mod X^(deg q + 1)
190            let reverted_q = DensePolynomial::from_coefficients_slice(
191                &(&reverted_divisor_inverse * &reverted_dividend).coeffs[..inversion_degree],
192            );
193            let q = DensePolynomial::from_coefficients_slice(
194                &Self::reverse_coeffs(&reverted_q, inversion_degree)
195                    .coeffs
196                    .as_slice(),
197            );
198
199            Some(q)
200        }
201    }
202
203    fn inverse_mod(poly: &DensePolynomial<F>, degree: usize) -> DensePolynomial<F> {
204        if poly.coeffs[0].is_zero() {
205            panic!("Cannot compute inverse of 0");
206        }
207
208        // Inverse mod X
209        let mut l = 1;
210        let mut current_inverse =
211            DensePolynomial::from_coefficients_slice(&[F::one() / poly.coeffs[0]]);
212
213        while l < degree {
214            // Lift inverse of poly mod X^l to inverse of poly mod X^2l
215            // (Also called Newton iteration in the reference book)
216            current_inverse = Self::hensel_lift(poly, &current_inverse, l);
217            l *= 2;
218        }
219        let max_coeff = min(current_inverse.coeffs.len(), degree + 1);
220
221        DensePolynomial::from_coefficients_slice(&current_inverse.coeffs[..max_coeff])
222    }
223
224    // Given inverse such that inverse * poly = 1 mod X^base_degree, return
225    // new_inverse such that new_inverse * poly = 1 mod X^(2*base_degree)
226    // More precisely, decomposing inverse * poly as 1 + c * X^base_degree (c is
227    // easily extracted from the coefficient list), then
228    // new_inverse = inverse - X^base_degree (c * inverse mod X^base_degree)
229    // For performance, poly is reduced mod X^2base_degree, and c mod X^base_degree
230    #[inline]
231    fn hensel_lift(
232        poly: &DensePolynomial<F>,
233        inverse: &DensePolynomial<F>,
234        base_degree: usize,
235    ) -> DensePolynomial<F> {
236        let poly_mod_x2deg = DensePolynomial::from_coefficients_slice(
237            &poly.coeffs[..min(2 * base_degree, poly.coeffs.len())],
238        ); // lowest degree guaranteeing >= base_degree coefficients of c
239        let one_plus_cxdeg = &poly_mod_x2deg * inverse;
240
241        if one_plus_cxdeg.degree() < base_degree {
242            // c = 0 thus we can simply return the inverse immediately
243            return inverse.clone();
244        }
245        let c = DensePolynomial::from_coefficients_slice(
246            &one_plus_cxdeg.coeffs[base_degree..min(2 * base_degree, one_plus_cxdeg.coeffs.len())],
247        ); // we keep c of degree <= deg
248
249        let mut new_inverse_coeffs = vec![F::zero(); 2 * base_degree];
250
251        // First fill the lower half
252        new_inverse_coeffs[..min(base_degree, inverse.coeffs.len())]
253            .clone_from_slice(&inverse.coeffs);
254
255        // Compute -inverse * c mod X^base_degree
256        let above_deg_coeffs: Vec<F> = (inverse * &c)
257            .coeffs
258            .iter()
259            .take(min(base_degree, inverse.degree() + &c.degree() + 1))
260            .map(|&x| -x)
261            .collect();
262
263        // Then fill the upper half
264        new_inverse_coeffs[base_degree..min(2 * base_degree, base_degree + above_deg_coeffs.len())]
265            .clone_from_slice(above_deg_coeffs.as_slice());
266
267        DensePolynomial::from_coefficients_vec(new_inverse_coeffs)
268    }
269
270    fn reverse_coeffs(poly: &DensePolynomial<F>, max_degree: usize) -> DensePolynomial<F> {
271        if poly.coeffs.len() > max_degree {
272            panic!(
273                "Polynomial too big (size {}) to be reverted in size {}",
274                poly.coeffs.len(),
275                max_degree
276            );
277        }
278
279        let mut rev_coeffs = vec![F::zero(); max_degree];
280        rev_coeffs[..poly.coeffs.len()].clone_from_slice(
281            &poly
282                .coeffs
283                .clone()
284                .into_iter()
285                .rev()
286                .collect::<Vec<_>>()
287                .as_slice(),
288        );
289
290        DensePolynomial::from_coefficients_vec(rev_coeffs)
291    }
292}
293
294impl<'a, F: 'a + FftField> DenseOrSparsePolynomial<'a, F> {
295    /// Construct `Evaluations` by evaluating a polynomial over the domain
296    /// `domain`.
297    pub fn evaluate_over_domain<D: EvaluationDomain<F>>(
298        poly: impl Into<Self>,
299        domain: D,
300    ) -> Evaluations<F, D> {
301        let poly = poly.into();
302        poly.eval_over_domain_helper(domain)
303    }
304
305    fn eval_over_domain_helper<D: EvaluationDomain<F>>(self, domain: D) -> Evaluations<F, D> {
306        let eval_sparse_poly = |s: &SparsePolynomial<F>| {
307            let evals = domain.elements().map(|elem| s.evaluate(&elem)).collect();
308            Evaluations::from_vec_and_domain(evals, domain)
309        };
310
311        match self {
312            SPolynomial(Cow::Borrowed(s)) => eval_sparse_poly(s),
313            SPolynomial(Cow::Owned(s)) => eval_sparse_poly(&s),
314            DPolynomial(Cow::Borrowed(d)) => {
315                if d.is_zero() {
316                    Evaluations::zero(domain)
317                } else {
318                    let mut chunks = d.coeffs.chunks(domain.size());
319                    let mut first = chunks.next().unwrap().to_vec();
320                    let offset = domain.coset_offset();
321                    // Reduce the coefficients of the polynomial mod X^domain.size()
322                    for (i, chunk) in chunks.enumerate() {
323                        if offset.is_one() {
324                            cfg_iter_mut!(first).zip(chunk).for_each(|(x, y)| *x += y);
325                        } else {
326                            let offset_power = offset.pow([((i + 1) * domain.size()) as u64]);
327                            cfg_iter_mut!(first)
328                                .zip(chunk)
329                                .for_each(|(x, y)| *x += offset_power * y);
330                        }
331                    }
332                    domain.fft_in_place(&mut first);
333                    Evaluations::from_vec_and_domain(first, domain)
334                }
335            },
336            DPolynomial(Cow::Owned(mut d)) => {
337                if d.is_zero() {
338                    Evaluations::zero(domain)
339                } else {
340                    let mut chunks = d.coeffs.chunks_mut(domain.size());
341                    let first = chunks.next().unwrap();
342                    let offset = domain.coset_offset();
343                    // Reduce the coefficients of the polynomial mod X^domain.size()
344                    for (i, chunk) in chunks.enumerate() {
345                        if offset.is_one() {
346                            cfg_iter_mut!(first).zip(chunk).for_each(|(x, y)| *x += y);
347                        } else {
348                            let offset_power = offset.pow([((i + 1) * domain.size()) as u64]);
349                            cfg_iter_mut!(first)
350                                .zip(chunk)
351                                .for_each(|(x, y)| *x += offset_power * y);
352                        }
353                    }
354                    domain.fft_in_place(&mut d.coeffs);
355                    Evaluations::from_vec_and_domain(d.coeffs, domain)
356                }
357            },
358        }
359    }
360}
ark_poly/polynomial/univariate/mod.rs

ark_poly/polynomial/univariate/
mod.rs
