automatica 1.0.0

Automatic control systems library
Documentation
use std::fmt::{Display, Formatter};
use std::ops::{Add, Div, Mul, Neg, Sub};

use polynomen::Poly;

use crate::traits::{Abs, Inv, One, Pow, Sign, Sqrt, Zero};

/// Polynomial wrapper for use in dependency crate
/// It wraps a Poly so that the interface of this crate does not leak the
/// internal implementation.
#[derive(Clone, Debug, PartialEq)]
pub struct PP<T>(pub(crate) Poly<PWrapper<T>>);

impl<T> PP<T>
where
    T: Clone + PartialEq + Zero,
{
    /// Create a PP from a slice of values
    pub(crate) fn new_from_coeffs(coeffs: &[T]) -> Self {
        Self(Poly::new_from_coeffs_iter(
            coeffs.iter().cloned().map(PWrapper),
        ))
    }

    /// Create a PP from a slice of wrapped values
    #[allow(dead_code)]
    pub(crate) fn new_from_wrapped(coeffs: &[PWrapper<T>]) -> Self {
        Self(Poly::new_from_coeffs(coeffs))
    }
}

impl<T> PP<T>
where
    T: Clone,
{
    /// Get the vector of coefficients without any wrapping
    pub(crate) fn to_unwrapped_vec(&self) -> Vec<T> {
        self.0
            .as_slice()
            .iter()
            .cloned()
            .map(|x| x.0)
            .collect::<Vec<_>>()
    }
}

macro_rules! roots {
    ($ty:ty) => {
        impl PP<$ty> {
            /// Calculate real roots
            pub(crate) fn real_roots(&self) -> Option<Vec<$ty>> {
                let p = Poly::new_from_coeffs_iter(self.0.as_slice().into_iter().map(|x| x.0));
                p.real_roots()
            }

            /// Calculate complex roots
            pub(crate) fn complex_roots(&self) -> Vec<($ty, $ty)> {
                let p = Poly::new_from_coeffs_iter(self.0.as_slice().into_iter().map(|x| x.0));
                p.complex_roots()
            }
        }
    };
}
roots!(f32);
roots!(f64);

impl<T> PP<T>
where
    T: Clone + Div<Output = T> + One + PartialEq + Zero,
{
    /// Return the monic polynomial and its original leading coefficient
    pub(crate) fn monic(&self) -> (Poly<PWrapper<T>>, PWrapper<T>) {
        self.0.monic()
    }
}

impl<T> Zero for PP<T>
where
    T: Clone + PartialEq + Zero,
{
    fn zero() -> Self {
        use polynomen::Zero;
        PP(Poly::new_from_coeffs(&[PWrapper::zero()]))
    }

    fn is_zero(&self) -> bool {
        self.0.is_zero()
    }
}

impl<T> One for PP<T>
where
    T: Clone + PartialEq + One + Zero,
{
    fn one() -> Self {
        use polynomen::One;
        PP(Poly::new_from_coeffs(&[PWrapper::one()]))
    }

    fn is_one(&self) -> bool {
        self.0.is_one()
    }
}

impl<T> Display for PP<T>
where
    T: Display + PartialOrd + Zero,
{
    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
        self.0.fmt(f)
    }
}

/// Wrapper that abstracts the traits of the `polynomen` crate.
/// It allows the mapping between this crate traits and the `polynomen` ones.
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
pub(crate) struct PWrapper<T>(pub(crate) T);

impl<T> polynomen::Abs for PWrapper<T>
where
    T: Abs,
{
    fn abs(&self) -> Self {
        Self(self.0.abs())
    }
}

impl<T> polynomen::Inv for PWrapper<T>
where
    T: Clone + Inv<Output = T>,
{
    fn inv(&self) -> Self {
        Self(self.0.clone().inv())
    }
}

impl<T> polynomen::Pow for PWrapper<T>
where
    T: Pow<T>,
{
    fn powf(&self, exp: Self) -> Self {
        Self(self.0.powf(exp.0))
    }

    fn powi(&self, exp: i32) -> Self {
        Self(self.0.powi(exp))
    }
}

impl<T> polynomen::Sign for PWrapper<T>
where
    T: Sign,
{
    fn signum(&self) -> Self {
        Self(self.0.signum())
    }
    fn is_sign_negative(&self) -> bool {
        self.0.is_sign_negative()
    }
}

impl<T> polynomen::Sqrt for PWrapper<T>
where
    T: Sqrt,
{
    fn sqrt(&self) -> Self {
        Self(self.0.sqrt())
    }
}

impl<T> polynomen::One for PWrapper<T>
where
    T: One,
{
    fn one() -> Self {
        Self(T::one())
    }

    fn is_one(&self) -> bool {
        self.0.is_one()
    }
}

impl<T> polynomen::Zero for PWrapper<T>
where
    T: Zero,
{
    fn zero() -> Self {
        Self(T::zero())
    }

    fn is_zero(&self) -> bool {
        self.0.is_zero()
    }
}

