automatica 1.0.0

//! Module that defines extension for complex numbers.

use std::{
    fmt::{Display, Formatter},
    ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub},
};

use crate::{
    traits::{Abs, Atan2, Cos, Exp, Hypot, Inv, One, Pow, Sin, Zero},
    wrappers::CWrapper,
};

/// Calculate the natural pulse of a complex number, it corresponds to its modulus.
///
/// # Arguments
///
/// * `c` - Complex number
///
/// # Example
/// ```
/// use automatica::{Complex, pulse};
/// let i = Complex::new(0., 1.);
/// assert_eq!(1., pulse(i));
/// ```
pub fn pulse<T: Hypot>(c: Complex<T>) -> T {
    c.re.hypot(c.im)
}

/// Calculate the damping of a complex number, it corresponds to the cosine of the
/// angle between the segment joining the complex number to the origin and the
/// real negative semiaxis.
///
/// By definition the damping of 0+0i is -1.
///
/// # Arguments
///
/// * `c` - Complex number
///
/// # Example
/// ```
/// use automatica::{Complex, damp};
/// let i = Complex::new(0., 1.);
/// assert_eq!(0., damp(i));
/// ```
pub fn damp<T>(c: Complex<T>) -> T
where
    T: Div<Output = T> + Neg<Output = T> + Hypot + One + Zero,
{
    let w = c.re.hypot(c.im);
    if w.is_zero() {
        // Handle the case where the pulse is zero to avoid division by zero.
        -T::one()
    } else {
        -c.re / w
    }
}

/// Complex number
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
pub struct Complex<T> {
    /// Real part
    pub re: T,
    /// Immaginary part
    pub im: T,
}

/// Convert tuple into complex number
impl<T> From<(T, T)> for Complex<T> {
    fn from((re, im): (T, T)) -> Self {
        Complex::new(re, im)
    }
}

/// Test equality between tuple and complex number
impl<T> PartialEq<(T, T)> for Complex<T>
where
    T: PartialEq,
{
    fn eq(&self, other: &(T, T)) -> bool {
        self.re == other.0 && self.im == other.1
    }
}

impl<T> Complex<T> {
    /// Create a complex number from its real and immaginary part
    ///
    /// # Arguments
    ///
    /// * `re` - Real part
    /// * `im` - Immaginary part
    pub fn new(re: T, im: T) -> Self {
        Self { re, im }
    }
}

impl<T> Complex<T>
where
    T: Clone + Cos + Mul<Output = T> + Sin,
{
    /// Create a complex number from polar form
    ///
    /// # Arguments
    ///
    /// * `rho` - Complex number modulus
    /// * `theta` - Complex number angle
    pub(crate) fn from_polar(rho: T, theta: T) -> Self {
        Self::new(rho.clone() * theta.cos(), rho * theta.sin())
    }
}

impl<T> Complex<T>
where
    T: Clone + Hypot,
{
    /// Norm of a complex number
    pub fn norm(&self) -> T {
        self.re.hypot(self.im.clone())
    }
}

impl<T> Complex<T>
where
    T: Atan2 + Clone,
{
    /// Angle of a complex number
    pub fn arg(&self) -> T {
        self.im.atan2(self.re.clone())
    }
}

impl<T> Pow<T> for Complex<T>
where
    T: Abs
        + Add<Output = T>
        + Atan2
        + Clone
        + Cos
        + Div<Output = T>
        + Hypot
        + Inv<Output = T>
        + Mul<Output = T>
        + Neg<Output = T>
        + One
        + PartialOrd
        + Pow<T>
        + Sin
        + Sub<Output = T>
        + Zero,
{
    fn powf(&self, exp: T) -> Self {
        let (rho, theta) = (self.norm(), self.arg());
        Self::from_polar(rho.powf(exp.clone()), theta * exp)
    }

    #[allow(clippy::cast_sign_loss)]
    fn powi(&self, exp: i32) -> Self {
        let (mut z, mut n) = if exp < 0 {
            (self.inv(), exp.wrapping_neg() as u32)
        } else {
            (self.clone(), exp as u32)
        };

        let mut y = Self::one();
        while n > 0 {
            if n % 2 == 1 {
                y = &y * &z;
            }
            z = &z * &z;
            n /= 2;
        }
        y
    }
}

impl<T> Exp for Complex<T>
where
    T: Clone + Cos + Exp + Mul<Output = T> + Sin,
{
    fn exp(&self) -> Self {
        Self::from_polar(self.re.exp(), self.im.clone())
    }
}

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

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

impl<T> One for Complex<T>
where
    T: One + Zero,
{
    fn one() -> Self {
        Self::new(T::one(), T::zero())
    }

    fn is_one(&self) -> bool {
        self.re.is_one() && self.im.is_zero()
    }
}

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

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

impl<T> AddAssign for Complex<T>
where
    T: Add<Output = T> + Clone,
{
    fn add_assign(&mut self, rhs: Self) {
        self.re = self.re.clone() + rhs.re;
        self.im = self.im.clone() + rhs.im;
    }
}

impl<'a, T> Add for &'a Complex<T>
where
    T: Add<Output = T> + Clone,
{
    type Output = Complex<T>;

    fn add(self, rhs: Self) -> Self::Output {
        Self::Output::new(
            self.re.clone() + rhs.re.clone(),
            self.im.clone() + rhs.im.clone(),
        )
    }
}

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

    fn add(self, rhs: T) -> Self::Output {
        Self::Output::new(self.re + rhs, self.im)
    }
}

