rustnomial/polynomial/
polynomial.rs

1use alloc::vec::Vec;
2use core::ops::{
3    Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Shl, ShlAssign, Shr,
4    ShrAssign, Sub, SubAssign,
5};
6
7use num::{One, Zero};
8
9use crate::numerics::{Abs, CanNegate, IsNegativeOne, IsPositive, TryFromUsizeContinuous};
10use crate::polynomial::find_roots::{find_roots, Roots};
11use crate::{
12    Degree, Derivable, Evaluable, FreeSizePolynomial, Integrable, Integral, MutablePolynomial,
13    SizedPolynomial, Term, TryAddError,
14};
15
16#[macro_export]
17macro_rules! polynomial {
18    ( $( $x:expr ),* ) => {
19        {
20            use $crate::Polynomial;
21            Polynomial::new(vec![$($x,)*])
22        }
23    };
24}
25
26#[derive(Debug, Clone)]
27/// A type that stores terms of a polynomial in a Vec.
28pub struct Polynomial<N> {
29    pub terms: Vec<N>,
30}
31
32pub(crate) fn first_nonzero_index<N>(coeffs: &[N]) -> usize
33where
34    N: Zero + Copy,
35{
36    for (degree, chunk) in coeffs.chunks_exact(4).enumerate() {
37        for (index, &val) in chunk.iter().enumerate() {
38            if !val.is_zero() {
39                return degree * 4 + index;
40            }
41        }
42    }
43
44    let mut len = coeffs.chunks_exact(4).len() * 4;
45    for &value in coeffs.chunks_exact(4).remainder().iter() {
46        if !value.is_zero() {
47            return len;
48        }
49        len += 1;
50    }
51
52    len
53}
54
55fn slice_mul<N>(lhs: &[N], rhs: &[N]) -> Vec<N>
56where
57    N: Mul<Output = N> + AddAssign + Copy + Zero,
58{
59    let rhs = &rhs[first_nonzero_index(&rhs)..];
60    let lhs = &lhs[first_nonzero_index(&lhs)..];
61    let mut terms = vec![N::zero(); rhs.len() + lhs.len() - 1];
62    for (index, &rterm) in rhs.iter().enumerate() {
63        if rterm.is_zero() {
64            continue;
65        }
66        for (&lterm, term) in lhs.iter().zip(terms[index..].iter_mut()) {
67            *term += rterm * lterm;
68        }
69    }
70    terms
71}
72
73fn vec_sub_w_scale<N>(lhs: &mut [N], lhs_degree: usize, rhs: &[N], rhs_deg: usize, rhs_scale: N)
74where
75    N: Copy + Mul<Output = N> + SubAssign,
76{
77    let loc = lhs.len() - lhs_degree - 1;
78    for (lhs_t, rhs_t) in lhs[loc..]
79        .iter_mut()
80        .zip(rhs[rhs.len() - rhs_deg - 1..].iter())
81    {
82        *lhs_t -= (*rhs_t) * rhs_scale;
83    }
84}
85
86pub(crate) fn degree<N>(coeffs: &[N]) -> Degree
87where
88    N: Zero + Copy,
89{
90    let index = first_nonzero_index(coeffs);
91    if index == coeffs.len() {
92        Degree::NegInf
93    } else {
94        Degree::Num(coeffs.len() - index - 1)
95    }
96}
97
98pub(crate) fn first_term<N>(poly_vec: &[N]) -> Term<N>
99where
100    N: Zero + Copy,
101{
102    for (degree, chunk) in poly_vec.chunks_exact(4).enumerate() {
103        for (index, &value) in chunk.iter().enumerate() {
104            if !value.is_zero() {
105                return Term::Term(value, poly_vec.len() - degree * 4 - index - 1);
106            }
107        }
108    }
109
110    let mut index = poly_vec.chunks_exact(4).len() * 4;
111    for &value in poly_vec.chunks_exact(4).remainder().iter() {
112        if !value.is_zero() {
113            return Term::Term(value, poly_vec.len() - index - 1);
114        }
115        index += 1;
116    }
117
118    Term::ZeroTerm
119}
120
121pub(crate) fn term_with_deg<N: Zero + Copy>(terms: &[N], degree: usize) -> Term<N> {
122    if degree < terms.len() {
123        Term::new(terms[terms.len() - degree - 1], degree)
124    } else {
125        Term::ZeroTerm
126    }
127}
128
129impl<N> Polynomial<N>
130where
131    N: Zero + Copy,
132{
133    /// Returns a `Polynomial` with the corresponding terms,
134    /// in order of ax^n + bx^(n-1) + ... + cx + d
135    ///
136    /// # Arguments
137    ///
138    /// * ` terms ` - A vector of constants, in decreasing order of degree.
139    ///
140    /// # Example
141    ///
142    /// ```
143    /// use rustnomial::Polynomial;
144    /// // Corresponds to 1.0x^2 + 4.0x + 4.0
145    /// let polynomial = Polynomial::new(vec![1.0, 4.0, 4.0]);
146    /// ```
147    pub fn new(terms: Vec<N>) -> Polynomial<N> {
148        let mut p = Polynomial { terms };
149        p.trim();
150        p
151    }
152
153    /// Reduces the size of the `Polynomial` in memory if the leading terms are zero.
154    ///
155    /// # Example
156    ///
157    /// ```
158    /// use rustnomial::Polynomial;
159    /// let mut polynomial = Polynomial::new(vec![1.0, 4.0, 4.0]);
160    /// polynomial.terms = vec![0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 4.0];
161    /// polynomial.trim();
162    /// assert_eq!(vec![1.0, 4.0, 4.0], polynomial.terms);
163    /// ```
164    pub fn trim(&mut self) {
165        let ind = first_nonzero_index(&self.terms);
166        if ind != 0 {
167            self.terms.drain(0..ind);
168        }
169    }
170
171    pub fn ordered_term_iter(&self) -> impl Iterator<Item = (N, usize)> + '_ {
172        let start = first_nonzero_index(&self.terms);
173        let terms = &self.terms[start..];
174        let deg = terms.len().saturating_sub(1);
175        terms.iter().enumerate().filter_map(move |(index, &coeff)| {
176            if coeff.is_zero() {
177                None
178            } else {
179                Some((coeff, deg - index))
180            }
181        })
182    }
183}
184
185impl<N: Copy + Zero> SizedPolynomial<N> for Polynomial<N> {
186    fn term_with_degree(&self, degree: usize) -> Term<N> {
187        term_with_deg(&self.terms, degree)
188    }
189
190    fn terms_as_vec(&self) -> Vec<(N, usize)> {
191        self.ordered_term_iter().collect()
192    }
193
194    /// Returns the degree of the `Polynomial` it is called on, corresponding to the
195    /// largest non-zero term.
196    ///
197    /// # Example
198    ///
199    /// ```
200    /// use rustnomial::{SizedPolynomial, Polynomial, Degree};
201    /// let polynomial = Polynomial::new(vec![1.0, 4.0, 4.0]);
202    /// assert_eq!(Degree::Num(2), polynomial.degree());
203    /// ```
204    fn degree(&self) -> Degree {
205        degree(&self.terms)
206    }
207
208    /// Returns a `Polynomial` with no terms.
209    ///
210    /// # Example
211    ///
212    /// ```
213    /// use rustnomial::{SizedPolynomial, Polynomial};
214    /// let zero = Polynomial::<i32>::zero();
215    /// assert!(zero.is_zero());
216    /// assert!(zero.ordered_term_iter().next().is_none());
217    /// assert!(zero.terms.is_empty());
218    /// ```
219    fn zero() -> Polynomial<N> {
220        Polynomial { terms: vec![] }
221    }
222
223    /// Sets self to zero.
224    ///
225    /// # Example
226    ///
227    /// ```
228    /// use rustnomial::{Polynomial, SizedPolynomial};
229    /// let mut non_zero = Polynomial::from(vec![0, 1]);
230    /// assert!(!non_zero.is_zero());
231    /// non_zero.set_to_zero();
232    /// assert!(non_zero.is_zero());
233    /// ```
234    fn set_to_zero(&mut self) {
235        self.terms.clear()
236    }
237}
238
239impl<N> MutablePolynomial<N> for Polynomial<N>
240where
241    N: Zero + Copy + AddAssign + SubAssign + CanNegate,
242{
243    fn try_add_term(&mut self, coeff: N, degree: usize) -> Result<(), TryAddError> {
244        Ok(self.add_term(coeff, degree))
245    }
246
247    fn try_sub_term(&mut self, coeff: N, degree: usize) -> Result<(), TryAddError> {
248        if self.terms.len() < degree + 1 {
249            if !N::can_negate() {
250                return Err(TryAddError::CanNotNegate);
251            }
252            let added_zeros = degree + 1 - self.terms.len();
253            self.terms
254                .splice(0..0, core::iter::repeat(N::zero()).take(added_zeros));
255        }
256        let index = self.terms.len() - degree - 1;
257        self.terms[index] -= coeff;
258
259        Ok(())
260    }
261}
262
263impl Polynomial<f64> {
264    /// Return the roots of the `Polynomial`.
265    ///
266    /// # Example
267    ///
268    /// ```
269    /// use rustnomial::{Polynomial, Roots, SizedPolynomial};
270    /// let zero = Polynomial::<f64>::zero();
271    /// assert_eq!(Roots::InfiniteRoots, zero.roots());
272    /// let constant = Polynomial::new(vec![1.]);
273    /// assert_eq!(Roots::NoRoots, constant.roots());
274    /// let monomial = Polynomial::new(vec![1.0, 0.,]);
275    /// assert_eq!(Roots::ManyRealRoots(vec![0.]), monomial.roots());
276    /// let binomial = Polynomial::new(vec![1.0, 2.0]);
277    /// assert_eq!(Roots::ManyRealRoots(vec![-2.0]), binomial.roots());
278    /// let trinomial = Polynomial::new(vec![1.0, 4.0, 4.0]);
279    /// assert_eq!(Roots::ManyRealRoots(vec![-2.0, -2.0]), trinomial.roots());
280    /// let quadnomial = Polynomial::new(vec![1.0, 6.0, 12.0, 8.0]);
281    /// assert_eq!(Roots::ManyRealRoots(vec![-2.0, -2.0, -2.0]), quadnomial.roots());
282    /// ```
283    pub fn roots(self) -> Roots<f64> {
284        find_roots(&self)
285    }
286}
287
288impl<N> FreeSizePolynomial<N> for Polynomial<N>
289where
290    N: Zero + Copy + AddAssign,
291{
292    /// Returns a `Polynomial` with the corresponding terms,
293    /// in order of ax^n + bx^(n-1) + ... + cx + d
294    ///
295    /// # Arguments
296    ///
297    /// * ` terms ` - A slice of (coefficient, degree) pairs.
298    ///
299    /// # Example
300    ///
301    /// ```
302    /// use rustnomial::{FreeSizePolynomial, Polynomial};
303    /// // Corresponds to 1.0x^2 + 4.0x + 4.0
304    /// let polynomial = Polynomial::from_terms(&[(1.0, 2), (4.0, 1), (4.0, 0)]);
305    /// assert_eq!(Polynomial::new(vec![1., 4., 4.]), polynomial);
306    /// ```
307    fn from_terms(terms: &[(N, usize)]) -> Self {
308        let mut a = Polynomial::zero();
309        for &(term, degree) in terms {
310            a.add_term(term, degree);
311        }
312        a
313    }
314
315    fn add_term(&mut self, coeff: N, degree: usize) {
316        if self.terms.len() < degree + 1 {
317            let added_zeros = degree + 1 - self.terms.len();
318            self.terms
319                .splice(0..0, core::iter::repeat(N::zero()).take(added_zeros));
320        }
321        let index = self.terms.len() - degree - 1;
322        self.terms[index] += coeff;
323    }
324}
325
326impl<N> Evaluable<N> for Polynomial<N>
327where
328    N: Zero + Copy + AddAssign + MulAssign + Mul<Output = N>,
329{
330    /// Returns the value of the `Polynomial` at the given point.
331    ///
332    /// # Example
333    ///
334    /// ```
335    /// use rustnomial::{Polynomial, Evaluable};
336    /// let a = Polynomial::new(vec![1, 2, 3, 4]);
337    /// assert_eq!(10, a.eval(1));
338    /// assert_eq!(1234, a.eval(10));
339    /// ```
340    fn eval(&self, point: N) -> N {
341        if let Some((&last, first)) = self.terms.split_last() {
342            if point.is_zero() {
343                return last;
344            }
345
346            // Equivalent to
347            // sum = (a(x^n) + b(x^(n-1)) + ... + z
348            // but instead of having explicit x^n, we do
349            // sum = (((ax) + b)x + ..)x + z
350            // which allows the same result with
351            // the same number of adds, and
352            // n - 1 multiplies vs 2n muls.
353            let mut sum = N::zero();
354            for &val in first.iter() {
355                sum += val;
356                sum *= point;
357            }
358            sum + last
359        } else {
360            N::zero()
361        }
362    }
363}
364
365impl<N> Derivable<N> for Polynomial<N>
366where
367    N: Zero + One + TryFromUsizeContinuous + Copy + MulAssign + SubAssign,
368{
369    /// Returns the derivative of the `Polynomial`.
370    ///
371    /// # Example
372    ///
373    /// ```
374    /// use rustnomial::{Polynomial, Derivable};
375    /// let polynomial = Polynomial::new(vec![4, 1, 5]);
376    /// assert_eq!(Polynomial::new(vec![8, 1]), polynomial.derivative());
377    /// ```
378    ///
379    /// # Errors
380    /// Will panic if `N` can not losslessly encode the numbers from 0 to the degree of `self`.
381    fn derivative(&self) -> Polynomial<N> {
382        let index = first_nonzero_index(&self.terms);
383        let mut degree = N::try_from_usize_cont(self.terms.len() - index)
384            .expect("Degree has no lossless representation in N.");
385        let mut terms = {
386            if let Some((_, terms)) = self.terms.split_at(index).1.split_last() {
387                terms.to_vec()
388            } else {
389                return Polynomial::zero();
390            }
391        };
392        for term in terms.iter_mut() {
393            degree -= N::one();
394            *term *= degree;
395        }
396        Polynomial { terms }
397    }
398}
399
400impl<N> Integrable<N, Polynomial<N>> for Polynomial<N>
401where
402    N: Zero
403        + One
404        + Copy
405        + DivAssign
406        + Mul<Output = N>
407        + MulAssign
408        + AddAssign
409        + TryFromUsizeContinuous
410        + SubAssign,
411{
412    /// Returns the integral of the `Polynomial`.
413    ///
414    /// # Example
415    ///
416    /// ```
417    /// use rustnomial::{Polynomial, Integrable};
418    /// let polynomial = Polynomial::new(vec![1.0, 2.0, 5.0]);
419    /// let integral = polynomial.integral();
420    /// assert_eq!(&Polynomial::new(vec![1.0/3.0, 1.0, 5.0, 0.0]), integral.inner());
421    /// ```
422    ///
423    /// # Errors
424    /// Will panic if `N` can not losslessly encode the numbers from 0 to the degree of self `self`.
425    fn integral(&self) -> Integral<N, Polynomial<N>> {
426        let index = first_nonzero_index(&self.terms);
427        let mut degree = N::try_from_usize_cont(self.terms.len() - index)
428            .expect("Degree can not be losslessly represented.");
429        let mut terms = self.terms[index..].to_vec();
430        for term in terms.iter_mut() {
431            *term /= degree;
432            degree -= N::one();
433        }
434        terms.push(N::zero());
435        Integral::new(Polynomial { terms })
436    }
437}
438
439impl<N> Polynomial<N>
440where
441    N: Mul<Output = N> + AddAssign + Copy + Zero + One,
442{
443    /// Raises the `Polynomial` to the power of exp, using exponentiation by squaring.
444    ///
445    /// # Example
446    ///
447    /// ```
448    /// use rustnomial::Polynomial;
449    /// let polynomial = Polynomial::new(vec![1.0, 2.0]);
450    /// let polynomial_sqr = polynomial.pow(2);
451    /// let polynomial_cub = polynomial.pow(3);
452    /// assert_eq!(polynomial.clone() * polynomial.clone(), polynomial_sqr);
453    /// assert_eq!(polynomial_sqr.clone() * polynomial.clone(), polynomial_cub);
454    /// ```
455    pub fn pow(&self, exp: usize) -> Polynomial<N> {
456        if exp == 0 {
457            Polynomial {
458                terms: vec![N::one(); 1],
459            }
460        } else if exp == 1 {
461            Polynomial::new(self.terms.clone())
462        } else if exp == 2 {
463            self * self
464        } else if exp % 2 == 0 {
465            self.pow(exp / 2).pow(2)
466        } else {
467            self * &self.pow(exp - 1)
468        }
469    }
470}
471
472impl<N> Polynomial<N>
473where
474    N: Copy + Zero + SubAssign + Mul<Output = N> + Div<Output = N>,
475{
476    /// Divides self by the given `Polynomial`, and returns the quotient and remainder.
477    pub fn div_mod(&self, rhs: &Polynomial<N>) -> (Polynomial<N>, Polynomial<N>) {
478        let zero = N::zero();
479
480        let (rhs_first, rhs_deg) = match first_term(&rhs.terms) {
481            Term::ZeroTerm => panic!("Can't divide polynomial by 0."),
482            Term::Term(coeff, deg) => (coeff, deg),
483        };
484
485        let (mut coeff, mut self_degree) = match first_term(&self.terms) {
486            Term::ZeroTerm => {
487                return (Polynomial::zero(), self.clone());
488            }
489            Term::Term(coeff, degree) => {
490                if degree < rhs_deg {
491                    return (Polynomial::zero(), self.clone());
492                }
493                (coeff, degree)
494            }
495        };
496
497        let mut remainder = self.terms.clone();
498        let mut div = vec![zero; self_degree - rhs_deg + 1];
499        let offset = self_degree;
500
501        while self_degree >= rhs_deg {
502            let scale = coeff / rhs_first;
503            vec_sub_w_scale(&mut remainder, self_degree, &rhs.terms, rhs_deg, scale);
504            div[offset - self_degree] = scale;
505            match first_term(&remainder) {
506                Term::ZeroTerm => break,
507                Term::Term(coeffx, degree) => {
508                    coeff = coeffx;
509                    self_degree = degree;
510                }
511            }
512        }
513
514        (Polynomial::new(div), Polynomial::new(remainder))
515    }
516}
517
518impl<N> Polynomial<N>
519where
520    N: Copy + Zero + SubAssign + Mul<Output = N> + Div<Output = N>,
521{
522    /// Divides self by the given `Polynomial`, and returns the quotient.
523    pub fn floor_div(&self, rhs: &Polynomial<N>) -> Polynomial<N> {
524        self.div_mod(rhs).0
525    }
526}
527
528impl<N> Rem<Polynomial<N>> for Polynomial<N>
529where
530    N: Copy + Zero + SubAssign + Mul<Output = N> + Div<Output = N>,
531{
532    type Output = Polynomial<N>;
533
534    /// Returns the remainder of dividing `self` by `rhs`.
535    fn rem(self, rhs: Polynomial<N>) -> Polynomial<N> {
536        let (rhs_first, rhs_deg) = match first_term(&rhs.terms) {
537            Term::ZeroTerm => panic!("Can't divide polynomial by 0."),
538            Term::Term(coeff, deg) => (coeff, deg),
539        };
540
541        let (mut scale, mut self_degree) = match first_term(&self.terms) {
542            Term::ZeroTerm => return self.clone(),
543            Term::Term(coeff, degree) => {
544                if degree < rhs_deg {
545                    return self.clone();
546                }
547                (coeff / rhs_first, degree)
548            }
549        };
550
551        let mut remainder = self.terms.clone();
552
553        while self_degree >= rhs_deg {
554            vec_sub_w_scale(&mut remainder, self_degree, &rhs.terms, rhs_deg, scale);
555            match first_term(&self.terms) {
556                Term::ZeroTerm => break,
557                Term::Term(coeff, degree) => {
558                    scale = coeff / rhs_first;
559                    self_degree = degree;
560                }
561            }
562        }
563
564        Polynomial::new(remainder)
565    }
566}
567
568impl<N> RemAssign<Polynomial<N>> for Polynomial<N>
569where
570    N: Copy + Zero + SubAssign + Mul<Output = N> + Div<Output = N>,
571{
572    /// Assign the remainder of dividing `self` by `rhs` to `self`.
573    fn rem_assign(&mut self, rhs: Polynomial<N>) {
574        let (rhs_first, rhs_deg) = match first_term(&rhs.terms) {
575            Term::ZeroTerm => panic!("Can't divide polynomial by 0."),
576            Term::Term(coeff, deg) => (coeff, deg),
577        };
578
579        let (mut scale, mut self_degree) = match first_term(&self.terms) {
580            Term::ZeroTerm => return,
581            Term::Term(coeff, degree) => {
582                if degree < rhs_deg {
583                    return;
584                }
585                (coeff / rhs_first, degree)
586            }
587        };
588
589        while self_degree >= rhs_deg {
590            vec_sub_w_scale(&mut self.terms, self_degree, &rhs.terms, rhs_deg, scale);
591            match first_term(&self.terms) {
592                Term::ZeroTerm => break,
593                Term::Term(coeff, degree) => {
594                    scale = coeff / rhs_first;
595                    self_degree = degree;
596                }
597            }
598        }
599    }
600}
601
602impl<N> PartialEq for Polynomial<N>
603where
604    N: PartialEq + Zero + Copy,
605{
606    /// Returns true if self and other have the same terms.
607    ///
608    /// # Example
609    ///
610    /// ```
611    /// use rustnomial::Polynomial;
612    /// let a = Polynomial::new(vec![1.0, 2.0]);
613    /// let b = Polynomial::new(vec![2.0, 2.0]);
614    /// let c = Polynomial::new(vec![1.0, 0.0]);
615    /// assert_ne!(a, b);
616    /// assert_ne!(a, c);
617    /// assert_eq!(a, b - c);
618    /// ```
619    fn eq(&self, other: &Self) -> bool {
620        self.ordered_term_iter().eq(other.ordered_term_iter())
621    }
622}
623
624impl<N> From<Vec<N>> for Polynomial<N>
625where
626    N: Copy + Zero,
627{
628    /// Returns a `SparsePolynomial` with the corresponding terms,
629    /// in order of ax^n + bx^(n-1) + ... + cx + d
630    ///
631    /// # Arguments
632    ///
633    /// * ` term_vec ` - A vector of constants, in decreasing order of degree.
634    ///
635    /// # Example
636    ///
637    /// ```
638    /// use rustnomial::{Polynomial};
639    /// // Corresponds to 1.0x^2 + 4.0x + 4.0
640    /// let polynomial = Polynomial::from(vec![1.0, 4.0, 4.0]);
641    /// let polynomial: Polynomial<f64> = vec![1.0, 4.0, 4.0].into();
642    /// ```
643    fn from(term_vec: Vec<N>) -> Self {
644        Polynomial::new(term_vec)
645    }
646}
647
648macro_rules! from_poly_a_to_b {
649    ($A:ty, $B:ty) => {
650        impl From<Polynomial<$A>> for Polynomial<$B> {
651            fn from(item: Polynomial<$A>) -> Self {
652                Polynomial::new(item.terms.into_iter().map(|x| x as $B).collect())
653            }
654        }
655    };
656}
657
658upcast!(from_poly_a_to_b);
659poly_from_str!(Polynomial);
660fmt_poly!(Polynomial);
661
662impl<N> Neg for Polynomial<N>
663where
664    N: Zero + Copy + Neg<Output = N>,
665{
666    type Output = Polynomial<N>;
667
668    fn neg(mut self) -> Polynomial<N> {
669        self.terms.iter_mut().for_each(|x| *x = -*x);
670        self
671    }
672}
673
674impl<N> Sub<Polynomial<N>> for Polynomial<N>
675where
676    N: Zero + Copy + Sub<Output = N> + SubAssign + Neg<Output = N>,
677{
678    type Output = Polynomial<N>;
679
680    fn sub(mut self, rhs: Polynomial<N>) -> Polynomial<N> {
681        self -= rhs;
682        self
683    }
684}
685
686impl<N> SubAssign<Polynomial<N>> for Polynomial<N>
687where
688    N: Neg<Output = N> + Sub<Output = N> + SubAssign + Copy + Zero,
689{
690    fn sub_assign(&mut self, mut rhs: Polynomial<N>) {
691        // This impl eats rhs's terms.
692        if rhs.terms.len() > self.terms.len() {
693            let offset = rhs.terms.len() - self.terms.len();
694            let (right, left) = rhs.terms.split_at_mut(offset);
695
696            right.iter_mut().for_each(|term| *term = -*term);
697            left.iter_mut()
698                .zip(&self.terms)
699                .for_each(|(term, &val)| *term = val - *term);
700            self.terms = rhs.terms;
701        } else {
702            let offset = self.terms.len() - rhs.terms.len();
703            self.terms[offset..]
704                .iter_mut()
705                .zip(rhs.terms)
706                .for_each(|(term, val)| *term -= val);
707        }
708    }
709}
710
711impl<N> Sub<&Polynomial<N>> for Polynomial<N>
712where
713    N: Zero + Copy + Sub<Output = N> + SubAssign + Neg<Output = N>,
714{
715    type Output = Polynomial<N>;
716
717    fn sub(mut self, rhs: &Polynomial<N>) -> Polynomial<N> {
718        self -= rhs;
719        self
720    }
721}
722
723impl<N> SubAssign<&Polynomial<N>> for Polynomial<N>
724where
725    N: Neg<Output = N> + Sub<Output = N> + SubAssign + Copy + Zero,
726{
727    fn sub_assign(&mut self, rhs: &Polynomial<N>) {
728        // This impl eats preprends the relevant terms from rhs into self.
729        if rhs.terms.len() > self.terms.len() {
730            let offset = rhs.terms.len() - self.terms.len();
731            let (right, left) = rhs.terms.split_at(offset);
732            self.terms
733                .iter_mut()
734                .zip(left)
735                .for_each(|(lhs, &rhs)| *lhs = *lhs - rhs);
736            self.terms.splice(0..0, right.iter().map(|coeff| -*coeff));
737        } else {
738            let offset = self.terms.len() - rhs.terms.len();
739            self.terms[offset..]
740                .iter_mut()
741                .zip(&rhs.terms)
742                .for_each(|(term, &val)| *term -= val);
743        }
744    }
745}
746
747impl<N> Add<Polynomial<N>> for Polynomial<N>
748where
749    N: Zero + Copy + AddAssign,
750{
751    type Output = Polynomial<N>;
752
753    fn add(self, rhs: Polynomial<N>) -> Polynomial<N> {
754        let (mut terms, small) = if rhs.terms.len() > self.terms.len() {
755            (rhs.terms, &self.terms)
756        } else {
757            (self.terms, &rhs.terms)
758        };
759
760        let offset = terms.len() - small.len();
761
762        for (index, &val) in terms[offset..].iter_mut().zip(small) {
763            *index += val;
764        }
765
766        Polynomial::new(terms)
767    }
768}
769
770impl<N: Copy + Zero + AddAssign> AddAssign<Polynomial<N>> for Polynomial<N> {
771    fn add_assign(&mut self, rhs: Polynomial<N>) {
772        let lhs = &self.terms;
773        let mut rhs = rhs.terms;
774
775        if rhs.len() > lhs.len() {
776            let offset = rhs.len() - lhs.len();
777            for (index, &val) in rhs[offset..].iter_mut().zip(lhs) {
778                *index += val;
779            }
780            self.terms = rhs;
781        } else {
782            let offset = lhs.len() - rhs.len();
783            for (index, val) in self.terms[offset..].iter_mut().zip(rhs) {
784                *index += val;
785            }
786        }
787    }
788}
789
790impl<N> Mul<Polynomial<N>> for Polynomial<N>
791where
792    N: Mul<Output = N> + AddAssign + Copy + Zero,
793{
794    type Output = Polynomial<N>;
795
796    fn mul(self, rhs: Polynomial<N>) -> Polynomial<N> {
797        Polynomial {
798            terms: slice_mul(&self.terms, &rhs.terms),
799        }
800    }
801}
802
803impl<N> Mul<&Polynomial<N>> for Polynomial<N>
804where
805    N: Mul<Output = N> + AddAssign + Copy + Zero,
806{
807    type Output = Polynomial<N>;
808
809    fn mul(self, rhs: &Polynomial<N>) -> Polynomial<N> {
810        Polynomial::new(slice_mul(&self.terms, &rhs.terms))
811    }
812}
813
814impl<N> Mul<Polynomial<N>> for &Polynomial<N>
815where
816    N: Mul<Output = N> + AddAssign + Copy + Zero,
817{
818    type Output = Polynomial<N>;
819
820    fn mul(self, rhs: Polynomial<N>) -> Polynomial<N> {
821        Polynomial {
822            terms: slice_mul(&self.terms, &rhs.terms),
823        }
824    }
825}
826
827impl<N> Mul<&Polynomial<N>> for &Polynomial<N>
828where
829    N: Mul<Output = N> + AddAssign + Copy + Zero,
830{
831    type Output = Polynomial<N>;
832
833    fn mul(self, rhs: &Polynomial<N>) -> Polynomial<N> {
834        Polynomial::new(slice_mul(&self.terms, &rhs.terms))
835    }
836}
837
838impl<N> MulAssign<Polynomial<N>> for Polynomial<N>
839where
840    N: Mul<Output = N> + AddAssign + Copy + Zero,
841{
842    fn mul_assign(&mut self, rhs: Polynomial<N>) {
843        self.terms = slice_mul(&self.terms, &rhs.terms);
844    }
845}
846
847impl<N> MulAssign<&Polynomial<N>> for Polynomial<N>
848where
849    N: Mul<Output = N> + AddAssign + Copy + Zero,
850{
851    fn mul_assign(&mut self, rhs: &Polynomial<N>) {
852        self.terms = slice_mul(&self.terms, &rhs.terms);
853    }
854}
855
856impl<N: Zero + Copy + Mul<Output = N>> Mul<N> for Polynomial<N> {
857    type Output = Polynomial<N>;
858
859    fn mul(self, rhs: N) -> Polynomial<N> {
860        Polynomial::new(self.terms.iter().map(|&x| x * rhs).collect())
861    }
862}
863
864impl<N: Copy + MulAssign> MulAssign<N> for Polynomial<N> {
865    fn mul_assign(&mut self, rhs: N) {
866        for p in self.terms.iter_mut() {
867            *p *= rhs;
868        }
869    }
870}
871
872impl<N> Div<N> for Polynomial<N>
873where
874    N: Zero + Copy + Div<Output = N>,
875{
876    type Output = Polynomial<N>;
877
878    fn div(self, rhs: N) -> Polynomial<N> {
879        Polynomial::new(self.terms.iter().map(|&x| x / rhs).collect())
880    }
881}
882
883impl<N: Copy + DivAssign> DivAssign<N> for Polynomial<N> {
884    fn div_assign(&mut self, rhs: N) {
885        for p in self.terms.iter_mut() {
886            *p /= rhs;
887        }
888    }
889}
890
891impl<N: Zero + Copy> Shl<i32> for Polynomial<N> {
892    type Output = Polynomial<N>;
893
894    fn shl(self, rhs: i32) -> Polynomial<N> {
895        if rhs < 0 {
896            self >> -rhs
897        } else {
898            let index = first_nonzero_index(&self.terms);
899            let mut terms = self.terms[index..].to_vec();
900            terms.extend(vec![N::zero(); rhs as usize]);
901            Polynomial { terms }
902        }
903    }
904}
905
906impl<N: Zero + Copy> ShlAssign<i32> for Polynomial<N> {
907    fn shl_assign(&mut self, rhs: i32) {
908        if rhs < 0 {
909            *self >>= -rhs;
910        } else {
911            self.terms.extend(vec![N::zero(); rhs as usize]);
912        }
913    }
914}
915
916impl<N: Zero + Copy> Shr<i32> for Polynomial<N> {
917    type Output = Polynomial<N>;
918
919    fn shr(self, rhs: i32) -> Polynomial<N> {
920        if rhs < 0 {
921            self << -rhs
922        } else {
923            let rhs = rhs as usize;
924            let index = first_nonzero_index(&self.terms);
925            Polynomial {
926                terms: if rhs > self.terms.len() {
927                    vec![]
928                } else {
929                    self.terms[index..self.terms.len() - rhs].to_vec()
930                },
931            }
932        }
933    }
934}
935
936impl<N: Zero + Copy> ShrAssign<i32> for Polynomial<N> {
937    fn shr_assign(&mut self, rhs: i32) {
938        if rhs < 0 {
939            *self <<= -rhs;
940        } else {
941            let rhs = rhs as usize;
942            if rhs > self.terms.len() {
943                self.terms = vec![];
944            } else {
945                self.terms = self.terms[..self.terms.len() - rhs].to_vec();
946            }
947        }
948    }
949}
950
951/// TODO:
952/// modulo floordiv
953#[cfg(test)]
954mod test {
955    use crate::{
956        polynomial, Degree, Derivable, Evaluable, Integrable, Polynomial, SizedPolynomial,
957    };
958
959    #[test]
960    fn test_polynomial_macro() {
961        assert_eq!(polynomial!(1, 2, 3, 4), Polynomial::new(vec![1, 2, 3, 4]));
962    }
963
964    #[test]
965    fn test_eval() {
966        let a = Polynomial::new(vec![1, 2, 3]);
967        assert_eq!(25 + 2 * 5 + 3, a.eval(5));
968    }
969
970    #[test]
971    fn test_derivative() {
972        let a = Polynomial::new(vec![1, 2, 3]);
973        let b = Polynomial::new(vec![2, 2]);
974        assert_eq!(b, a.derivative());
975
976        let a = Polynomial::new(vec![0, 1, 2, 3]);
977        assert_eq!(b, a.derivative());
978
979        let a = Polynomial::new(vec![1, 2, 3, 4]);
980        let b = Polynomial::new(vec![3, 4, 3]);
981        assert_eq!(b, a.derivative());
982
983        assert_eq!(Polynomial::<i32>::zero(), Polynomial::zero().derivative());
984    }
985
986    #[test]
987    fn test_integral() {
988        let a = Polynomial::new(vec![3, 2, 1]);
989        let b = Polynomial::new(vec![1, 1, 1, 0]);
990        assert_eq!(&b, a.integral().inner());
991    }
992
993    #[test]
994    fn test_integral_eval() {
995        let a = Polynomial::new(vec![3, 2, 1]);
996        assert_eq!(3, a.integral().eval(0, 1));
997    }
998
999    #[test]
1000    fn test_integral_const_substitute() {
1001        let a = Polynomial::new(vec![3, 2, 1]);
1002        let b = Polynomial::new(vec![1, 1, 1, 5]);
1003        assert_eq!(b, a.integral().replace_c(5));
1004    }
1005
1006    #[test]
1007    fn test_add_lhs_bigger() {
1008        let a = Polynomial::new(vec![1, 2, 3]);
1009        let b = Polynomial::new(vec![1, 2, 3, 4]);
1010        let c = Polynomial::new(vec![1, 3, 5, 7]);
1011        assert_eq!(c, b + a);
1012    }
1013
1014    #[test]
1015    fn test_add_rhs_bigger() {
1016        let a = Polynomial::new(vec![1, 2, 3]);
1017        let b = Polynomial::new(vec![1, 2, 3, 4]);
1018        let c = Polynomial::new(vec![1, 3, 5, 7]);
1019        assert_eq!(c, a + b);
1020    }
1021
1022    #[test]
1023    fn test_add_lhs_bigger_assign() {
1024        let a = Polynomial::new(vec![1, 2, 3]);
1025        let mut b = Polynomial::new(vec![1, 2, 3, 4]);
1026        b += a;
1027        let c = Polynomial::new(vec![1, 3, 5, 7]);
1028        assert_eq!(c, b);
1029    }
1030
1031    #[test]
1032    fn test_add_rhs_bigger_assign() {
1033        let mut a = Polynomial::new(vec![1, 2, 3]);
1034        let b = Polynomial::new(vec![1, 2, 3, 4]);
1035        a += b;
1036        let c = Polynomial::new(vec![1, 3, 5, 7]);
1037        assert_eq!(c, a);
1038    }
1039
1040    #[test]
1041    fn test_sub_lhs_bigger() {
1042        let a = Polynomial::new(vec![2, 3, 4]);
1043        let b = Polynomial::new(vec![1, 2, 3, 4]);
1044        let c = Polynomial::new(vec![1, 0, 0, 0]);
1045        assert_eq!(c, b - a);
1046    }
1047
1048    #[test]
1049    fn test_sub_rhs_bigger() {
1050        let a = Polynomial::new(vec![2, 3, 4]);
1051        let b = Polynomial::new(vec![1, 2, 3, 4]);
1052        let c = Polynomial::new(vec![-1, 0, 0, 0]);
1053        assert_eq!(c, a - b);
1054    }
1055
1056    #[test]
1057    fn test_sub_lhs_bigger_assign() {
1058        let a = Polynomial::new(vec![2, 3, 4]);
1059        let mut b = Polynomial::new(vec![1, 2, 3, 4]);
1060        b -= a;
1061        let c = Polynomial::new(vec![1, 0, 0, 0]);
1062        assert_eq!(c, b);
1063    }
1064
1065    #[test]
1066    fn test_sub_rhs_bigger_assign() {
1067        let mut a = Polynomial::new(vec![2, 3, 4]);
1068        let b = Polynomial::new(vec![1, 2, 3, 4]);
1069        a -= b;
1070        let c = Polynomial::new(vec![-1, 0, 0, 0]);
1071        assert_eq!(c, a);
1072    }
1073
1074    #[test]
1075    fn test_negate() {
1076        let a = Polynomial::new(vec![1, 2, 3, 0, -5]);
1077        let c = Polynomial::new(vec![-1, -2, -3, 0, 5]);
1078        assert_eq!(c, -a);
1079    }
1080
1081    #[test]
1082    fn test_mul_poly() {
1083        let a = Polynomial::new(vec![1, 2]);
1084        let b = a.clone();
1085        let c = Polynomial::new(vec![1, 4, 4]);
1086        assert_eq!(c, a * b);
1087    }
1088
1089    #[test]
1090    fn test_mul_assign_poly() {
1091        let mut a = Polynomial::new(vec![1, 2]);
1092        let b = a.clone();
1093        a *= b;
1094        let c = Polynomial::new(vec![1, 4, 4]);
1095        assert_eq!(c, a);
1096    }
1097
1098    #[test]
1099    fn test_mul_num() {
1100        let a = Polynomial::new(vec![1, 2]);
1101        let c = Polynomial::new(vec![10, 20]);
1102        assert_eq!(c, a * 10);
1103    }
1104
1105    #[test]
1106    fn test_mul_assign_num() {
1107        let mut a = Polynomial::new(vec![1, 2]);
1108        a *= 10;
1109        let c = Polynomial::new(vec![10, 20]);
1110        assert_eq!(c, a);
1111    }
1112
1113    #[test]
1114    fn test_equality() {
1115        let a = Polynomial::new(vec![1, 2]);
1116        let mut c = Polynomial::new(vec![0, 0, 0, 1, 2]);
1117        c.terms = vec![0, 0, 0, 1, 2];
1118
1119        assert_eq!(c, a);
1120
1121        c.terms = vec![1, 2, 0, 0, 0];
1122
1123        assert_ne!(c, a);
1124    }
1125
1126    #[test]
1127    fn test_equality_first_match() {
1128        let a = Polynomial::new(vec![1, 2]);
1129        let b = Polynomial::new(vec![1, 0]);
1130        assert_ne!(a, b);
1131    }
1132
1133    #[test]
1134    fn test_equality_different() {
1135        let a = Polynomial::new(vec![1, 2]);
1136        let b = Polynomial::new(vec![3, 7, 4]);
1137        assert_ne!(a, b);
1138    }
1139
1140    #[test]
1141    fn test_shl_pos() {
1142        let a = Polynomial::new(vec![1, 2]);
1143        let c = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1144        assert_eq!(c, a << 5);
1145    }
1146
1147    #[test]
1148    fn test_shl_assign_pos() {
1149        let mut a = Polynomial::new(vec![1, 2]);
1150        a <<= 5;
1151        let c = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1152        assert_eq!(c, a);
1153    }
1154
1155    #[test]
1156    fn test_shl_neg() {
1157        let a = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1158        let c = Polynomial::new(vec![1, 2]);
1159        assert_eq!(c, a << -5);
1160    }
1161
1162    #[test]
1163    fn test_shl_assign_neg() {
1164        let mut a = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1165        a <<= -5;
1166        let c = Polynomial::new(vec![1, 2]);
1167        assert_eq!(c, a);
1168    }
1169
1170    #[test]
1171    fn test_shr_pos() {
1172        let a = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1173        let c = Polynomial::new(vec![1, 2]);
1174        assert_eq!(c, a >> 5);
1175    }
1176
1177    #[test]
1178    fn test_shr_assign_pos() {
1179        let mut a = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1180        a >>= 5;
1181        let c = Polynomial::new(vec![1, 2]);
1182        assert_eq!(c, a);
1183    }
1184
1185    #[test]
1186    fn test_shr_neg() {
1187        let a = Polynomial::new(vec![1, 2]);
1188        let c = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1189        assert_eq!(c, a >> -5);
1190    }
1191
1192    #[test]
1193    fn test_shr_assign_neg() {
1194        let mut a = Polynomial::new(vec![1, 2]);
1195        a >>= -5;
1196        let c = Polynomial::new(vec![1, 2, 0, 0, 0, 0, 0]);
1197        assert_eq!(c, a);
1198    }
1199
1200    #[test]
1201    fn test_shr_to_zero() {
1202        let a = Polynomial::new(vec![1, 2]);
1203        assert_eq!(Polynomial::zero(), a >> 5);
1204    }
1205
1206    #[test]
1207    fn test_shr_assign_to_zero() {
1208        let mut a = Polynomial::new(vec![1, 2]);
1209        a >>= 5;
1210        assert_eq!(Polynomial::zero(), a);
1211    }
1212
1213    #[test]
1214    fn test_exp() {
1215        let a = &Polynomial::new(vec![1, 2]);
1216        let mut b = a.clone();
1217        assert_eq!(Polynomial::new(vec![1]), a.pow(0));
1218        for i in 1..10 {
1219            assert_eq!(b, a.pow(i));
1220            b *= a;
1221        }
1222    }
1223
1224    #[test]
1225    fn test_trim() {
1226        let input_ouput_vecs = vec![
1227            (vec![0, 1, 2], vec![1, 2]),
1228            (vec![0; 10000], vec![]),
1229            (vec![], vec![]),
1230        ];
1231        for (test, expected) in input_ouput_vecs.into_iter() {
1232            let a = Polynomial::new(test);
1233            assert_eq!(expected, a.terms)
1234        }
1235    }
1236    #[test]
1237    fn test_degree() {
1238        let a = Polynomial::new(vec![0, 0, 0, -1, -2, 3]);
1239        assert_eq!(Degree::Num(2), a.degree());
1240    }
1241}
rustnomial/polynomial/polynomial.rs

rustnomial/polynomial/
polynomial.rs
