polynomen 1.0.0

//! # Polynomials Library
//!
//! Polynomial implementation
//! * builder from coefficients or roots
//! * degree
//! * extend by adding 0 coefficients to higher order terms
//! * transformation to monic form
//! * rounding to zero for small coefficients
//! * coefficient indexing
//! * zero and unit polynomials
//! * arithmetic operations between polynomials (addition, subtraction,
//!   multiplication, division, reminder, negation)
//! * arithmetic operations with floats (addition, subtraction,
//!   multiplication, division)
//! * multiplication using fast fourier transform
//! * polynomial exponentiation
//! * differentiation and integration
//! * polynomial evaluation
//! * evaluation of polynomial ratios that reduces overflows
//! * greatest common divisor between two polynomials
//! * polynomial norms (l1, l2, l∞)
//! * roots finding (real and complex) using eigenvalues of the companion matrix or iterative method

#![warn(
    rustdoc::missing_crate_level_docs,
    missing_debug_implementations,
    missing_docs,
    unreachable_pub
)]

#[cfg(test)]
#[macro_use]
extern crate approx;

use std::{
    fmt::{Debug, Formatter},
    ops::{Add, Div, Index, IndexMut, Mul, Neg, Sub},
};

pub mod arithmetic;
mod calculus;
mod complex;
mod convex_hull;
mod eval;
mod fft;
mod gcd;
mod norms;
mod roots;
mod traits;

pub use eval::eval_poly_ratio;
pub use traits::*;

/// Polynomial object
///
/// Contains the vector of coefficients form the lowest to the highest degree
///
/// `p(x) = c0 + c1*x + c2*x^2 + ...`
#[derive(Debug, PartialEq, Clone)]
pub struct Poly<T> {
    coeffs: Vec<T>,
}

/// Macro shortcut to crate a polynomial from its coefficients.
///
/// # Example
/// ```
/// #[macro_use] extern crate polynomen;
/// use polynomen::Poly;
/// let p1 = poly!(1, 2, 3);
/// let p2 = Poly::new_from_coeffs(&[1, 2, 3]);
/// assert_eq!(p1, p2);
/// ```
#[macro_export]
macro_rules! poly {
    ($($c:expr),+ $(,)*) => {
        $crate::Poly::new_from_coeffs(&[$($c,)*])
    };
}

impl<T> Poly<T> {
    /// Length of the polynomial coefficients
    fn len(&self) -> usize {
        self.coeffs.len()
    }

    /// Return the coefficients of the polynomial as a slice
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let c = &[1., 2., 3.];
    /// let p = Poly::new_from_coeffs(c);
    /// assert_eq!(c, p.as_slice());
    /// ```
    #[must_use]
    pub fn as_slice(&self) -> &[T] {
        self.as_ref()
    }
}

impl<T: Clone> Poly<T> {
    /// Vector copy of the polynomial's coefficients
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs(&[1., 2., 3.]);
    /// assert_eq!(vec![1., 2., 3.], p.coeffs());
    /// ```
    #[must_use]
    pub fn coeffs(&self) -> Vec<T> {
        self.coeffs.clone()
    }
}

impl<T: Clone + PartialEq + Zero> Poly<T> {
    /// Create a new polynomial given a slice of real coefficients.
    /// It trims any leading zeros in the high order coefficients.
    ///
    /// # Arguments
    ///
    /// * `coeffs` - slice of coefficients
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs(&[1., 2., 3.]);
    /// ```
    pub fn new_from_coeffs(coeffs: &[T]) -> Self {
        let p = Self {
            coeffs: coeffs.into(),
        };
        p.trim()
    }

    /// Create a new polynomial given a iterator of real coefficients.
    /// It trims any leading zeros in the high order coefficients.
    ///
    /// # Arguments
    ///
    /// * `coeffs` - iterator of coefficients
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs_iter(1..4);
    /// ```
    pub fn new_from_coeffs_iter<II>(coeffs: II) -> Self
    where
        II: IntoIterator<Item = T>,
    {
        let p = Self {
            coeffs: coeffs.into_iter().collect(),
        };
        p.trim()
    }