impl<T> Add<&T> for Complex<T>
where
    T: Add<Output = T> + Clone,
{
    type Output = Self;

    fn add(self, rhs: &T) -> Self::Output {
        Self::Output::new(self.re + rhs.clone(), self.im)
    }
}

impl<T> Div for Complex<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,
{
    type Output = Self;

    fn div(self, rhs: Self) -> Self::Output {
        let (re, im) = complex_division::compdiv(
            CWrapper(self.re),
            CWrapper(self.im),
            CWrapper(rhs.re),
            CWrapper(rhs.im),
        );
        Complex::new(re.0, im.0)
    }
}

impl<T> Div<T> for Complex<T>
where
    T: Clone + Div<Output = T>,
{
    type Output = Self;

    fn div(self, rhs: T) -> Self::Output {
        Complex::new(self.re / rhs.clone(), self.im / rhs)
    }
}

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

    fn neg(self) -> Self::Output {
        Self::Output::new(-self.re, -self.im)
    }
}

impl<T> Mul for Complex<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + Sub<Output = T>,
{
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        Self::Output::new(
            self.re.clone() * rhs.re.clone() - self.im.clone() * rhs.im.clone(),
            self.re * rhs.im + self.im * rhs.re,
        )
    }
}

impl<T> Mul<&Self> for Complex<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + Sub<Output = T>,
{
    type Output = Self;

    fn mul(self, rhs: &Self) -> Self::Output {
        Self::Output::new(
            self.re.clone() * rhs.re.clone() - self.im.clone() * rhs.im.clone(),
            self.re * rhs.im.clone() + self.im * rhs.re.clone(),
        )
    }
}
impl<T> Mul for &Complex<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + Sub<Output = T>,
{
    type Output = Complex<T>;

    fn mul(self, rhs: Self) -> Self::Output {
        Self::Output::new(
            self.re.clone() * rhs.re.clone() - self.im.clone() * rhs.im.clone(),
            self.re.clone() * rhs.im.clone() + self.im.clone() * rhs.re.clone(),
        )
    }
}

impl<T> Mul<T> for Complex<T>
where
    T: Clone + Mul<Output = T>,
{
    type Output = Self;

    fn mul(self, rhs: T) -> Self::Output {
        Self::Output::new(self.re * rhs.clone(), self.im * rhs)
    }
}

impl<T> MulAssign for Complex<T>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + Sub<Output = T>,
{
    fn mul_assign(&mut self, rhs: Self) {
        let re = self.re.clone() * rhs.re.clone() - self.im.clone() * rhs.im.clone();
        let im = self.re.clone() * rhs.im + self.im.clone() * rhs.re;
        self.re = re;
        self.im = im;
    }
}

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

    fn sub(self, rhs: T) -> Self::Output {
        Self::Output::new(self.re - rhs, self.im)
    }
}

impl<T> Inv for &Complex<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,
{
    type Output = Complex<T>;
    fn inv(self) -> Self::Output {
        let (re, im) =
            complex_division::compinv(CWrapper(self.re.clone()), CWrapper(self.im.clone()));
        Complex::new(re.0, im.0)
    }
}

impl<T> Display for Complex<T>
where
    T: Display + PartialOrd + Zero,
{
    fn fmt(&self, f: &mut Formatter) -> std::fmt::Result {
        self.re.fmt(f)?;
        if !f.sign_plus() && self.im >= T::zero() {
            write!(f, " +")?;
        } else {
            write!(f, " ")?;
        }
        self.im.fmt(f)?;
        write!(f, "i")
    }
}

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

    use crate::Complex;

    #[test]
    fn pulse_damp() {
        let c = Complex::new(4., 3.);
        assert_relative_eq!(5., pulse(c));
        assert_relative_eq!(-0.8, damp(c));

        let i = Complex::new(0., 1.);
        assert_relative_eq!(1., pulse(i));
        assert_relative_eq!(0., damp(i));

        let zero = Complex::new(0., 0.);
        assert_relative_eq!(0., pulse(zero));
        assert_relative_eq!(-1., damp(zero));
    }

    #[test]
    fn tuple_conversion() {
        assert_eq!(Complex::new(1, 2), Complex::from((1, 2)));
    }

    #[test]
    fn float_power() {
        let c = Complex::new(-4., 0.).powf(0.5);
        assert_relative_eq!(0., c.re);
        assert_relative_eq!(2., c.im);
    }

    #[test]
    fn positive_integer_power() {
        assert_eq!(Complex::new(-0.25, 0.), Complex::new(0., 2.).powi(-2));
        assert_eq!(Complex::new(0., -64.), Complex::new(0., 4.).powi(3));
    }

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

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

    #[test]
    fn add_assign() {
        let mut c = Complex::new(3., -5.);
        c += Complex::new(10., -38.);
        assert_eq!(Complex::new(13., -43.), c);
    }

    #[test]
    fn mul_assign() {
        let mut c = Complex::new(3., 0.);
        c *= Complex::new(0., -2.);
        assert_eq!(Complex::new(0., -6.), c);
    }

    #[test]
    fn display() {
        assert_eq!("3 +4i", format!("{}", Complex::new(3., 4.)));
        assert_eq!("-5.00 -0.10i", format!("{:0.2}", Complex::new(-5., -0.1)));
    }
}