impl<N, T> Add<PWrapper<T>> for PWrapper<N>
where
    N: Add<T, Output = N>,
{
    type Output = Self;

    fn add(self, rhs: PWrapper<T>) -> Self::Output {
        Self(self.0 + rhs.0)
    }
}

impl<'a, N, T> Add<&'a PWrapper<T>> for PWrapper<N>
where
    T: 'a,
    N: 'a + Add<&'a T, Output = N>,
{
    type Output = Self;

    fn add(self, rhs: &'a PWrapper<T>) -> Self::Output {
        Self(self.0 + &rhs.0)
    }
}

impl<T> Sub for PWrapper<T>
where
    T: Sub<Output = T>,
{
    type Output = Self;

    fn sub(self, rhs: Self) -> Self::Output {
        Self(self.0 - rhs.0)
    }
}

impl<T> Mul for PWrapper<T>
where
    T: Mul<Output = T>,
{
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        Self(self.0 * rhs.0)
    }
}

impl<'a, T> Mul<&'a PWrapper<T>> for PWrapper<T>
where
    T: 'a + Mul<&'a T, Output = T>,
{
    type Output = Self;

    fn mul(self, rhs: &'a PWrapper<T>) -> Self::Output {
        Self(self.0 * &rhs.0)
    }
}

impl<T> Div for PWrapper<T>
where
    T: Div<Output = T>,
{
    type Output = Self;

    fn div(self, rhs: Self) -> Self::Output {
        Self(self.0 / rhs.0)
    }
}

impl<T> Neg for PWrapper<T>
where
    T: Neg<Output = T>,
{
    type Output = Self;

    fn neg(self) -> Self::Output {
        Self(-self.0)
    }
}

impl<T> Display for PWrapper<T>
where
    T: Display,
{
    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
        self.0.fmt(f)
    }
}

/// Wrapper that abstracts the traits of the `complex-division` crate.
/// It allows the mapping between this crate traits and the `complex-division` ones.
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
pub(crate) struct CWrapper<T>(pub(crate) T);

impl<T> Neg for CWrapper<T>
where
    T: Neg<Output = T>,
{
    type Output = Self;

    fn neg(self) -> Self::Output {
        Self(-self.0)
    }
}

impl<T> Add for CWrapper<T>
where
    T: Add<Output = T>,
{
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        Self(self.0 + rhs.0)
    }
}

impl<T> Mul for CWrapper<T>
where
    T: Mul<Output = T>,
{
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        Self(self.0 * rhs.0)
    }
}

impl<T> Div for CWrapper<T>
where
    T: Div<Output = T>,
{
    type Output = Self;

    fn div(self, rhs: Self) -> Self::Output {
        Self(self.0 / rhs.0)
    }
}

impl<T> Sub for CWrapper<T>
where
    T: Sub<Output = T>,
{
    type Output = Self;

    fn sub(self, rhs: Self) -> Self::Output {
        Self(self.0 - rhs.0)
    }
}

/// Forward `Number` trait to the traits in this package.
impl<T> complex_division::Number for CWrapper<T>
where
    T: Abs
        + Add<Output = T>
        + Clone
        + Div<Output = T>
        + Inv<Output = T>
        + Mul<Output = T>
        + Neg<Output = T>
        + PartialOrd
        + Sub<Output = T>
        + Zero,
{
    fn abs(&self) -> Self {
        Self(self.0.abs())
    }

    fn inv(&self) -> Self {
        Self(self.0.clone().inv())
    }

    fn is_zero(&self) -> bool {
        self.0.is_zero()
    }
}

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

    #[test]
    fn new_wrapped() {
        let p = PP::new_from_wrapped(&[PWrapper(1.0_f32), PWrapper(2.)]);
        assert_eq!(Some(vec![-0.5]), p.real_roots());
    }

    #[test]
    fn pp_zero() {
        assert!(PP::<f32>::zero().is_zero());
    }

    #[test]
    fn pp_one() {
        assert!(PP::<f64>::one().is_one());
    }

    #[test]
    fn pwrapper_abs() {
        use polynomen::Abs;
        assert_eq!(PWrapper(4.), PWrapper(-4.).abs());
    }

    #[test]
    fn pwrapper_inv() {
        use polynomen::Inv;
        assert_eq!(PWrapper(-0.25), PWrapper(-4.).inv());
    }

    #[test]
    fn pwrapper_powf() {
        use polynomen::Pow;
        assert_eq!(PWrapper(81.), PWrapper(3.).powf(PWrapper(4.)));
    }

    #[test]
    fn pwrapper_zero() {
        use polynomen::Zero;
        assert!(PWrapper::<f32>::zero().is_zero());
    }

    #[test]
    fn pwrapper_one() {
        use polynomen::One;
        assert!(PWrapper::<f64>::one().is_one());
    }
}