    /// Trim the zeros coefficients of high degree terms.
    /// It will not leave an empty `coeffs` vector: zero poly is returned.
    #[must_use]
    fn trim(mut self) -> Self {
        self.trim_mut();
        self
    }

    /// Trim the zeros coefficients of high degree terms.
    /// It will not leave an empty `coeffs` vector: zero poly is returned.
    /// Modification in place
    fn trim_mut(&mut self) {
        if let Some(p) = self.coeffs.iter().rposition(|c| c != &T::zero()) {
            let new_length = p + 1;
            debug_assert!(new_length > 0);
            self.coeffs.truncate(new_length);
        } else {
            self.coeffs.resize(1, T::zero());
        }
        debug_assert!(!self.coeffs.is_empty());
    }

    /// Degree of the polynomial
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs(&[1., 2., 3.]);
    /// assert_eq!(Some(2), p.degree());
    /// ```
    #[must_use]
    pub fn degree(&self) -> Option<usize> {
        debug_assert!(
            !self.coeffs.is_empty(),
            "Degree is not defined on empty polynomial"
        );
        if self.is_zero() {
            None
        } else {
            Some(self.coeffs.len() - 1)
        }
    }

    /// Extend the polynomial coefficients with 0 to the given degree in place.
    /// It does not truncate the polynomial.
    ///
    /// # Arguments
    ///
    /// * `degree` - Degree of the new highest coefficient.
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let mut p = Poly::new_from_coeffs(&[1, 2, 3]);
    /// p.extend(5);
    /// assert_eq!(vec![1, 2, 3, 0, 0, 0], p.coeffs());
    /// ```
    pub fn extend(&mut self, degree: usize) {
        match self.degree() {
            None => self.coeffs.resize(degree + 1, T::zero()),
            Some(d) if degree > d => self.coeffs.resize(degree + 1, T::zero()),
            Some(_) => (),
        };
        debug_assert!(!self.coeffs.is_empty());
    }
}

impl<T: Clone + PartialEq + Zero> Poly<T> {
    /// Truncate the polynomial to the given degree.
    /// No modifications if `degree` is higher than or equal to the degree
    /// of the polynomial
    pub(crate) fn truncate(&mut self, degree: usize) {
        if Some(degree) < self.degree() {
            self.coeffs.truncate(degree + 1);
        }
    }
}

impl<T> Poly<T> {
    /// Apply a method to every coefficient of the polynomial and reduce to a single value.
    ///
    /// # Arguments
    ///
    /// * `map` - Method applied to every coefficient.
    /// * `reduce` - Method that reduce the coefficients to a single value.
    ///
    /// # Panics
    ///
    /// The method panics if the polynomial is empty, it should be guaranteed that
    /// the polynomial is never empty.
    fn map_reduce<M, R, U>(&self, map: M, reduce: R) -> U
    where
        M: FnMut(&T) -> U,
        R: FnMut(U, U) -> U,
    {
        self.coeffs.iter().map(map).reduce(reduce).unwrap()
    }
}

impl<T: Clone + Div<Output = T> + One + PartialEq + Zero> Poly<T> {
    /// Return the monic polynomial and the leading coefficient.
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs(&[1., 2., 10.]);
    /// let (p2, c) = p.monic();
    /// assert_eq!(Poly::new_from_coeffs(&[0.1, 0.2, 1.]), p2);
    /// assert_eq!(10., c);
    /// ```
    #[must_use]
    pub fn monic(&self) -> (Self, T) {
        let lc = self.leading_coeff();
        let monic_poly = self / lc.clone();
        (monic_poly, lc)
    }
}

