automatica 1.0.0

use polynomen::Poly;

/// Fused multiplication and addition trait
pub trait MulAdd<T = Self> {
    /// Result of fused multiplication and addition
    type Output;
    /// Fused multiplication and addition with reduced loss of precision
    fn mul_add(self, mul: T, add: T) -> Self::Output;
}

/// Multiplicative inversion of a number trait
pub trait Inv {
    /// Result of multiplicative inversion
    type Output;
    /// Multiplicative inversion of a number
    fn inv(self) -> Self::Output;
}

/// Absolute value trait
pub trait Abs {
    /// Absolute value
    fn abs(&self) -> Self;
}

/// Numeric constants
pub trait Const {
    /// Greek pi
    fn pi() -> Self;
    /// Twice Greek pi
    fn tau() -> Self;
}

/// Additive identity trait
pub trait Zero {
    /// Additive identity
    fn zero() -> Self;
    /// Check if `self` is the additive identity
    fn is_zero(&self) -> bool;
}

/// Multiplicative identity trait
pub trait One {
    /// Multiplicative identity
    fn one() -> Self;
    /// Check if `self` is the multiplicative identity
    fn is_one(&self) -> bool;
}

/// Cosine of a number trait
pub trait Cos {
    /// Cosine of a number
    fn cos(&self) -> Self;
}

/// Sine of a number trait
pub trait Sin {
    /// Sine of a number
    fn sin(&self) -> Self;
}

/// Tangent of a number trait
pub trait Tan {
    /// Tangent of a number
    fn tan(&self) -> Self;
}

/// Exponent of a number trait
pub trait Pow<T> {
    /// Exponent of a number
    fn powf(&self, exp: T) -> Self;
    /// Exponent of a number with integer exponent
    fn powi(&self, exp: i32) -> Self;
}

/// Exponential trait
pub trait Exp {
    /// Exponential function
    fn exp(&self) -> Self;
}

/// Sign of a number trait
pub trait Sign {
    /// Return 1 for positive, 0 for zero, -1 for negative
    /// If the number has +0.0 and -0.0, it does not return 0
    fn signum(&self) -> Self;
    /// Return true if the number si negative
    fn is_sign_negative(&self) -> bool;
}

/// Square root of a number trait
pub trait Sqrt {
    /// Square root of a number
    fn sqrt(&self) -> Self;
}

/// Logarithm of a number trait
pub trait Log {
    //// Natural logarithm of a number
    //fn ln(&self) -> Self;
    /// Base 10 logarithm of a number
    fn log10(&self) -> Self;
}

/// Trait for rounding
pub trait Floor {
    /// Round the number to the higher integer less than the number
    fn floor(&self) -> Self;
}

/// Trait for hypotenuse
pub trait Hypot {
    /// Calculate the hypotenuse right angle triangle of catheti `x` and `y`.
    fn hypot(&self, y: Self) -> Self;
}

/// Trait for arctangent
pub trait Atan2 {
    /// Arctangent of catheti `self` and `other`
    fn atan2(&self, other: Self) -> Self;
}

/// Infinity representation
pub trait Infinity {
    /// Positive infinity
    fn infinity() -> Self;
}

/// Trait to handle degrees
pub trait Degree {
    /// Convert radians to degrees
    fn to_degrees(self) -> Self;
}

/// Type conversion
pub trait NumCast: Sized {
    /// Convert from type `T`
    fn from(n: usize) -> Option<Self>;
}

/// Maximum of two numbers trait
pub trait Max {
    /// Maximum of two numbers
    fn max(&self, other: &Self) -> Self;
}

/// Equality trait.
pub trait RelativeEq {
    /// Equality between numbers.
    fn relative_eq(&self, b: &Self, epsilon: &Self) -> bool;
}

