//! Rational functions
//!
//! Ratio whose numerator and denominator are polynomials.
//!
//! ```text
//! b_n*x^n + b_(n-1)*x^(n-1) + ... + b_1*x + b_0
//! f(x) = ---------------------------------------------
//! a_m*x^m + a_(m-1)*x^(m-1) + ... + a_1*x + a_0
//! ```
use std::{
fmt,
fmt::{Debug, Display, Formatter},
ops::{Add, Div, Mul},
};
use polynomen::Poly;
use crate::{
wrappers::{PWrapper, PP},
One, Polynomial, Zero,
};
mod arithmetic;
/// Rational function
#[derive(Clone, Debug, PartialEq)]
pub struct Rf<T> {
/// Rational function numerator
num: PP<T>,
/// Rational function denominator
den: PP<T>,
}
impl<T> Rf<T>
where
T: Clone + PartialEq + Zero,
{
/// Create a new rational function given its numerator and denominator
///
/// # Arguments
///
/// * `num` - Rational function numerator
/// * `den` - Rational function denominator
///
/// # Example
/// ```
/// use automatica::Rf;
/// let rf = Rf::new([1., 2.], [-4., 6., -2.]);
/// ```
#[must_use]
pub fn new<R, S>(num: R, den: S) -> Self
where
R: AsRef<[T]>,
S: AsRef<[T]>,
{
Self {
num: PP::new_from_coeffs(num.as_ref()),
den: PP::new_from_coeffs(den.as_ref()),
}
}
#[must_use]
pub(crate) fn new_pp(num: PP<T>, den: PP<T>) -> Self {
Self { num, den }
}
#[must_use]
pub(crate) fn new_from_poly(num: Poly<PWrapper<T>>, den: Poly<PWrapper<T>>) -> Self {
Self {
num: PP(num),
den: PP(den),
}
}
/// Extract rational function numerator
///
/// # Example
/// ```
/// use automatica::{Polynomial, Rf};
/// let num = [1., 2.];
/// let rf = Rf::new(num.clone(), [-4., 6., -2.]);
/// assert_eq!(num, rf.num().as_slice());
/// ```
#[must_use]
pub fn num(&self) -> impl Polynomial<T> {
self.num.to_unwrapped_vec()
}
#[must_use]
pub(crate) fn num_int(&self) -> &PP<T> {
&self.num
}
/// Extract rational function denominator
///
/// # Example
/// ```
/// use automatica::{Polynomial, Rf};
/// let den = [-4., 6., -2.];
/// let rf = Rf::new([1., 2.], den.clone());
/// assert_eq!(den, rf.den().as_slice());
/// ```
#[must_use]
pub fn den(&self) -> impl Polynomial<T> {
self.den.to_unwrapped_vec()
}
#[must_use]
pub(crate) fn den_int(&self) -> &PP<T> {
&self.den
}
}
impl<T: Clone + PartialEq + Zero> Rf<T> {
/// Calculate the relative degree between denominator and numerator.
///
/// # Example
/// ```
/// use automatica::{Inv, Rf};
/// let rf = Rf::new([1., 2.], [-4., 6., -2.]);
/// let expected = rf.relative_degree();
/// assert_eq!(expected, 1);
/// assert_eq!(rf.inv().relative_degree(), -1);
/// ```
#[must_use]
#[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)]
pub fn relative_degree(&self) -> i32 {
match (self.den.0.degree(), self.num.0.degree()) {
(Some(d), Some(n)) => d as i32 - n as i32,
(Some(d), None) => d as i32,
(None, Some(n)) => -(n as i32),
_ => 0,
}
}
}
macro_rules! rf_roots {
($ty:ty) => {
impl Rf<$ty> {
/// Calculate the poles of the rational function
#[must_use]
pub fn real_poles(&self) -> Option<Vec<$ty>> {
self.den.real_roots()
}
/// Calculate the poles of the rational function
#[must_use]
pub fn complex_poles(&self) -> Vec<($ty, $ty)> {
self.den.complex_roots()
}
/// Calculate the zeros of the rational function
#[must_use]
pub fn real_zeros(&self) -> Option<Vec<$ty>> {
self.num.real_roots()
}
/// Calculate the zeros of the rational function
#[must_use]
pub fn complex_zeros(&self) -> Vec<($ty, $ty)> {
self.num.complex_roots()
}
}
};
}
rf_roots!(f32);
rf_roots!(f64);
impl<T> Rf<T>
where
T: Clone + Div<Output = T> + One + PartialEq + Zero,
{
/// Normalization of rational function. If the denominator is zero the same
/// rational function is returned.
///
/// from:
/// ```text
/// b_n*z^n + b_(n-1)*z^(n-1) + ... + b_1*z + b_0
/// G(z) = ---------------------------------------------
/// a_n*z^n + a_(n-1)*z^(n-1) + ... + a_1*z + a_0
/// ```
/// to:
/// ```text
/// b'_n*z^n + b'_(n-1)*z^(n-1) + ... + b'_1*z + b'_0
/// G(z) = -------------------------------------------------
/// z^n + a'_(n-1)*z^(n-1) + ... + a'_1*z + a'_0
/// ```
///
/// # Example
/// ```
/// use automatica::Rf;
/// let rf = Rf::new([1., 2.], [-4., 6., -2.]);
/// let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
/// assert_eq!(expected, rf.normalize());
/// ```
#[must_use]
pub fn normalize(&self) -> Self {
if self.den.0.is_zero() {
return self.clone();
}
let (den, an) = self.den.monic();
let num = PP(&self.num.0 / an);
Self::new_pp(num, PP(den))
}
/// In place normalization of rational function. If the denominator is zero
/// no operation is done.
///
/// from:
/// ```text
/// b_n*z^n + b_(n-1)*z^(n-1) + ... + b_1*z + b_0
/// G(z) = ---------------------------------------------
/// a_n*z^n + a_(n-1)*z^(n-1) + ... + a_1*z + a_0
/// ```
/// to:
/// ```text
/// b'_n*z^n + b'_(n-1)*z^(n-1) + ... + b'_1*z + b'_0
/// G(z) = -------------------------------------------------
/// z^n + a'_(n-1)*z^(n-1) + ... + a'_1*z + a'_0
/// ```
///
/// # Example
/// ```
/// use automatica::Rf;
/// let mut rf = Rf::new([1., 2.], [-4., 6., -2.]);
/// rf.normalize_mut();
/// let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
/// assert_eq!(expected, rf);
/// ```
pub fn normalize_mut(&mut self) {
if self.den.0.is_zero() {
return;
}
let an = self.den.0.monic_mut();
self.num.0.div_mut(&an);
}
}
impl<T> Rf<T>
where
T: Clone + PartialEq,
{
/// Evaluate the rational function.
///
/// # Arguments
///
/// * `s` - Value at which the rational function is evaluated.
///
/// # Example
/// ```
/// use automatica::{Complex as C, Rf};
/// let rf = Rf::new([1., 2., 3.], [-4., -3., 1.]);
/// assert_eq!(-8.5, rf.eval_by_val(3.));
/// assert_eq!(C::new(0.64, -0.98), rf.eval_by_val(C::new(0., 2.0)));
/// ```
pub fn eval_by_val<N>(&self, s: N) -> N
where
N: Add<T, Output = N> + Clone + Div<Output = N> + Mul<Output = N> + Zero,
{
let res = self.num.0.eval_by_val(PWrapper(s.clone())) / self.den.0.eval_by_val(PWrapper(s));
res.0
}
}
impl<T> Rf<T>
where
T: Clone + PartialEq,
{
/// Evaluate the rational function.
///
/// # Arguments
///
/// * `s` - Value at which the rational function is evaluated.
///
/// # Example
/// ```
/// use automatica::{Complex as C, Rf};
/// let rf = Rf::new([1., 2., 3.], [-4., -3., 1.]);
/// assert_eq!(-8.5, rf.eval(&3.));
/// assert_eq!(C::new(0.64, -0.98), rf.eval(&C::new(0., 2.0)));
/// ```
pub fn eval<N>(&self, s: &N) -> N
where
N: Add<T, Output = N> + Clone + Div<Output = N> + Mul<N, Output = N> + Zero,
{
// let x = PWrapper(s.clone());
// let res = self.num.0.eval(&x) / self.den.0.eval(&x);
let res = self.num.0.eval_by_val(PWrapper(s.clone()))
/ self.den.0.eval_by_val(PWrapper(s.clone()));
res.0
}
}
/// Implementation of rational function printing
impl<T> Display for Rf<T>
where
T: Clone + Display + One + PartialEq + PartialOrd + Zero,
{
fn fmt(&self, f: &mut Formatter) -> fmt::Result {
let (s_num, s_den) = if let Some(precision) = f.precision() {
let num = format!("{poly:.prec$}", poly = self.num.0, prec = precision);
let den = format!("{poly:.prec$}", poly = self.den.0, prec = precision);
(num, den)
} else {
let num = format!("{}", self.num.0);
let den = format!("{}", self.den.0);
(num, den)
};
let length = s_num.len().max(s_den.len());
let dash = "\u{2500}".repeat(length);
write!(f, "{}\n{}\n{}", s_num, dash, s_den)
}
}
#[cfg(test)]
mod tests {
use crate::{Complex, Inv};
use super::*;
#[test]
#[allow(clippy::float_cmp)]
fn rational_function_creation() {
let num = [1., 2., 3.];
let den = [-4.2, -3.12, 0.0012];
let rf = Rf::new(num, den);
assert_eq!(num, rf.num().as_slice());
assert_eq!(den, rf.den().as_slice());
}
#[test]
fn relative_degree() {
let rf = Rf::new([1., 2.], [-4., 6., -2.]);
let expected = rf.relative_degree();
assert_eq!(expected, 1);
assert_eq!(rf.inv().relative_degree(), -1);
assert_eq!(-1, Rf::new([1., 1.], Poly::zero()).relative_degree());
assert_eq!(1, Rf::new(Poly::zero(), [1., 1.]).relative_degree());
assert_eq!(
0,
Rf::<f32>::new(Poly::zero(), Poly::zero()).relative_degree()
);
}
#[test]
fn evaluation() {
let rf = Rf::new([-0.75, 0.25], [0.75, 0.75, 1.]);
let res = rf.eval(&Complex::new(0., 0.9));
assert_abs_diff_eq!(0.429, res.re, epsilon = 0.001);
assert_abs_diff_eq!(1.073, res.im, epsilon = 0.001);
}
#[test]
fn evaluation_by_value() {
let rf = Rf::new([-0.75, 0.25], [0.75, 0.75, 1.]);
let res1 = rf.eval(&Complex::new(0., 0.9));
let res2 = rf.eval_by_val(Complex::new(0., 0.9));
assert_eq!(res1, res2);
}
#[test]
fn poles() {
let rf = Rf::new([1.0_f64], [6., -5., 1.]);
assert_eq!(Some(vec![2., 3.]), rf.real_poles());
}
#[test]
fn complex_poles() {
let rf = Rf::new([1.0_f32], [10., -6., 1.]);
assert_eq!(
vec![Complex::new(3., -1.), Complex::new(3., 1.)],
rf.complex_poles()
);
}
#[test]
fn zeros() {
let rf = Rf::new([1.0_f64], [6., -5., 1.]);
assert_eq!(None, rf.real_zeros());
}
#[test]
fn complex_zeros() {
let rf = Rf::new([3.25_f32, 3., 1.], [10., -3., 1.]);
assert_eq!(
vec![Complex::new(-1.5, -1.), Complex::new(-1.5, 1.)],
rf.complex_zeros()
);
}
#[test]
fn print() {
let rf = Rf::new(Poly::one(), Poly::new_from_roots(&[-1.]));
assert_eq!(
"1\n\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\n1 +1s",
format!("{}", rf)
);
let rf2 = Rf::new([1.123], [0.987, -1.321]);
assert_eq!(
"1.12\n\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\n0.99 -1.32s",
format!("{:.2}", rf2)
);
}
#[test]
fn normalization() {
let rf = Rf::new([1., 2.], [-4., 6., -2.]);
let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
assert_eq!(expected, rf.normalize());
let rf2 = Rf::new([1.], [0.]);
assert_eq!(rf2, rf2.normalize());
}
#[test]
fn rf_normalization_mutable() {
let mut rf = Rf::new([1., 2.], [-4., 6., -2.]);
rf.normalize_mut();
let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
assert_eq!(expected, rf);
let mut rf2 = Rf::new([1.], [0.]);
let rf3 = rf2.clone();
rf2.normalize_mut();
assert_eq!(rf2, rf3);
}
#[test]
fn eval_rational_function() {
let s_num = [-1., 1.];
let s_den = [0., 1.];
let s = Rf::new(s_num, s_den);
let p = Poly::new_from_coeffs(&[1., 2., 3.]);
let r = p.eval(&s);
let expected = Rf::new([3., -8., 6.], [0., 0., 1.]);
assert_eq!(expected, r);
}
}