series/
poly.rs

1use crate::traits::{AsSlice, KaratsubaMul};
2use crate::util::{trim_end, trim_start};
3use crate::{Coeff, IntoIter, Iter, PolynomialIn, PolynomialInParts};
4use crate::{PolynomialSlice, Series};
5
6use std::ops::{
7    Add, AddAssign, Div, DivAssign, Index, Mul, MulAssign, Neg, Range,
8    RangeFrom, RangeFull, RangeInclusive, RangeTo, RangeToInclusive, Sub,
9    SubAssign,
10};
11use std::{convert, iter};
12
13use num_traits::{One, Zero};
14
15/// Laurent polynomial in a single variable
16#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
17#[derive(PartialEq, Eq, Debug, Clone, Hash, Ord, PartialOrd)]
18pub struct Polynomial<C: Coeff> {
19    pub(crate) min_pow: Option<isize>,
20    pub(crate) coeffs: Vec<C>,
21}
22
23impl<C: Coeff> Polynomial<C> {
24    /// Create a new Laurent polynomial
25    ///
26    /// # Example
27    ///
28    /// This creates a Laurent polynomial starting at power -1 with
29    /// coefficients 1, 2, 3. In other words, the polynomial x^-1 + 2 + 3*x.
30    /// ```rust
31    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
32    /// ```
33    pub fn new(min_pow: isize, coeffs: Vec<C>) -> Polynomial<C> {
34        let mut res = Polynomial {
35            min_pow: Some(min_pow),
36            coeffs,
37        };
38        res.trim();
39        res
40    }
41
42    /// Turn into a polynomial in a named variable
43    ///
44    /// # Example
45    ///
46    /// ```rust
47    /// let s = series::Polynomial::new(-1, vec!(1,2,3)).in_var("x");
48    /// assert_eq!(s.var(), &"x");
49    /// ```
50    pub fn in_var<Var>(self, var: Var) -> PolynomialIn<Var, C> {
51        PolynomialIn { var, poly: self }
52    }
53
54    /// Get the leading power of the polynomial variable
55    ///
56    /// For vanishing polynomials `None` is returned
57    ///
58    /// # Example
59    ///
60    /// ```rust
61    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
62    /// assert_eq!(p.min_pow(), Some(-1));
63    /// let p = series::Polynomial::new(-1, vec![0]);
64    /// assert_eq!(p.min_pow(), None);
65    /// ```
66    pub fn min_pow(&self) -> Option<isize> {
67        self.min_pow
68    }
69
70    /// Get the highest power of the polynomial variable
71    ///
72    /// For vanishing polynomials `None` is returned
73    ///
74    /// # Example
75    ///
76    /// ```rust
77    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
78    /// assert_eq!(p.max_pow(), Some(1));
79    /// let p = series::Polynomial::new(-1, vec![0]);
80    /// assert_eq!(p.max_pow(), None);
81    /// ```
82    pub fn max_pow(&self) -> Option<isize> {
83        self.min_pow.map(|c| c - 1 + self.len() as isize)
84    }
85
86    /// Get the difference between the highest and the lowest power of
87    /// the polynomial variable
88    ///
89    /// # Example
90    ///
91    /// ```rust
92    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
93    /// assert_eq!(p.len(), 3);
94    /// ```
95    pub fn len(&self) -> usize {
96        self.coeffs.len()
97    }
98
99    /// Check if the polynomial is zero
100    ///
101    /// # Example
102    ///
103    /// ```rust
104    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
105    /// assert!(!p.is_empty());
106    ///
107    /// let p = series::Polynomial::new(-1, vec!(0));
108    /// assert!(p.is_empty());
109    /// ```
110    pub fn is_empty(&self) -> bool {
111        self.coeffs.is_empty()
112    }
113
114    /// Iterator over the polynomial powers and coefficients.
115    ///
116    /// # Example
117    ///
118    /// ```rust
119    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
120    /// let mut iter = p.iter();
121    /// assert_eq!(iter.next(), Some((-1, &1)));
122    /// assert_eq!(iter.next(), Some((0, &2)));
123    /// assert_eq!(iter.next(), Some((1, &3)));
124    /// assert_eq!(iter.next(), None);
125    /// ```
126    pub fn iter(&self) -> Iter<C> {
127        (self.min_pow().unwrap_or(0)..).zip(self.coeffs.iter())
128    }
129
130    /// Turn a polynomial into a series with the given cutoff
131    ///
132    /// # Example
133    ///
134    /// ```rust
135    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
136    /// let s = series::Series::with_cutoff(-1..5, vec!(1,2,3));
137    /// assert_eq!(p.cutoff_at(5), s);
138    /// ```
139    ///
140    /// # Panics
141    ///
142    /// Panics if the cutoff power is lower than the starting power
143    ///
144    pub fn cutoff_at(self, cutoff_pow: isize) -> Series<C> {
145        Series::with_cutoff(
146            self.min_pow.unwrap_or(cutoff_pow)..cutoff_pow,
147            self.coeffs,
148        )
149    }
150
151    fn as_empty_slice(&self) -> PolynomialSlice<'_, C> {
152        PolynomialSlice::new(0, &self.coeffs[self.len()..])
153    }
154
155    /// Try to get the coefficient of the polynomial variable to the
156    /// given power. Returns [None] if `pow` is less than
157    /// [min_pow](Self::min_pow) or greater than
158    /// [max_pow](Self::max_pow).
159    ///
160    /// # Example
161    ///
162    /// ```rust
163    /// let p = series::Polynomial::new(-1, vec!(1,0,3));
164    /// assert_eq!(p.try_coeff(-5), None);
165    /// assert_eq!(p.try_coeff(-2), None);
166    /// assert_eq!(p.try_coeff(-1), Some(&1));
167    /// assert_eq!(p.try_coeff(0), Some(&0));
168    /// assert_eq!(p.try_coeff(1), Some(&3));
169    /// assert_eq!(p.try_coeff(2), None);
170    /// assert_eq!(p.try_coeff(5), None);
171    /// ```
172    pub fn try_coeff(&self, pow: isize) -> Option<&C> {
173        self.as_slice(..).try_coeff(pow)
174    }
175
176    /// Apply a function to a specific coefficient
177    ///
178    /// `f(c)` is applied to the coefficient `c` of the variable to
179    /// the power `pow`.
180    ///
181    /// # Example
182    ///
183    /// ```rust
184    /// let mut p = series::Polynomial::new(-1, vec![1,2,3]);
185    /// p.apply_at(0, |c| *c = 0);
186    /// assert_eq!(p.coeff(0), &0);
187    ///
188    /// // We can remove existing terms and add new ones
189    /// p.apply_at(-1, |c| *c = 0);
190    /// assert_eq!(p.min_pow(), Some(1));
191    /// p.apply_at(-3, |c| *c = 1);
192    /// assert_eq!(p.min_pow(), Some(-3));
193    /// assert_eq!(p.coeff(-3), &1);
194    /// p.apply_at(10, |c| *c = 1);
195    /// assert_eq!(p.max_pow(), Some(10));
196    /// assert_eq!(p.coeff(10), &1);
197    /// ```
198    pub fn apply_at<F: FnOnce(&mut C)>(&mut self, pow: isize, f: F) {
199        let Some(min_pow) = self.min_pow() else {
200            return self.apply_at_new(pow, f);
201        };
202        if pow < min_pow {
203            return self.apply_at_new_front(pow, f);
204        } else if pow >= min_pow + self.len() as isize {
205            return self.apply_at_new_back(pow, f);
206        }
207        let index = (pow - min_pow) as usize;
208        f(&mut self.coeffs[index]);
209        if (index == 0 || 1 + index == self.len())
210            && self.coeffs[index].is_zero()
211        {
212            self.trim()
213        }
214    }
215
216    fn apply_at_new<F: FnOnce(&mut C)>(&mut self, pow: isize, f: F) {
217        assert!(self.is_empty());
218        let c = gen_from_zero(f);
219        if c.is_zero() {
220            return;
221        };
222        self.min_pow = Some(pow);
223        self.coeffs = vec![c];
224    }
225
226    fn apply_at_new_back<F: FnOnce(&mut C)>(&mut self, pow: isize, f: F) {
227        let c = gen_from_zero(f);
228        if c.is_zero() {
229            return;
230        };
231        let new_len = (pow - self.min_pow().unwrap()) as usize;
232        self.coeffs.resize_with(new_len, || C::zero());
233        self.coeffs.push(c)
234    }
235
236    fn apply_at_new_front<F: FnOnce(&mut C)>(&mut self, pow: isize, f: F) {
237        let c = gen_from_zero(f);
238        if c.is_zero() {
239            return;
240        };
241        let nnew = (self.min_pow().unwrap() - pow) as usize;
242        self.coeffs.push(c);
243        self.coeffs
244            .resize_with(self.len() + (nnew - 1), || C::zero());
245        self.coeffs.rotate_right(nnew);
246        self.min_pow = Some(pow);
247    }
248
249    /// Transform all coefficients
250    ///
251    /// `f(p, c)` is applied to each monomial, where `p` is the power
252    /// of the variable and `c` a mutable reference to the
253    /// coefficient. `p` takes all values in the range
254    /// `min_pow()..=max_pow()`.
255    ///
256    /// # Example
257    ///
258    /// Replace each coefficient by its square
259    /// ```rust
260    /// let mut s = series::Polynomial::new(-1, vec!(1,2,3,4));
261    /// s.for_each(|_, c| *c *= *c);
262    /// assert_eq!(s.coeff(-1), &1);
263    /// assert_eq!(s.coeff(0), &4);
264    /// assert_eq!(s.coeff(1), &9);
265    /// assert_eq!(s.coeff(2), &16);
266    /// ```
267    pub fn for_each<F>(&mut self, mut f: F)
268    where
269        F: FnMut(isize, &mut C),
270    {
271        let Some(min_pow) = self.min_pow else { return };
272        for (n, c) in &mut self.coeffs.iter_mut().enumerate() {
273            f(min_pow + n as isize, c)
274        }
275        self.trim();
276    }
277}
278
279impl<C: 'static + Coeff + Send + Sync> Polynomial<C> {
280    /// Get the coefficient of the polynomial variable to the
281    /// given power.
282    ///
283    /// # Example
284    ///
285    /// ```rust
286    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
287    /// assert_eq!(p.coeff(-5), &0);
288    /// assert_eq!(p.coeff(-2), &0);
289    /// assert_eq!(p.coeff(-1), &1);
290    /// assert_eq!(p.coeff(0), &2);
291    /// assert_eq!(p.coeff(1), &3);
292    /// assert_eq!(p.coeff(2), &0);
293    /// assert_eq!(p.coeff(5), &0);
294    /// ```
295    pub fn coeff(&self, pow: isize) -> &C {
296        self.as_slice(..).coeff(pow)
297    }
298}
299
300fn gen_from_zero<C: Zero, F: FnOnce(&mut C)>(f: F) -> C {
301    let mut c = C::zero();
302    f(&mut c);
303    c
304}
305
306impl<C: Coeff> Default for Polynomial<C> {
307    fn default() -> Self {
308        Self {
309            min_pow: None,
310            coeffs: vec![],
311        }
312    }
313}
314
315impl<'a, C: 'a + Coeff> AsSlice<'a, Range<isize>> for Polynomial<C> {
316    type Output = PolynomialSlice<'a, C>;
317
318    fn as_slice(&'a self, r: Range<isize>) -> Self::Output {
319        if let Some(min_pow) = self.min_pow() {
320            let start = (r.start - min_pow) as usize;
321            let end = (r.end - min_pow) as usize;
322            PolynomialSlice::new(r.start, &self.coeffs[start..end])
323        } else {
324            self.as_empty_slice()
325        }
326    }
327}
328
329impl<'a, C: 'a + Coeff> AsSlice<'a, RangeInclusive<isize>> for Polynomial<C> {
330    type Output = PolynomialSlice<'a, C>;
331
332    fn as_slice(&'a self, r: RangeInclusive<isize>) -> Self::Output {
333        if let Some(min_pow) = self.min_pow() {
334            let (start, end) = r.into_inner();
335            let ustart = (start - min_pow) as usize;
336            let end = (end - min_pow) as usize;
337            PolynomialSlice::new(start, &self.coeffs[ustart..=end])
338        } else {
339            self.as_empty_slice()
340        }
341    }
342}
343
344impl<'a, C: 'a + Coeff> AsSlice<'a, RangeToInclusive<isize>> for Polynomial<C> {
345    type Output = PolynomialSlice<'a, C>;
346
347    fn as_slice(&'a self, r: RangeToInclusive<isize>) -> Self::Output {
348        if let Some(min_pow) = self.min_pow() {
349            let end = (r.end - min_pow) as usize;
350            PolynomialSlice::new(min_pow, &self.coeffs[..=end])
351        } else {
352            self.as_empty_slice()
353        }
354    }
355}
356
357impl<'a, C: 'a + Coeff> AsSlice<'a, RangeFrom<isize>> for Polynomial<C> {
358    type Output = PolynomialSlice<'a, C>;
359
360    fn as_slice(&'a self, r: RangeFrom<isize>) -> Self::Output {
361        if let Some(min_pow) = self.min_pow() {
362            let start = r.start - min_pow;
363            PolynomialSlice::new(r.start, &self.coeffs[(start as usize)..])
364        } else {
365            self.as_empty_slice()
366        }
367    }
368}
369
370impl<'a, C: 'a + Coeff> AsSlice<'a, RangeTo<isize>> for Polynomial<C> {
371    type Output = PolynomialSlice<'a, C>;
372
373    fn as_slice(&'a self, r: RangeTo<isize>) -> Self::Output {
374        if let Some(min_pow) = self.min_pow() {
375            let end = (r.end - min_pow) as usize;
376            PolynomialSlice::new(min_pow, &self.coeffs[..end])
377        } else {
378            self.as_empty_slice()
379        }
380    }
381}
382
383impl<'a, C: 'a + Coeff> AsSlice<'a, RangeFull> for Polynomial<C> {
384    type Output = PolynomialSlice<'a, C>;
385
386    fn as_slice(&'a self, r: RangeFull) -> Self::Output {
387        if let Some(min_pow) = self.min_pow() {
388            PolynomialSlice::new(min_pow, &self.coeffs[r])
389        } else {
390            self.as_empty_slice()
391        }
392    }
393}
394
395impl<C: Coeff> convert::From<Series<C>> for Polynomial<C> {
396    fn from(s: Series<C>) -> Self {
397        Polynomial::new(s.min_pow, s.coeffs)
398    }
399}
400
401impl<C: Coeff> Index<isize> for Polynomial<C> {
402    type Output = C;
403
404    /// Get the coefficient of the polynomial variable to the
405    /// given power.
406    ///
407    /// # Panics
408    ///
409    /// Panics if the index is smaller than the leading power or
410    /// bigger than the highest power
411    ///
412    /// # Example
413    ///
414    /// ```rust
415    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
416    /// assert_eq!(p[-1], 1);
417    /// assert_eq!(p[0], 2);
418    /// assert_eq!(p[1], 3);
419    /// assert!(std::panic::catch_unwind(|| p[-2]).is_err());
420    /// assert!(std::panic::catch_unwind(|| p[2]).is_err());
421    /// ```
422    fn index(&self, index: isize) -> &Self::Output {
423        &self.coeffs[(index - self.min_pow.unwrap()) as usize]
424    }
425}
426
427impl<C: Coeff> std::iter::IntoIterator for Polynomial<C> {
428    type Item = (isize, C);
429    type IntoIter = crate::IntoIter<C>;
430
431    /// Consuming iterator over the polynomial powers and coefficients.
432    ///
433    /// # Example
434    ///
435    /// ```rust
436    /// let p = series::Polynomial::new(-1, vec!(1,2,3));
437    /// let mut iter = p.into_iter();
438    /// assert_eq!(iter.next(), Some((-1, 1)));
439    /// assert_eq!(iter.next(), Some((0, 2)));
440    /// assert_eq!(iter.next(), Some((1, 3)));
441    /// assert_eq!(iter.next(), None);
442    /// ```
443    fn into_iter(self) -> IntoIter<C> {
444        (self.min_pow().unwrap_or(0)..).zip(self.coeffs)
445    }
446}
447
448impl<C: Coeff> Polynomial<C> {
449    fn trim(&mut self) {
450        let zero = C::zero();
451        trim_end(&mut self.coeffs, &zero);
452        if self.coeffs.is_empty() {
453            self.min_pow = None;
454        } else {
455            let min_pow_shift = trim_start(&mut self.coeffs, &zero);
456            if let Some(min_pow) = self.min_pow.as_mut() {
457                *min_pow += min_pow_shift as isize
458            }
459        }
460    }
461
462    fn extend_min(&mut self, extend: usize) {
463        debug_assert!(self.min_pow.is_some());
464        let to_insert = iter::repeat_with(C::zero).take(extend);
465        self.coeffs.splice(0..0, to_insert);
466        if let Some(min_pow) = self.min_pow.as_mut() {
467            *min_pow -= extend as isize
468        }
469    }
470
471    fn extend_max(&mut self, extend: usize) {
472        let to_insert = iter::repeat_with(C::zero).take(extend);
473        self.coeffs.extend(to_insert);
474    }
475}
476
477impl<C: Coeff + Neg<Output = C>> Neg for Polynomial<C> {
478    type Output = Polynomial<C>;
479
480    /// Compute -p for a Laurent polynomial p
481    ///
482    /// # Example
483    ///
484    /// ```rust
485    /// let p = series::Polynomial::new(-3, vec!(1.,0.,-3.));
486    /// let minus_p = series::Polynomial::new(-3, vec!(-1.,0.,3.));
487    /// assert_eq!(-p, minus_p);
488    /// ```
489    fn neg(self) -> Self::Output {
490        let neg_coeff = self.coeffs.into_iter().map(|c| -c).collect();
491        Polynomial::new(self.min_pow.unwrap_or(0), neg_coeff)
492    }
493}
494
495impl<'a, C: Coeff> Neg for &'a Polynomial<C>
496where
497    PolynomialSlice<'a, C>: Neg,
498{
499    type Output = <PolynomialSlice<'a, C> as Neg>::Output;
500
501    /// Compute -p for a Laurent polynomial p
502    ///
503    /// # Example
504    ///
505    /// ```rust
506    /// let p = series::Polynomial::new(-3, vec!(1.,0.,-3.));
507    /// let minus_p = series::Polynomial::new(-3, vec!(-1.,0.,3.));
508    /// assert_eq!(-&p, minus_p);
509    /// ```
510    fn neg(self) -> Self::Output {
511        self.as_slice(..).neg()
512    }
513}
514
515impl<'a, C: Coeff + Clone> AddAssign<&'a Polynomial<C>> for Polynomial<C>
516where
517    for<'c> C: AddAssign<&'c C>,
518{
519    /// Set p = p + q for two Laurent polynomials p and q
520    ///
521    /// # Example
522    ///
523    /// ```rust
524    /// use series::Polynomial;
525    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
526    /// let q = Polynomial::new(-1, vec!(3., 4., 5.));
527    /// let res = Polynomial::new(-3, vec!(1.,0.,0.,4.,5.));
528    /// p += &q;
529    /// assert_eq!(res, p);
530    /// ```
531    ///
532    fn add_assign(&mut self, other: &'a Polynomial<C>) {
533        self.add_assign(other.as_slice(..))
534    }
535}
536
537impl<'a, C: Coeff + Clone> AddAssign<PolynomialSlice<'a, C>> for Polynomial<C>
538where
539    for<'c> C: AddAssign<&'c C>,
540{
541    fn add_assign(&mut self, other: PolynomialSlice<'a, C>) {
542        if other.min_pow().is_none() {
543            return;
544        }
545        if self.min_pow().is_none() {
546            self.coeffs = other.coeffs.to_owned();
547            self.min_pow = other.min_pow();
548            return;
549        }
550        let min_pow = self.min_pow().unwrap();
551        let other_min_pow = other.min_pow().unwrap();
552        if other_min_pow < min_pow {
553            self.extend_min((min_pow - other_min_pow) as usize);
554        }
555        let max_pow = self.max_pow().unwrap();
556        let other_max_pow = other.max_pow().unwrap();
557        if other_max_pow > max_pow {
558            self.extend_max((other_max_pow - max_pow) as usize);
559        }
560        debug_assert!(self.min_pow().unwrap() <= other_min_pow);
561        debug_assert!(self.max_pow().unwrap() >= other_max_pow);
562        let min_pow = self.min_pow().unwrap();
563        for (pow, coeff) in other.iter() {
564            self.coeffs[(pow - min_pow) as usize] += coeff;
565        }
566        self.trim();
567    }
568}
569
570impl<C: Coeff> AddAssign<Polynomial<C>> for Polynomial<C>
571where
572    for<'c> C: AddAssign<&'c C>,
573    C: Clone + AddAssign,
574{
575    /// Set p = p + q for two polynomial p and q
576    ///
577    /// # Example
578    ///
579    /// ```rust
580    /// use series::Polynomial;
581    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
582    /// let q = Polynomial::new(-1, vec!(3., 4., 5.));
583    /// let res = Polynomial::new(-3, vec!(1.,0.,0.,4.,5.));
584    /// p += q;
585    /// assert_eq!(res, p);
586    /// ```
587    fn add_assign(&mut self, other: Polynomial<C>) {
588        //TODO: code duplication with AddAssign<PolynomialSlice>
589        if other.min_pow().is_none() {
590            return;
591        }
592        if self.min_pow().is_none() {
593            self.coeffs = other.coeffs.to_owned();
594            self.min_pow = other.min_pow();
595            return;
596        }
597        let min_pow = self.min_pow().unwrap();
598        let other_min_pow = other.min_pow().unwrap();
599        if other_min_pow < min_pow {
600            self.extend_min((min_pow - other_min_pow) as usize);
601        }
602        let max_pow = self.max_pow().unwrap();
603        let other_max_pow = other.max_pow().unwrap();
604        if other_max_pow > max_pow {
605            self.extend_max((other_max_pow - max_pow) as usize);
606        }
607        debug_assert!(self.min_pow() <= other.min_pow());
608        debug_assert!(self.max_pow() >= other.max_pow());
609        let min_pow = self.min_pow().unwrap();
610        for (pow, coeff) in other.into_iter() {
611            self.coeffs[(pow - min_pow) as usize] += coeff;
612        }
613        self.trim();
614    }
615}
616
617impl<C: Coeff + Clone, Rhs> Add<Rhs> for Polynomial<C>
618where
619    Polynomial<C>: AddAssign<Rhs>,
620{
621    type Output = Polynomial<C>;
622
623    fn add(mut self, other: Rhs) -> Self::Output {
624        self += other;
625        self
626    }
627}
628
629impl<'a, C: Coeff> SubAssign<&'a Polynomial<C>> for Polynomial<C>
630where
631    for<'c> &'c Polynomial<C>: Neg<Output = Polynomial<C>>,
632    Polynomial<C>: AddAssign<Polynomial<C>>,
633{
634    /// Set p = p - q for two polynomial p and q
635    ///
636    /// # Example
637    ///
638    /// ```rust
639    /// use series::Polynomial;
640    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
641    /// let res = Polynomial::new(0, vec!());
642    /// p -= &p.clone();
643    /// assert_eq!(res, p);
644    /// ```
645    fn sub_assign(&mut self, other: &'a Polynomial<C>) {
646        *self += -other;
647    }
648}
649
650impl<'a, C: Coeff> SubAssign<PolynomialSlice<'a, C>> for Polynomial<C>
651where
652    for<'c> PolynomialSlice<'c, C>: Neg<Output = Polynomial<C>>,
653    Polynomial<C>: AddAssign<Polynomial<C>>,
654{
655    fn sub_assign(&mut self, other: PolynomialSlice<'a, C>) {
656        *self += -other;
657    }
658}
659
660impl<C: Coeff> SubAssign<Polynomial<C>> for Polynomial<C>
661where
662    Polynomial<C>: AddAssign + Neg<Output = Polynomial<C>>,
663{
664    /// Set p = p - q for two polynomial p and q
665    ///
666    /// # Example
667    ///
668    /// ```rust
669    /// use series::Polynomial;
670    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
671    /// let res = Polynomial::new(0, vec!());
672    /// p -= p.clone();
673    /// assert_eq!(res, p);
674    /// ```
675    fn sub_assign(&mut self, other: Polynomial<C>) {
676        *self += -other;
677    }
678}
679
680impl<C: Coeff, T> Sub<T> for Polynomial<C>
681where
682    Polynomial<C>: SubAssign<T>,
683{
684    type Output = Polynomial<C>;
685
686    fn sub(mut self, other: T) -> Self::Output {
687        self -= other;
688        self
689    }
690}
691
692impl<'a, C: Coeff + Clone + AddAssign> MulAssign<&'a Polynomial<C>>
693    for Polynomial<C>
694where
695    Polynomial<C>: MulAssign<PolynomialSlice<'a, C>>,
696{
697    /// Set p = p * q for two polynomials p,q
698    ///
699    /// # Example
700    ///
701    /// ```rust
702    /// use series::Polynomial;
703    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
704    /// p *= &p.clone();
705    /// let res = Polynomial::new(-6, vec!(1.,0.,-6.,0.,9.));
706    /// assert_eq!(res, p);
707    /// ```
708    fn mul_assign(&mut self, other: &'a Polynomial<C>) {
709        self.mul_assign(other.as_slice(..))
710    }
711}
712
713impl<'a, C: Coeff> MulAssign<PolynomialSlice<'a, C>> for Polynomial<C>
714where
715    for<'b> PolynomialSlice<'b, C>:
716        Mul<PolynomialSlice<'a, C>, Output = Polynomial<C>>,
717{
718    fn mul_assign(&mut self, other: PolynomialSlice<'a, C>) {
719        let prod = self.as_slice(..) * other;
720        *self = prod;
721    }
722}
723
724impl<C: Coeff> MulAssign for Polynomial<C>
725where
726    for<'a> Polynomial<C>: MulAssign<&'a Polynomial<C>>,
727{
728    /// Set p = p * q for two polynomials p,q
729    ///
730    /// # Example
731    ///
732    /// ```rust
733    /// use series::Polynomial;
734    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
735    /// p *= &p.clone();
736    /// let res = Polynomial::new(-6, vec!(1.,0.,-6.,0.,9.));
737    /// assert_eq!(res, p);
738    /// ```
739    fn mul_assign(&mut self, other: Polynomial<C>) {
740        *self *= &other
741    }
742}
743
744impl<C: Coeff> MulAssign<C> for Polynomial<C>
745where
746    for<'a> C: MulAssign<&'a C>,
747{
748    /// Multiply each monomial by a factor
749    ///
750    /// # Example
751    ///
752    /// ```rust
753    /// use series::Polynomial;
754    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
755    /// p *= 2.;
756    /// let res = Polynomial::new(-3, vec!(2.,0.,-6.));
757    /// assert_eq!(res, p);
758    /// ```
759    fn mul_assign(&mut self, other: C) {
760        self.mul_assign(&other)
761    }
762}
763
764impl<'a, C: Coeff> MulAssign<&'a C> for Polynomial<C>
765where
766    C: MulAssign<&'a C>,
767{
768    /// Multiply each monomial by a factor
769    ///
770    /// # Example
771    ///
772    /// ```rust
773    /// use series::Polynomial;
774    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
775    /// p *= &2.;
776    /// let res = Polynomial::new(-3, vec!(2.,0.,-6.));
777    /// assert_eq!(res, p);
778    /// ```
779    fn mul_assign(&mut self, other: &'a C) {
780        for c in &mut self.coeffs {
781            *c *= other
782        }
783    }
784}
785
786impl<C: Coeff> DivAssign<C> for Polynomial<C>
787where
788    for<'a> C: DivAssign<&'a C>,
789{
790    /// Divide each monomial by a factor
791    ///
792    /// # Example
793    ///
794    /// ```rust
795    /// use series::Polynomial;
796    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
797    /// p /= 2.;
798    /// let res = Polynomial::new(-3, vec!(0.5,0.,-1.5));
799    /// assert_eq!(res, p);
800    /// ```
801    fn div_assign(&mut self, other: C) {
802        self.div_assign(&other)
803    }
804}
805
806impl<'a, C: Coeff> DivAssign<&'a C> for Polynomial<C>
807where
808    C: DivAssign<&'a C>,
809{
810    /// Divide each monomial by a factor
811    ///
812    /// # Example
813    ///
814    /// ```rust
815    /// use series::Polynomial;
816    /// let mut p = Polynomial::new(-3, vec!(1.,0.,-3.));
817    /// p /= &2.;
818    /// let res = Polynomial::new(-3, vec!(0.5,0.,-1.5));
819    /// assert_eq!(res, p);
820    /// ```
821    fn div_assign(&mut self, other: &'a C) {
822        for c in &mut self.coeffs {
823            *c /= other
824        }
825    }
826}
827
828impl<C: Coeff> AddAssign<C> for Polynomial<C>
829where
830    for<'c> Polynomial<C>: AddAssign<&'c C>,
831{
832    fn add_assign(&mut self, other: C) {
833        self.add_assign(&other)
834    }
835}
836
837impl<'a, C: Coeff> AddAssign<&'a C> for Polynomial<C>
838where
839    C: Clone + AddAssign<&'a C>,
840{
841    fn add_assign(&mut self, other: &'a C) {
842        match self.min_pow() {
843            None => {
844                self.min_pow = Some(0);
845                self.coeffs.push(other.clone())
846            }
847            Some(min_pow) => {
848                if min_pow > 0 {
849                    self.extend_min(min_pow as usize);
850                } else if self.max_pow().unwrap() < 0 {
851                    self.extend_max((-self.max_pow().unwrap()) as usize);
852                }
853                debug_assert!(self.min_pow().unwrap() <= 0);
854                debug_assert!(0 <= self.max_pow().unwrap());
855                let min_pow = self.min_pow().unwrap();
856                self.coeffs[-min_pow as usize] += other;
857            }
858        }
859        self.trim();
860    }
861}
862
863impl<C: Coeff> SubAssign<C> for Polynomial<C>
864where
865    Polynomial<C>: AddAssign<C>,
866    C: Neg<Output = C> + AddAssign,
867{
868    fn sub_assign(&mut self, other: C) {
869        self.add_assign(-other);
870    }
871}
872
873impl<'c, C: Coeff> SubAssign<&'c C> for Polynomial<C>
874where
875    Polynomial<C>: AddAssign<C>,
876    C: AddAssign,
877    &'c C: Neg<Output = C>,
878{
879    fn sub_assign(&mut self, other: &'c C) {
880        self.add_assign(-other);
881    }
882}
883
884impl<C: Coeff> Mul for Polynomial<C>
885where
886    Polynomial<C>: MulAssign,
887{
888    type Output = Polynomial<C>;
889
890    fn mul(mut self, other: Polynomial<C>) -> Self::Output {
891        self *= other;
892        self
893    }
894}
895
896impl<'a, C: Coeff> Mul<&'a Polynomial<C>> for Polynomial<C>
897where
898    Polynomial<C>: MulAssign<PolynomialSlice<'a, C>>,
899{
900    type Output = Polynomial<C>;
901
902    fn mul(self, other: &'a Polynomial<C>) -> Self::Output {
903        self * other.as_slice(..)
904    }
905}
906
907impl<'a, C: Coeff> Mul<PolynomialSlice<'a, C>> for Polynomial<C>
908where
909    Polynomial<C>: MulAssign<PolynomialSlice<'a, C>>,
910{
911    type Output = Polynomial<C>;
912
913    fn mul(mut self, other: PolynomialSlice<'a, C>) -> Self::Output {
914        self *= other;
915        self
916    }
917}
918
919impl<C: Coeff> Mul<C> for Polynomial<C>
920where
921    for<'c> C: MulAssign<&'c C>,
922{
923    type Output = Polynomial<C>;
924
925    fn mul(mut self, other: C) -> Self::Output {
926        self *= &other;
927        self
928    }
929}
930
931impl<'a, C: Coeff> Mul<&'a C> for Polynomial<C>
932where
933    for<'c> C: MulAssign<&'c C>,
934{
935    type Output = Polynomial<C>;
936
937    fn mul(mut self, other: &'a C) -> Self::Output {
938        self *= other;
939        self
940    }
941}
942
943impl<C: Coeff> Div<C> for Polynomial<C>
944where
945    for<'c> C: DivAssign<&'c C>,
946{
947    type Output = Polynomial<C>;
948
949    fn div(mut self, other: C) -> Self::Output {
950        self /= &other;
951        self
952    }
953}
954
955impl<'a, C: Coeff> Div<&'a C> for Polynomial<C>
956where
957    for<'c> C: DivAssign<&'c C>,
958{
959    type Output = Polynomial<C>;
960
961    fn div(mut self, other: &'a C) -> Self::Output {
962        self /= other;
963        self
964    }
965}
966
967impl<'a, C: Coeff, T> Mul<T> for &'a Polynomial<C>
968where
969    PolynomialSlice<'a, C>: Mul<T, Output = Polynomial<C>>,
970{
971    type Output = Polynomial<C>;
972
973    fn mul(self, other: T) -> Self::Output {
974        self.as_slice(..) * other
975    }
976}
977
978impl<'a, C: Coeff, T> Div<T> for &'a Polynomial<C>
979where
980    PolynomialSlice<'a, C>: Div<T, Output = Polynomial<C>>,
981{
982    type Output = Polynomial<C>;
983
984    fn div(self, other: T) -> Self::Output {
985        self.as_slice(..) / other
986    }
987}
988
989impl<'a, C: Coeff, T> Add<T> for &'a Polynomial<C>
990where
991    PolynomialSlice<'a, C>: Add<T, Output = Polynomial<C>>,
992{
993    type Output = Polynomial<C>;
994
995    fn add(self, other: T) -> Self::Output {
996        self.as_slice(..) + other
997    }
998}
999
1000impl<'a, C: Coeff, T> Sub<T> for &'a Polynomial<C>
1001where
1002    PolynomialSlice<'a, C>: Sub<T, Output = Polynomial<C>>,
1003{
1004    type Output = Polynomial<C>;
1005
1006    fn sub(self, other: T) -> Self::Output {
1007        self.as_slice(..) - other
1008    }
1009}
1010
1011impl<'a, C: Coeff + Clone> From<PolynomialSlice<'a, C>> for Polynomial<C> {
1012    fn from(s: PolynomialSlice<'a, C>) -> Polynomial<C> {
1013        Polynomial::new(s.min_pow.unwrap_or(0), s.coeffs.to_vec())
1014    }
1015}
1016
1017impl<C: Coeff> Zero for Polynomial<C>
1018where
1019    Polynomial<C>: Add<Output = Polynomial<C>>,
1020{
1021    fn zero() -> Self {
1022        Self {
1023            min_pow: None,
1024            coeffs: vec![],
1025        }
1026    }
1027
1028    fn is_zero(&self) -> bool {
1029        self.coeffs.is_empty()
1030    }
1031}
1032
1033impl<C: Coeff> One for Polynomial<C>
1034where
1035    Polynomial<C>: Mul<Output = Polynomial<C>>,
1036{
1037    fn one() -> Self {
1038        Self {
1039            min_pow: Some(0),
1040            coeffs: vec![C::one()],
1041        }
1042    }
1043
1044    fn is_one(&self) -> bool {
1045        self.min_pow == Some(0)
1046            && self.coeffs.len() == 1
1047            && self.coeffs[0].is_one()
1048    }
1049}
1050
1051impl<C: Coeff, Var> From<PolynomialIn<Var, C>> for Polynomial<C> {
1052    fn from(source: PolynomialIn<Var, C>) -> Self {
1053        let PolynomialInParts {
1054            var: _,
1055            min_pow,
1056            coeffs,
1057        } = source.into();
1058        Self { min_pow, coeffs }
1059    }
1060}
1061
1062/// Data parts of a polynomial
1063///
1064/// # Example
1065///
1066/// ```rust
1067/// // destructure a polynomial
1068/// let p = series::Polynomial::new(-1, vec![1,2,3]);
1069/// let series::PolynomialParts{min_pow, coeffs} = p.into();
1070/// assert_eq!(min_pow, Some(-1));
1071/// assert_eq!(coeffs, vec![1,2,3]);
1072/// ```
1073#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
1074#[derive(PartialEq, Eq, Debug, Clone, Hash, Ord, PartialOrd)]
1075pub struct PolynomialParts<C> {
1076    pub min_pow: Option<isize>,
1077    pub coeffs: Vec<C>,
1078}
1079
1080impl<C: Coeff> From<Polynomial<C>> for PolynomialParts<C> {
1081    fn from(p: Polynomial<C>) -> Self {
1082        PolynomialParts {
1083            min_pow: p.min_pow,
1084            coeffs: p.coeffs,
1085        }
1086    }
1087}
1088
1089impl<'a, 'b, C: Coeff> KaratsubaMul<&'b Polynomial<C>> for &'a Polynomial<C>
1090where
1091    C: Clone,
1092    for<'c> C: AddAssign,
1093    for<'c> Polynomial<C>:
1094        AddAssign<&'c Polynomial<C>> + SubAssign<&'c Polynomial<C>>,
1095    Polynomial<C>: AddAssign<Polynomial<C>> + SubAssign<Polynomial<C>>,
1096    for<'c> PolynomialSlice<'c, C>: Add<Output = Polynomial<C>>,
1097    for<'c> &'c C: Mul<Output = C>,
1098{
1099    type Output = Polynomial<C>;
1100
1101    fn karatsuba_mul(
1102        self,
1103        rhs: &'b Polynomial<C>,
1104        min_size: usize,
1105    ) -> Self::Output {
1106        self.as_slice(..).karatsuba_mul(rhs.as_slice(..), min_size)
1107    }
1108}
1109
1110impl<'a, 'b, C: Coeff> KaratsubaMul<PolynomialSlice<'b, C>>
1111    for &'a Polynomial<C>
1112where
1113    C: Clone,
1114    for<'c> C: AddAssign,
1115    for<'c> Polynomial<C>:
1116        AddAssign<&'c Polynomial<C>> + SubAssign<&'c Polynomial<C>>,
1117    Polynomial<C>: AddAssign<Polynomial<C>> + SubAssign<Polynomial<C>>,
1118    for<'c> PolynomialSlice<'c, C>: Add<Output = Polynomial<C>>,
1119    for<'c> &'c C: Mul<Output = C>,
1120{
1121    type Output = Polynomial<C>;
1122
1123    fn karatsuba_mul(
1124        self,
1125        rhs: PolynomialSlice<'b, C>,
1126        min_size: usize,
1127    ) -> Self::Output {
1128        self.as_slice(..).karatsuba_mul(rhs, min_size)
1129    }
1130}
1131
1132// dubious helpers trait that only serve to prevent obscure
1133// compiler errors (rust 1.36.0)
1134pub(crate) trait MulHelper<'a, 'b, C: Coeff> {
1135    fn add_prod(
1136        &mut self,
1137        a: PolynomialSlice<'a, C>,
1138        b: PolynomialSlice<'b, C>,
1139        min_karatsuba_size: usize,
1140    );
1141
1142    fn add_prod_naive(
1143        &mut self,
1144        a: PolynomialSlice<'a, C>,
1145        b: PolynomialSlice<'b, C>,
1146    );
1147
1148    fn add_prod_karatsuba(
1149        &mut self,
1150        a: PolynomialSlice<'a, C>,
1151        b: PolynomialSlice<'b, C>,
1152        min_karatsuba_size: usize,
1153    );
1154
1155    fn add_prod_unchecked(
1156        &mut self,
1157        a: PolynomialSlice<'a, C>,
1158        b: PolynomialSlice<'b, C>,
1159        min_karatsuba_size: usize,
1160    );
1161
1162    fn resize_to_fit(
1163        &mut self,
1164        a: PolynomialSlice<'a, C>,
1165        b: PolynomialSlice<'b, C>,
1166    );
1167}
1168
1169impl<'a, 'b, C: Coeff> MulHelper<'a, 'b, C> for Polynomial<C>
1170where
1171    C: 'a + 'b + Clone,
1172    for<'c> C: AddAssign,
1173    for<'c> Polynomial<C>:
1174        AddAssign<&'c Polynomial<C>> + SubAssign<&'c Polynomial<C>>,
1175    Polynomial<C>: AddAssign<Polynomial<C>> + SubAssign<Polynomial<C>>,
1176    for<'c> PolynomialSlice<'c, C>:
1177        Add<Output = Polynomial<C>> + Mul<Output = Polynomial<C>>,
1178    for<'c> &'c C: Mul<Output = C>,
1179{
1180    fn add_prod(
1181        &mut self,
1182        a: PolynomialSlice<'a, C>,
1183        b: PolynomialSlice<'b, C>,
1184        min_karatsuba_size: usize,
1185    ) {
1186        if a.min_pow().is_none() || b.min_pow().is_none() {
1187            return;
1188        }
1189        self.resize_to_fit(a, b);
1190        self.add_prod_unchecked(a, b, min_karatsuba_size);
1191        self.trim();
1192    }
1193
1194    fn resize_to_fit(
1195        &mut self,
1196        a: PolynomialSlice<'a, C>,
1197        b: PolynomialSlice<'b, C>,
1198    ) {
1199        debug_assert_ne!(a.min_pow(), None);
1200        debug_assert_ne!(b.min_pow(), None);
1201        let prod_min_pow = a.min_pow().unwrap() + b.min_pow().unwrap();
1202        match self.min_pow() {
1203            Some(min_pow) => {
1204                if min_pow > prod_min_pow {
1205                    let num_missing = (min_pow - prod_min_pow) as usize;
1206                    let to_prepend = iter::repeat(C::zero()).take(num_missing);
1207                    self.coeffs.splice(0..0, to_prepend);
1208                    self.min_pow = Some(prod_min_pow);
1209                }
1210            }
1211            None => {
1212                self.min_pow = Some(prod_min_pow);
1213            }
1214        }
1215        let prod_len = a.len() + b.len() - 1;
1216        if self.len() < prod_len {
1217            self.coeffs.resize(prod_len, C::zero());
1218        }
1219    }
1220
1221    fn add_prod_unchecked(
1222        &mut self,
1223        a: PolynomialSlice<'a, C>,
1224        b: PolynomialSlice<'b, C>,
1225        min_karatsuba_size: usize,
1226    ) {
1227        if std::cmp::min(a.len(), b.len()) < min_karatsuba_size {
1228            self.add_prod_naive(a, b);
1229        } else {
1230            // TODO: split a or b if it's too long?
1231            self.add_prod_karatsuba(a, b, min_karatsuba_size);
1232        }
1233    }
1234
1235    fn add_prod_karatsuba(
1236        &mut self,
1237        a: PolynomialSlice<'a, C>,
1238        b: PolynomialSlice<'b, C>,
1239        min_karatsuba_size: usize,
1240    ) {
1241        let mid = ((std::cmp::min(a.len(), b.len()) + 1) / 2) as isize;
1242        let (a_low, mut a_high) = a.split_at(a.min_pow().unwrap() + mid);
1243        let (b_low, mut b_high) = b.split_at(b.min_pow().unwrap() + mid);
1244        if let Some(min_pow) = a_high.min_pow.as_mut() {
1245            *min_pow -= mid
1246        }
1247        if let Some(min_pow) = b_high.min_pow.as_mut() {
1248            *min_pow -= mid
1249        }
1250        let t = a_low + a_high;
1251        let mut u = b_low + b_high;
1252        if let Some(min_pow) = u.min_pow.as_mut() {
1253            *min_pow += mid
1254        }
1255        self.add_prod_unchecked(
1256            t.as_slice(..),
1257            u.as_slice(..),
1258            min_karatsuba_size,
1259        );
1260        let mut t = a_low * b_low;
1261        *self += &t;
1262        if let Some(min_pow) = t.min_pow.as_mut() {
1263            *min_pow += mid
1264        }
1265        *self -= t;
1266        let mut t = a_high * b_high;
1267        if let Some(min_pow) = t.min_pow.as_mut() {
1268            *min_pow += mid
1269        }
1270        *self -= &t;
1271        if let Some(min_pow) = t.min_pow.as_mut() {
1272            *min_pow += mid
1273        }
1274        *self += t;
1275    }
1276
1277    fn add_prod_naive(
1278        &mut self,
1279        a: PolynomialSlice<'a, C>,
1280        b: PolynomialSlice<'b, C>,
1281    ) {
1282        let min_pow = self.min_pow().unwrap();
1283        for (i, a) in a.iter() {
1284            for (j, b) in b.iter() {
1285                self.coeffs[(i + j - min_pow) as usize] += a * b;
1286            }
1287        }
1288    }
1289}
series/poly.rs

series/
poly.rs