macro_rules! impl_trait_float {
    ($t:ty, $id:ident, $z:expr, $o:expr) => {
        impl MulAdd<$t> for $t {
            type Output = Self;
            fn mul_add(self, mul: $t, add: $t) -> Self::Output {
                self.mul_add(mul, add)
            }
        }
        impl Abs for $t {
            fn abs(&self) -> Self {
                <$t>::abs(*self)
            }
        }
        impl Cos for $t {
            fn cos(&self) -> Self {
                <$t>::cos(*self)
            }
        }
        impl Exp for $t {
            fn exp(&self) -> Self {
                <$t>::exp(*self)
            }
        }
        impl Inv for $t {
            type Output = Self;
            fn inv(self) -> Self::Output {
                <$t>::recip(self)
            }
        }
        impl Log for $t {
            fn log10(&self) -> Self {
                <$t>::log10(*self)
            }
        }
        impl Max for $t {
            fn max(&self, other: &Self) -> Self {
                <$t>::max(*self, *other)
            }
        }
        impl Pow<$t> for $t {
            fn powf(&self, exp: $t) -> Self {
                <$t>::powf(*self, exp)
            }
            fn powi(&self, exp: i32) -> Self {
                <$t>::powi(*self, exp)
            }
        }
        impl Sign for $t {
            fn signum(&self) -> Self {
                <$t>::signum(*self)
            }
            fn is_sign_negative(&self) -> bool {
                <$t>::is_sign_negative(*self)
            }
        }
        impl Sin for $t {
            fn sin(&self) -> Self {
                <$t>::sin(*self)
            }
        }
        impl Tan for $t {
            fn tan(&self) -> Self {
                <$t>::tan(*self)
            }
        }
        impl Hypot for $t {
            fn hypot(&self, y: Self) -> Self {
                <$t>::hypot(*self, y)
            }
        }
        impl Atan2 for $t {
            fn atan2(&self, other: Self) -> Self {
                <$t>::atan2(*self, other)
            }
        }
        impl Degree for $t {
            fn to_degrees(self) -> Self {
                <$t>::to_degrees(self)
            }
        }
        impl Floor for $t {
            fn floor(&self) -> Self {
                <$t>::floor(*self)
            }
        }
        impl Sqrt for $t {
            fn sqrt(&self) -> Self {
                <$t>::sqrt(*self)
            }
        }
        impl Infinity for $t {
            fn infinity() -> Self {
                std::$id::INFINITY
            }
        }
        impl Zero for $t {
            fn zero() -> Self {
                $z
            }
            fn is_zero(&self) -> bool {
                *self == $z
            }
        }
        impl One for $t {
            fn one() -> Self {
                $o
            }
            #[allow(clippy::float_cmp)]
            fn is_one(&self) -> bool {
                *self == $o
            }
        }
        impl Const for $t {
            fn pi() -> Self {
                std::$id::consts::PI
            }
            fn tau() -> Self {
                std::$id::consts::TAU
            }
        }
        impl NumCast for $t {
            fn from(n: usize) -> Option<Self> {
                if n == 0 {
                    return Some(0.);
                }
                // Safe cast $source is either 4 or 8 bytes.
                let size: u32 = std::mem::size_of::<usize>() as u32 * 8;
                let lz = n.leading_zeros();
                let tz = n.trailing_zeros();
                // Leading bit is not used for sign in unsigned type.
                // Since the first one of the mantissa is implicit, the actual
                // number of representable bits is #mantissa + 1.
                // Integer types cannot be converted in subnormal floats.
                if size - (lz + tz) > <$t>::MANTISSA_DIGITS + 1 {
                    None
                } else {
                    Some(n as $t)
                }
            }
        }
        impl RelativeEq for $t {
            fn relative_eq(&self, b: &Self, epsilon: &Self) -> bool {
                // https://floating-point-gui.de/errors/comparison/
                let a = self;
                let abs_a = a.abs();
                let abs_b = b.abs();
                let diff = (a - b).abs();
                if a == b {
                    // Handle exact equality.
                    true
                } else if a == &0. || b == &0. || (abs_a + abs_b < <$t>::MIN_POSITIVE) {
                    // Handle case where at least one number is close to zero.
                    diff < (epsilon * <$t>::MIN_POSITIVE)
                } else {
                    diff / (abs_a + abs_b).min(<$t>::MAX) < *epsilon
                }
            }
        }
    };
}
impl_trait_float!(f32, f32, 0.0, 1.0);
impl_trait_float!(f64, f64, 0.0, 1.0);