impl<T: Clone + Div<Output = T> + One + PartialEq + Zero> Poly<T> {
    /// Return the monic polynomial and the leading coefficient,
    /// it mutates the polynomial in place.
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let mut p = Poly::new_from_coeffs(&[1., 2., 10.]);
    /// let c = p.monic_mut();
    /// assert_eq!(Poly::new_from_coeffs(&[0.1, 0.2, 1.]), p);
    /// assert_eq!(10., c);
    /// ```
    pub fn monic_mut(&mut self) -> T {
        let lc = self.leading_coeff();
        self.div_mut(&lc);
        lc
    }
}

impl<T: Clone + One> Poly<T> {
    /// Return the leading coefficient of the polynomial.
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs(&[1., 2., 10.]);
    /// let c = p.leading_coeff();
    /// assert_eq!(10., c);
    /// ```
    #[must_use]
    pub fn leading_coeff(&self) -> T {
        self.coeffs.last().unwrap_or(&T::one()).clone()
    }
}

impl<T: Add<Output = T> + Clone + Mul<Output = T> + Neg<Output = T> + One + PartialEq + Zero>
    Poly<T>
{
    /// Create a new polynomial given a slice of real roots
    /// It trims any leading zeros in the high order coefficients.
    ///
    /// # Arguments
    ///
    /// * `roots` - slice of roots
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_roots(&[1., 2., 3.]);
    /// ```
    pub fn new_from_roots(roots: &[T]) -> Self {
        let p = roots.iter().fold(Self::one(), |acc, r| {
            acc * Self {
                coeffs: vec![-r.clone(), T::one()],
            }
        });
        p.trim()
    }

    /// Create a new polynomial given an iterator of real roots
    /// It trims any leading zeros in the high order coefficients.
    ///
    /// # Arguments
    ///
    /// * `roots` - iterator of roots
    ///
    /// # Example
    /// ```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_roots_iter((1..4));
    /// ```
    pub fn new_from_roots_iter<II>(roots: II) -> Self
    where
        II: IntoIterator<Item = T>,
    {
        let p = roots.into_iter().fold(Self::one(), |acc, r| {
            acc * Self {
                coeffs: vec![-r, T::one()],
            }
        });
        p.trim()
    }
}

impl<T: Abs + Clone + PartialEq + PartialOrd + Zero> Poly<T> {
    /// Round off to zero coefficients smaller than `atol`.
    ///
    /// # Arguments
    ///
    /// * `atol` - Absolute tolerance (should be positive)
    ///
    /// # Example
    ///```
    /// use polynomen::Poly;
    /// let p = Poly::new_from_coeffs(&[1., 0.002, 1., -0.0001]);
    /// let actual = p.roundoff(&0.01);
    /// let expected = Poly::new_from_coeffs(&[1., 0., 1.]);
    /// assert_eq!(expected, actual);
    ///```
    pub fn roundoff(&self, atol: &T) -> Self {
        let atol = atol.abs();
        let new_coeff = self
            .coeffs
            .iter()
            .map(|c| if c.abs() < atol { T::zero() } else { c.clone() });
        Poly::new_from_coeffs_iter(new_coeff)
    }

    /// Round off to zero coefficients smaller than `atol` in place.
    ///
    /// # Arguments
    ///
    /// * `atol` - Absolute tolerance (should be positive)
    ///
    /// # Example
    ///```
    /// use polynomen::Poly;
    /// let mut p = Poly::new_from_coeffs(&[1., 0.002, 1., -0.0001]);
    /// p.roundoff_mut(&0.01);
    /// let expected = Poly::new_from_coeffs(&[1., 0., 1.]);
    /// assert_eq!(expected, p);
    ///```
    pub fn roundoff_mut(&mut self, atol: &T) {
        let atol = atol.abs();
        for c in &mut self.coeffs {
            if c.abs() < atol {
                *c = T::zero();
            }
        }
        self.trim_mut();
    }
}

