automatica 1.0.0

//! Arithmetic module for rational functions

use std::ops::{Add, Div, Mul, Neg, Sub};

use crate::{
    rational_function::Rf,
    wrappers::{PWrapper, PP},
    Inv, One, Zero,
};

/// Implementation of rational function methods
impl<T: Clone + PartialEq> Rf<T> {
    /// Compute the reciprocal of a rational function in place.
    pub fn inv_mut(&mut self) {
        std::mem::swap(&mut self.num, &mut self.den);
    }
}

impl<T> Inv for &Rf<T>
where
    T: Clone + PartialEq + Zero,
{
    type Output = Rf<T>;

    /// Compute the reciprocal of a rational function.
    fn inv(self) -> Self::Output {
        Self::Output::new_pp(self.den.clone(), self.num.clone())
    }
}

impl<T: Clone + PartialEq> Inv for Rf<T> {
    type Output = Self;

    /// Compute the reciprocal of a rational function.
    fn inv(mut self) -> Self::Output {
        self.inv_mut();
        self
    }
}

/// Implementation of rational function negation.
/// Negative sign is transferred to the numerator.
impl<T> Neg for &Rf<T>
where
    T: Clone + Neg<Output = T> + PartialEq + Zero,
{
    type Output = Rf<T>;

    fn neg(self) -> Self::Output {
        Self::Output::new_pp(PP(-&self.num.0), self.den.clone())
    }
}

/// Implementation of rational function negation.
/// Negative sign is transferred to the numerator.
impl<T> Neg for Rf<T>
where
    T: Clone + Neg<Output = T> + PartialEq,
{
    type Output = Self;

    fn neg(mut self) -> Self::Output {
        self.num = PP(-self.num.0);
        self
    }
}

/// Implementation of rational function addition
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T> Add for &Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Rf<T>;

    fn add(self, rhs: Self) -> Self::Output {
        if self.is_zero() {
            return rhs.clone();
        }
        if rhs.is_zero() {
            return self.clone();
        }
        let (num, den) = if self.den == rhs.den {
            (&self.num.0 + &rhs.num.0, self.den.0.clone())
        } else {
            (
                &self.num.0 * &rhs.den.0 + &self.den.0 * &rhs.num.0,
                &self.den.0 * &rhs.den.0,
            )
        };
        Self::Output::new_pp(PP(num), PP(den))
    }
}

/// Implementation of rational function addition
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T> Add for Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Self;

    fn add(mut self, rhs: Self) -> Self {
        if self.is_zero() {
            return rhs;
        }
        if rhs.is_zero() {
            return self;
        }
        if self.den == rhs.den {
            self.num.0 = self.num.0 + rhs.num.0;
        } else {
            // first modify numerator.
            self.num.0 = &self.num.0 * &rhs.den.0 + &self.den.0 * &rhs.num.0;
            self.den.0 = self.den.0 * rhs.den.0;
        }
        self
    }
}

/// Implementation of rational function addition
impl<T> Add<T> for Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + PartialEq + Zero,
{
    type Output = Self;

    fn add(mut self, rhs: T) -> Self {
        self.num.0 = self.num.0 + &self.den.0 * PWrapper(rhs);
        self
    }
}

/// Implementation of rational function addition
impl<T> Add<&T> for Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + PartialEq + Zero,
{
    type Output = Self;

    fn add(mut self, rhs: &T) -> Self {
        self.num.0 = self.num.0 + &self.den.0 * PWrapper(rhs.clone());
        self
    }
}

/// Implementation of rational function subtraction
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T> Sub for &Rf<T>
where
    T: Add<Output = T>
        + Clone
        + Mul<Output = T>
        + Neg<Output = T>
        + PartialEq
        + Sub<Output = T>
        + Zero
        + One,
{
    type Output = Rf<T>;

    fn sub(self, rhs: Self) -> Self::Output {
        if self.is_zero() {
            return -rhs.clone();
        }
        if rhs.is_zero() {
            return self.clone();
        }
        let (num, den) = if self.den == rhs.den {
            (&self.num.0 - &rhs.num.0, self.den.0.clone())
        } else {
            (
                &self.num.0 * &rhs.den.0 - &self.den.0 * &rhs.num.0,
                &self.den.0 * &rhs.den.0,
            )
        };
        Self::Output::new_pp(PP(num), PP(den))
    }
}

/// Implementation of rational function subtraction
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T> Sub for Rf<T>
where
    T: Add<Output = T>
        + Clone
        + Mul<Output = T>
        + Neg<Output = T>
        + One
        + PartialEq
        + Sub<Output = T>
        + Zero,
{
    type Output = Self;

    fn sub(mut self, rhs: Self) -> Self {
        if self.is_zero() {
            return -rhs;
        }
        if rhs.is_zero() {
            return self;
        }
        if self.den == rhs.den {
            self.num.0 = self.num.0 - rhs.num.0;
        } else {
            // first modify numerator.
            self.num.0 = &self.num.0 * &rhs.den.0 - &self.den.0 * &rhs.num.0;
            self.den.0 = self.den.0 * rhs.den.0;
        }
        self
    }
}