/// Trait for generic results that require a polynomial
pub trait Polynomial<T> {
    /// Get the copy of the coefficients as a vector
    fn coeffs(&self) -> Vec<T>;
    /// Create a `Self` with the give coefficients as a slice
    fn new_from_coeffs(coeffs: &[T]) -> Self;
    /// Get the coefficients as a slice
    fn as_slice(&self) -> &[T];
}

impl<T> Polynomial<T> for Poly<T>
where
    T: Clone + PartialEq + polynomen::Zero,
{
    fn coeffs(&self) -> Vec<T> {
        self.coeffs()
    }

    fn new_from_coeffs(coeffs: &[T]) -> Self {
        Poly::new_from_coeffs(coeffs)
    }

    fn as_slice(&self) -> &[T] {
        self.as_slice()
    }
}

impl<T> Polynomial<T> for Vec<T>
where
    T: Clone,
{
    fn coeffs(&self) -> Vec<T> {
        self.clone()
    }
    fn new_from_coeffs(coeffs: &[T]) -> Self {
        Vec::from(coeffs)
    }
    fn as_slice(&self) -> &[T] {
        self.as_slice()
    }
}

impl<T> Zero for Poly<T>
where
    T: PartialEq + polynomen::Zero,
{
    fn zero() -> Self {
        <Self as polynomen::Zero>::zero()
    }

    fn is_zero(&self) -> bool {
        <Self as polynomen::Zero>::is_zero(self)
    }
}

impl<T> One for Poly<T>
where
    T: PartialEq + polynomen::One,
{
    fn one() -> Self {
        <Self as polynomen::One>::one()
    }

    fn is_one(&self) -> bool {
        <Self as polynomen::One>::is_one(self)
    }
}

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

    #[test]
    fn numeric_numcast() {
        let a: f32 = NumCast::from(12).unwrap();
        assert_eq!(12., a);
        let a: f64 = NumCast::from(123_432).unwrap();
        assert_eq!(123_432., a);
        let a: Option<f32> = NumCast::from(8_000_000_001);
        assert_eq!(None, a);
    }

    #[test]
    fn numeric_numcast_f32_limit() {
        // f32 mantissa is 24 bit long, so anything longer than 25 bits of
        // significant bits is not representable as f32.
        let a: f32 = NumCast::from(0b1000_0000_0000_0000_0000_0000_1).unwrap();
        assert_eq!(1. + 24.0_f32.exp2(), a);
        let a: Option<f32> = NumCast::from(0b1000_0000_0000_0000_0000_0000_01);
        assert_eq!(None, a);
    }

    #[test]
    fn numeric_numcast_f64_limit() {
        // f64 mantissa is 53 bit long, so anything longer than 54 bits of
        // significant bits is not representable as f64.
        let a: f64 =
            NumCast::from(0b1000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_01)
                .unwrap();
        assert_eq!(1. + 53.0_f64.exp2(), a);
        let a: Option<f64> =
            NumCast::from(0b1000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_001);
        assert_eq!(None, a);
    }

    #[test]
    fn polynomial_trait_for_poly() {
        let p = <Poly<i32> as Polynomial<i32>>::new_from_coeffs(&[0, 1, 2]);
        let c = Polynomial::<i32>::coeffs(&p);
        let s = Polynomial::<i32>::as_slice(&p);
        assert_eq!(s, &c);
    }

    #[test]
    fn polynomial_trait_for_vec() {
        let p = <Vec<i32> as Polynomial<i32>>::new_from_coeffs(&[0, 1, 2]);
        let c = Polynomial::<i32>::coeffs(&p);
        let s = Polynomial::<i32>::as_slice(&p);
        assert_eq!(s, &c);
    }

    #[test]
    fn polynomial_trait_one() {
        assert!(<Poly<u8> as One>::one().is_one());
    }

    #[test]
    fn polynomial_trait_zero() {
        assert!(<Poly<u8> as Zero>::zero().is_zero());
    }
}