/// Implement read only indexing of polynomial returning its coefficients.
///
/// # Panics
///
/// Panics for out of bounds access.
///
/// # Example
/// ```
/// use polynomen::Poly;
/// let p = Poly::new_from_coeffs(&[0, 1, 2, 3]);
/// assert_eq!(2, p[2]);
/// ```
impl<T> Index<usize> for Poly<T> {
    type Output = T;

    fn index(&self, i: usize) -> &Self::Output {
        &self.coeffs[i]
    }
}

/// Implement mutable indexing of polynomial returning its coefficients.
/// Misuse of this method may invalidate some polynomial invariant.
/// If some coefficient is zeroed a call to `roundoff` of `roundoff_mut` may be needed.
///
/// # Panics
///
/// Panics for out of bounds access.
///
/// # Example
/// ```
/// use polynomen::Poly;
/// let mut p = Poly::new_from_coeffs(&[0, 1, 2, 3]);
/// p[2] = 4;
/// assert_eq!(4, p[2]);
/// ```
impl<T> IndexMut<usize> for Poly<T> {
    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
        &mut self.coeffs[i]
    }
}

/// Implementation of the additive identity for polynomials
///
/// # Example
/// ```
/// use polynomen::{Poly, Zero};
/// let zero = Poly::<u8>::zero();
/// assert!(zero.is_zero());
/// ```
impl<T: PartialEq + Zero> Zero for Poly<T> {
    fn zero() -> Self {
        // The polynomial is never empty.
        Self {
            coeffs: vec![T::zero()],
        }
    }

    fn is_zero(&self) -> bool {
        self.coeffs.len() == 1 && self.coeffs[0] == T::zero()
    }
}

/// Implementation of the multiplicative identity for polynomials
///
/// # Example
/// ```
/// use polynomen::{One, Poly};
/// let one = Poly::<u8>::one();
/// assert!(one.is_one());
/// ```
impl<T: One + PartialEq> One for Poly<T> {
    fn one() -> Self {
        // The polynomial is never empty.
        Self {
            coeffs: vec![T::one()],
        }
    }

    fn is_one(&self) -> bool {
        self.coeffs.len() == 1 && self.coeffs[0] == T::one()
    }
}

/// Implement printing of polynomial
///
/// # Example
/// ```
/// use polynomen::Poly;
/// let p = Poly::new_from_coeffs(&[0, 1, 2, 0, 3]);
/// assert_eq!("1s +2s^2 +3s^4", format!("{}", p));
/// ```
macro_rules! display {
    ($trait:path) => {
        impl<T: $trait + PartialOrd + Zero> $trait for Poly<T> {
            fn fmt(&self, f: &mut Formatter) -> std::fmt::Result {
                debug_assert!(!self.coeffs.is_empty());
                if self.len() == 1 {
                    return self[0].fmt(f);
                }

                let iter = self
                    .coeffs
                    .iter()
                    .enumerate()
                    .filter(|(_, x)| !x.is_zero())
                    .enumerate();
                for (i, (n, c)) in iter {
                    match (i, f.sign_plus(), c < &T::zero()) {
                        (0, _, _) => (),
                        (_, false, false) => write!(f, " +")?,
                        (_, _, _) => write!(f, " ")?,
                    }
                    if n == 0 {
                        c.fmt(f)?;
                    } else if n == 1 {
                        c.fmt(f)?;
                        write!(f, "s")?;
                    } else {
                        c.fmt(f)?;
                        write!(f, "s^")?;
                        write!(f, "{}", n)?;
                    }
                }
                write!(f, "")
            }
        }
    };
}

display!(std::fmt::Binary);
display!(std::fmt::Display);
display!(std::fmt::LowerExp);
display!(std::fmt::LowerHex);
display!(std::fmt::Octal);
display!(std::fmt::UpperExp);
display!(std::fmt::UpperHex);