/// Implementation of rational function multiplication
impl<T> Mul for &Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Rf<T>;

    fn mul(self, rhs: Self) -> Self::Output {
        if self.is_zero() || rhs.is_zero() {
            return Self::Output::zero();
        }
        let num = &self.num.0 * &rhs.num.0;
        let den = &self.den.0 * &rhs.den.0;
        Self::Output::new_pp(PP(num), PP(den))
    }
}

/// Implementation of rational function multiplication
impl<T> Mul for Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Self;

    fn mul(mut self, rhs: Self) -> Self {
        if self.is_zero() || rhs.is_zero() {
            return Self::Output::zero();
        }
        self.num.0 = self.num.0 * rhs.num.0;
        self.den.0 = self.den.0 * rhs.den.0;
        self
    }
}

/// Implementation of rational function multiplication
impl<T> Mul<&Rf<T>> for Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Self;

    fn mul(mut self, rhs: &Rf<T>) -> Self {
        if self.is_zero() || rhs.is_zero() {
            return Self::Output::zero();
        }
        self.num.0 = self.num.0 * &rhs.num.0;
        self.den.0 = self.den.0 * &rhs.den.0;
        self
    }
}

/// Implementation of rational function division
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T> Div for &Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Rf<T>;

    fn div(self, rhs: Self) -> Self::Output {
        match (self.is_zero(), rhs.is_zero()) {
            (true, false) => return Self::Output::zero(),
            (false, true) => return Self::Output::zero().inv(),
            _ => (),
        };
        let num = &self.num.0 * &rhs.den.0;
        let den = &self.den.0 * &rhs.num.0;
        Self::Output::new_pp(PP(num), PP(den))
    }
}

/// Implementation of rational function division
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T> Div for Rf<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
{
    type Output = Self;

    fn div(mut self, rhs: Self) -> Self {
        match (self.is_zero(), rhs.is_zero()) {
            (true, false) => return Self::Output::zero(),
            (false, true) => return Self::Output::zero().inv(),
            _ => (),
        };
        self.num.0 = self.num.0 * rhs.den.0;
        self.den.0 = self.den.0 * rhs.num.0;
        self
    }
}