// TODO: this trait implementation works from Rust 1.41.
// It is similar to the method .coeffs().
// I keep it commented if the will be more features that require newer
// compiler version I will decomment it.
// /// Conversion from `Poly` to a `Vec` containing its coefficients.
// impl<T> From<Poly<T>> for Vec<T> {
//     fn from(poly: Poly<T>) -> Self {
//         poly.coeffs
//     }
// }

/// View the `Poly` coefficients as slice.
impl<T> AsRef<[T]> for Poly<T> {
    fn as_ref(&self) -> &[T] {
        self.coeffs.as_ref()
    }
}

/// Calculate the complex roots of the quadratic equation x^2 + b*x + c = 0.
///
/// # Arguments
///
/// * `b` - first degree coefficient
/// * `c` - zero degree coefficient
///
/// # Example
///```
/// let actual = polynomen::complex_quadratic_roots(0., 1.);
/// assert_eq!(((0., -1.), (0., 1.)), actual);
///```
pub fn complex_quadratic_roots<T>(b: T, c: T) -> ((T, T), (T, T))
where
    T: Abs
        + Add<Output = T>
        + Clone
        + Div<Output = T>
        + Inv
        + Mul<Output = T>
        + Neg<Output = T>
        + One
        + PartialOrd
        + Pow
        + Sign
        + Sqrt
        + Sub<Output = T>
        + Zero,
{
    let (r1, r2) = roots::complex_quadratic_roots_impl(b, c);
    ((r1.re.x, r1.im.x), (r2.re.x, r2.im.x))
}

/// Calculate the real roots of the quadratic equation x^2 + b*x + c = 0.
///
/// # Arguments
///
/// * `b` - first degree coefficient
/// * `c` - zero degree coefficient
///
/// # Example
///```
/// use polynomen;
/// let actual = polynomen::real_quadratic_roots(-2., 1.);
/// assert_eq!(Some((1., 1.)), actual);
///```
pub fn real_quadratic_roots<T>(b: T, c: T) -> Option<(T, T)>
where
    T: Add<Output = T>
        + Clone
        + Div<Output = T>
        + Mul<Output = T>
        + Neg<Output = T>
        + One
        + Pow
        + Sign
        + Sqrt
        + Sub<Output = T>
        + Zero,
{
    roots::real_quadratic_roots_impl(b, c)
}

#[cfg(test)]
mod tests {
    use proptest::prelude::*;

    use super::*;

    #[test]
    fn poly_formatting() {
        assert_eq!("0", format!("{}", Poly::<i16>::zero()));
        assert_eq!("0", format!("{}", Poly::<u16>::new_from_coeffs(&[])));
        assert_eq!("1 +2s^3 -4s^4", format!("{}", poly!(1, 0, 0, 2, -4)));
        assert_eq!("1.235", format!("{:.3}", Poly::new_from_coeffs(&[1.23456])));
        let p = poly!(1.2345, -5.4321, 13.1234);
        assert_eq!("+1.23 -5.43s +13.12s^2", format!("{:+.2}", &p));
        assert_eq!("1.23 -5.43s +13.12s^2", format!("{:.2}", &p));
        assert_eq!("1.2345e0 -5.4321e0s +1.31234e1s^2", format!("{:e}", &p));
    }

    #[test]
    fn poly_creation_coeffs() {
        let c = [4.3, 5.32];
        for (c1, c2) in c.iter().zip(Poly::new_from_coeffs(&c).coeffs) {
            assert_relative_eq!(*c1, c2);
        }

        let c2 = [0., 1., 1., 0., 0., 0.];
        for (i, c) in c2[..3].iter().zip(Poly::new_from_coeffs(&c2).coeffs) {
            assert_relative_eq!(*i, c);
        }

        let zero: [f64; 1] = [0.];
        for (z, c) in zero.iter().zip(poly!(0., 0.).coeffs) {
            assert_relative_eq!(*z, c);
        }

        let int = [1, 2, 3, 4, 5];
        assert_eq!(int, Poly::new_from_coeffs(&int).coeffs.as_slice());

        let float = [0.1_f32, 0.34, 3.43];
        for (f, c) in float.iter().zip(Poly::new_from_coeffs(&float).coeffs) {
            assert_relative_eq!(*f, c);
        }

        assert_eq!(
            Poly::new_from_coeffs(&[1, 2, 3]),
            Poly::new_from_coeffs_iter(1..=3)
        );
    }

    #[test]
    fn coeffs() {
        let int = [1, 2, 3, 4, 5];
        let p = Poly::new_from_coeffs(&int);
        assert_eq!(int, p.coeffs().as_slice());
    }

    #[test]
    fn as_slice() {
        let int = [1, 2, 3, 4, 5];
        let p = Poly::new_from_coeffs(&int);
        assert_eq!(int, p.as_slice());
    }

    #[test]
    fn poly_creation_roots() {
        assert_eq!(poly!(4.0_f64, 4., 1.), Poly::new_from_roots(&[-2., -2.]));

        assert_eq!(poly!(4, 4, 1), Poly::new_from_roots(&[-2, -2]));

        assert!(vec![-2., -2.]
            .iter()
            .zip(
                Poly::new_from_roots(&[-2.0_f64, -2.])
                    .real_roots()
                    .unwrap()
                    .iter()
            )
            .map(|(x, y): (&f64, &f64)| (x - y).abs())
            .all(|x| x < 0.000_001));

        assert!(vec![1.0_f32, 2., 3.]
            .iter()
            .zip(
                Poly::new_from_roots(&[1.0_f32, 2., 3.])
                    .real_roots()
                    .unwrap()
                    .iter()
            )
            .map(|(x, y): (&f32, &f32)| (x - y).abs())
            .all(|x| x < 0.000_01));

        assert_eq!(
            poly!(0., -2., 1., 1.),
            Poly::new_from_roots(&[-0., -2., 1.])
        );

        let v = vec![1, 2, -3];
        assert_eq!(Poly::new_from_roots(&v), Poly::new_from_roots_iter(v));
    }

    #[test]
    fn len() {
        let p = Poly::new_from_coeffs(&[1., 2., 3.]);
        assert_eq!(3, p.len());
    }

    #[test]
    fn degree() {
        let p = Poly::new_from_coeffs(&[1., 2., 3.]);
        assert_eq!(Some(2), p.degree());

        let p2 = Poly::new_from_coeffs(&[0.]);
        assert_eq!(None, p2.degree());
    }

    #[test]
    fn extend_less() {
        let mut p1 = poly!(3, 4, 2);
        let p2 = p1.clone();
        p1.extend(1);
        assert_eq!(p1, p2);
    }

    #[test]
    fn extend_more() {
        let mut p1 = poly!(3, 4, 2);
        let p2 = Poly {
            coeffs: vec![3, 4, 2, 0, 0, 0, 0],
        };
        p1.extend(6);
        assert_eq!(p1, p2);
    }

    #[test]
    fn extend_zero() {
        let mut p1 = Poly::<u32>::zero();
        let p2 = Poly {
            coeffs: vec![0, 0, 0, 0],
        };
        p1.extend(3);
        assert_eq!(p1, p2);
    }

    #[test]
    fn truncate() {
        let mut p1 = poly!(3, 4, 2);
        let p2 = poly!(3, 4);
        p1.truncate(1);
        assert_eq!(p1, p2);
    }

    #[test]
    fn indexing() {
        assert_abs_diff_eq!(3., poly!(1., 3.)[1], epsilon = 0.);

        let mut p = Poly::new_from_roots(&[1., 4., 5.]);
        p[2] = 3.;
        assert_eq!(poly!(-20., 29., 3., 1.), p);
    }

    #[test]
    fn float_coeffs_identities() {
        assert!(Poly::<f64>::zero().is_zero());
        assert!(Poly::<f64>::one().is_one());

        assert!(Poly::<f32>::zero().is_zero());
        assert!(Poly::<f32>::one().is_one());
    }