impl<T> Zero for Rf<T>
where
    T: Clone + One + PartialEq + Zero,
{
    fn zero() -> Self {
        Self::new_pp(PP::zero(), PP::one())
    }

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

impl<T> polynomen::Zero for Rf<T>
where
    T: Clone + One + PartialEq + Zero,
{
    fn zero() -> Self {
        Self::new_pp(PP::zero(), PP::one())
    }

    fn is_zero(&self) -> bool {
        Zero::is_zero(&self.num.0) && !Zero::is_zero(&self.den.0)
    }
}

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

    #[test]
    fn rf_inversion() {
        let num1 = [1., 2., 3.];
        let den1 = [-4.2, -3.12, 0.0012];
        let rf1 = Rf::new(num1, den1);
        let num2 = [-4.2, -3.12, 0.0012];
        let den2 = [1., 2., 3.];
        let mut rf2 = Rf::new(num2, den2);

        assert_eq!(rf2, (&rf1).inv());
        rf2.inv_mut();
        assert_eq!(rf2, rf1);

        assert_eq!(rf2.inv(), rf1.inv());
    }

    #[test]
    fn rf_neg() {
        let rf1 = Rf::new([1., 2.], [1., 5.]);
        let rf2 = Rf::new([-1., -2.], [1., 5.]);
        assert_eq!(-&rf1, rf2);
        assert_eq!(rf1, -(-(&rf1)));
    }

    #[test]
    fn add_references() {
        let rf1 = Rf::new([1., 2.], [1., 5.]);
        let rf2 = Rf::new([3.], [1., 5.]);
        let actual = &rf1 + &rf2;
        let expected = Rf::new([4., 2.], [1., 5.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, &expected + &Rf::zero());
        assert_eq!(expected, &Rf::zero() + &expected);
    }

    #[test]
    fn add_values() {
        let rf1 = Rf::new([1., 2.], [3., -4.]);
        let rf2 = Rf::new([3.], [1., 5.]);
        let actual = rf1 + rf2;
        let expected = Rf::new([10., -5., 10.], [3., 11., -20.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, expected.clone() + Rf::zero());
        assert_eq!(expected, Rf::<f32>::zero() + expected.clone());
    }

    #[test]
    fn add_multiple_values() {
        let rf1 = Rf::new([1., 2.], [3., -4.]);
        let rf2 = Rf::new([3.], [1., 5.]);
        let rf3 = Rf::new([0., 4.], [3., 11., -20.]);
        let actual = &(&rf1 + &rf2) + &rf3;
        let expected = Rf::new([10., -1., 10.], [3., 11., -20.]);
        assert_eq!(expected, actual);

        let actual2 = (rf1 + rf2) + rf3;
        assert_eq!(actual, actual2);
    }

    #[test]
    fn add_scalar_value() {
        let rf = Rf::new([1., 2.], [3., -4.]);
        let actual = rf + 1.;
        let expected = Rf::new([4., -2.], [3., -4.]);
        assert_eq!(expected, actual);
    }

    #[test]
    #[allow(clippy::op_ref)]
    fn add_scalar_reference() {
        let rf = Rf::new([1., 2.], [3., -4.]);
        let actual = rf + &2.;
        let expected = Rf::new([7., -6.], [3., -4.]);
        assert_eq!(expected, actual);
    }

    #[test]
    fn sub_references() {
        let rf1 = Rf::new([-1., 9.], [4., -1.]);
        let rf2 = Rf::new([3.], [4., -1.]);
        let actual = &rf1 - &rf2;
        let expected = Rf::new([-4., 9.], [4., -1.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, &expected - &Rf::zero());
        assert_eq!(-&expected, &Rf::zero() - &expected);
    }

    #[test]
    fn sub_values() {
        let rf1 = Rf::new([1., -2.], [4., -4.]);
        let rf2 = Rf::new([2.], [2., 3.]);
        let actual = rf1 - rf2;
        let expected = Rf::new([-6., 7., -6.], [8., 4., -12.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, expected.clone() - Rf::zero());
        assert_eq!(-&expected, Rf::zero() - expected);
    }

    #[test]
    fn sub_multiple_values() {
        let rf1 = Rf::new([1., -2.], [4., -4.]);
        let rf2 = Rf::new([2.], [2., 3.]);
        let rf3 = Rf::new([0., 2.], [8., 4., -12.]);
        let actual = &(&rf1 - &rf2) - &rf3;
        let expected = Rf::new([-6., 5., -6.], [8., 4., -12.]);
        assert_eq!(expected, actual);

        let actual2 = (rf1 - rf2) - rf3;
        assert_eq!(actual, actual2);
    }

    #[test]
    fn mul_references() {
        let rf1 = Rf::new([1., 2., 3.], [1., 5.]);
        let rf2 = Rf::new([3.], [1., 6., 5.]);
        let actual = &rf1 * &rf2;
        let expected = Rf::new([3., 6., 9.], [1., 11., 35., 25.]);
        assert_eq!(expected, actual);
        assert_eq!(Rf::zero(), &expected * &Rf::zero());
        assert_eq!(Rf::zero(), &Rf::zero() * &expected);
    }

    #[test]
    fn mul_values() {
        let rf1 = Rf::new([1., 2., 3.], [1., 5.]);
        let rf2 = Rf::new([-5.], [1., 6., 5.]);
        let actual = rf1 * rf2;
        let expected = Rf::new([-5., -10., -15.], [1., 11., 35., 25.]);
        assert_eq!(expected, actual);
        assert_eq!(Rf::zero(), expected.clone() * Rf::zero());
        assert_eq!(Rf::zero(), Rf::zero() * expected);
    }

    #[test]
    fn mul_value_reference() {
        let rf1 = Rf::new([1., 2., 3.], [1., 5.]);
        let rf2 = Rf::new([-5.], [1., 6., 5.]);
        let actual = rf1 * &rf2;
        let expected = Rf::new([-5., -10., -15.], [1., 11., 35., 25.]);
        assert_eq!(expected, actual);
        assert_eq!(Rf::zero(), expected.clone() * &Rf::zero());
        assert_eq!(Rf::zero(), Rf::zero() * &expected);
    }

    #[test]
    fn div_values() {
        let rf1 = Rf::new([1.0_f32, 2., 3.], [1., 5.]);
        let rf2 = Rf::new([3.], [1., 6., 5.]);
        let actual = rf2 / rf1;
        let expected = Rf::new([3., 15.], [1., 8., 20., 28., 15.]);
        assert_eq!(expected, actual);
        assert_eq!(Rf::zero(), &Rf::zero() / &expected);
        assert!((&expected / &Rf::zero()).eval(&1.).is_infinite());
        assert!((&Rf::<f32>::zero() / &Rf::zero()).eval(&1.).is_nan());
    }

    #[test]
    #[allow(clippy::eq_op)]
    fn div_references() {
        let rf1 = Rf::new([1.0_f64, 2., 3.], [1., 5.]);
        let actual = &rf1 / &rf1;
        let expected = Rf::new([1., 7., 13., 15.], [1., 7., 13., 15.]);
        assert_eq!(expected, actual);
        assert_eq!(Rf::zero(), Rf::zero() / expected.clone());
        assert!((expected / Rf::zero()).eval(&1.).is_infinite());
        assert!((Rf::<f32>::zero() / Rf::zero()).eval(&1.).is_nan());
    }

    #[test]
    fn zero_rf() {
        assert!(Rf::<f32>::zero().is_zero());
        assert!(!Rf::new([0.], [0.]).is_zero());
    }

    #[test]
    fn zero_rf_polynomen() {
        assert!(<Rf::<f32> as polynomen::Zero>::zero().is_zero());
    }
}