    #[test]
    fn integer_coeffs_identities() {
        assert!(Poly::<i8>::zero().is_zero());
        assert!(Poly::<i8>::one().is_one());

        assert!(Poly::<u8>::zero().is_zero());
        assert!(Poly::<u8>::one().is_one());

        assert!(Poly::<i16>::zero().is_zero());
        assert!(Poly::<i16>::one().is_one());

        assert!(Poly::<u16>::zero().is_zero());
        assert!(Poly::<u16>::one().is_one());

        assert!(Poly::<i32>::zero().is_zero());
        assert!(Poly::<i32>::one().is_one());

        assert!(Poly::<u32>::zero().is_zero());
        assert!(Poly::<u32>::one().is_one());

        assert!(Poly::<i64>::zero().is_zero());
        assert!(Poly::<i64>::one().is_one());

        assert!(Poly::<u64>::zero().is_zero());
        assert!(Poly::<u64>::one().is_one());

        assert!(Poly::<i128>::zero().is_zero());
        assert!(Poly::<i128>::one().is_one());

        assert!(Poly::<u128>::zero().is_zero());
        assert!(Poly::<u128>::one().is_one());

        assert!(Poly::<isize>::zero().is_zero());
        assert!(Poly::<isize>::one().is_one());

        assert!(Poly::<usize>::zero().is_zero());
        assert!(Poly::<usize>::one().is_one());
    }

    proptest! {
        #[test]
        fn qc_leading_coefficient(c: i8) {
            prop_assume!(c != 0);
            assert_eq!(c, poly!(1, -5, c).leading_coeff());
        }
    }

    #[test]
    fn monic_poly() {
        let p = poly!(-3., 6., 9.);
        let (p2, c) = p.monic();
        assert_relative_eq!(9., c);
        assert_relative_eq!(1., p2.leading_coeff());
    }

    #[test]
    fn monic_mutable_poly() {
        let mut p = poly!(-3., 6., 9.);
        let c = p.monic_mut();
        assert_relative_eq!(9., c);
        assert_relative_eq!(1., p.leading_coeff());
    }

    #[test]
    fn conversion_into_slice() {
        assert_eq!(&[3, -6, 8], poly!(3, -6, 8).as_ref());
    }

    #[test]
    fn round_off_coefficients() {
        let p = Poly::new_from_coeffs(&[1., 0.002, 1., -0.0001]);
        let actual = p.roundoff(&0.01);
        let expected = Poly::new_from_coeffs(&[1., 0., 1.]);
        assert_eq!(expected, actual);
    }

    #[test]
    fn round_off_zero() {
        let zero = Poly::zero();
        assert_eq!(zero, zero.roundoff(&0.001));
    }

    #[test]
    fn round_off_returns_zero() {
        let p = Poly::new_from_coeffs(&[0.0032, 0.002, -0.0023, -0.0001]);
        let actual = p.roundoff(&-0.01);
        assert_eq!(Poly::zero(), actual);
    }

    #[test]
    fn round_off_coefficients_mut() {
        let mut p = Poly::new_from_coeffs(&[1., 0.002, 1., -0.0001]);
        p.roundoff_mut(&0.01);
        let expected = Poly::new_from_coeffs(&[1., 0., 1.]);
        assert_eq!(expected, p);
    }
}

mod compile_fail_test {
    /// ```compile_fail
    /// use polynomen::{poly, Eval};
    /// let p = poly!(1.0e200, 2., 3.);
    /// p.eval(5.0_f32);
    /// ```
    #[allow(dead_code)]
    fn a() {}

    /// ``` compile_fail
    /// use polynomen::{poly, Eval};
    /// let p = poly!(1.5, 2., 3.);
    /// assert_eq!(86, p.eval(5));
    /// ```
    #[allow(dead_code)]
    fn b() {}
}