complex-division 1.0.1

//! # Library for the calculation of complex number division.
//!
//! This library implements the division between complex numbers with floating
//! point representation.
//!
//! The algorithm is aimed to avoid overflows and underflows.
//!
//! A dedicated method is available for complex number inversion.
//! It is guaranteed that:
//!
//! ```text
//! compdiv(1., 0., c, d) = compinv(c, d)
//! ```
//!
//! ## Minimum compiler version
//!
//! The minimum `rustc` compiler version is 1.56.
//!
//! ## References
//!
//! Michael Baudin, Robert L. Smith, A Robust Complex Division in Scilab, 2012, arXiv:1210.4539v2 [cs.MS]

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

/// Division between complex numbers that avoids overflows.
///
/// # Arguments
///
/// * `a` - Dividend real part
/// * `b` - Dividend imaginary part
/// * `c` - Divisor real part
/// * `d` - Divisor imaginary part
///
/// Michael Baudin, Robert L. Smith, A Robust Complex Division in Scilab, 2012, arXiv:1210.4539v2 [cs.MS]
pub fn compdiv<T: Number>(a: T, b: T, c: T, d: T) -> (T, T) {
    compdiv_impl(a, b, c, d)
}

/// Division between complex numbers that avoids overflows.
///
/// # Arguments
///
/// * `a` - Dividend real part
/// * `b` - Dividend imaginary part
/// * `c` - Divisor real part
/// * `d` - Divisor imaginary part
///
/// Michael Baudin, Robert L. Smith, A Robust Complex Division in Scilab, 2012, arXiv:1210.4539v2 [cs.MS]
#[allow(clippy::many_single_char_names)]
fn compdiv_impl<T: Number>(a: T, b: T, c: T, d: T) -> (T, T) {
    if d.abs() <= c.abs() {
        let (e, f) = compdiv_internal(a, b, c, d);
        (e, f)
    } else {
        // Real and imaginary parts shall be swapped.
        let (e, f) = compdiv_internal(b, a, d, c);
        // And the imaginary part shall change sign.
        (e, -f)
    }
}

/// Implementation of division between complex numbers.
///
/// # Arguments
///
/// * `a` - Dividend real part if Im{divisor} <= Re{divisor} else imaginary
/// * `b` - Dividend imaginary part if Im{divisor} <= Re{divisor} else real
/// * `c` - Divisor real part if Im{divisor} <= Re{divisor} else imaginary
/// * `d` - Divisor imaginary part if Im{divisor} <= Re{divisor} else real
///
/// Michael Baudin, Robert L. Smith, A Robust Complex Division in Scilab, 2012, arXiv:1210.4539v2 [cs.MS]
#[allow(clippy::many_single_char_names)]
fn compdiv_internal<T: Number>(a: T, b: T, c: T, d: T) -> (T, T) {
    let r = d.clone() / c.clone();
    let t = (c.clone() + d.clone() * r.clone()).inv();
    if r.is_zero() {
        let e = (a.clone() + d.clone() * (b.clone() / c.clone())) * t.clone();
        let f = (b - d * (a / c)) * t;
        (e, f)
    } else {
        let e = (a.clone() + b.clone() * r.clone()) * t.clone();
        let f = (b - a * r) * t;
        (e, f)
    }
}

/// Inversion of complex numbers that avoids overflows.
///
/// # Arguments
///
/// * `c` - Complex number real part
/// * `d` - Complex number imaginary part
///
/// It is implemented using the division algorithm, replacing the numerator with
/// 1+0i and simplify terms.
/// Michael Baudin, Robert L. Smith, A Robust Complex Division in Scilab, 2012, arXiv:1210.4539v2 [cs.MS]
pub fn compinv<T: Number>(c: T, d: T) -> (T, T) {
    compinv_impl(c, d)
}

/// Inversion of complex numbers that avoids overflows.
///
/// # Arguments
///
/// * `c` - Complex number real part
/// * `d` - Complex number imaginary part
///
/// It is implemented using the division algorithm, replacing the numerator with
/// 1+0i and simplify terms.
/// Michael Baudin, Robert L. Smith, A Robust Complex Division in Scilab, 2012, arXiv:1210.4539v2 [cs.MS]
#[allow(clippy::many_single_char_names)]
fn compinv_impl<T: Number>(c: T, d: T) -> (T, T) {
    if d.abs() <= c.abs() {
        let r = d.clone() / c.clone();
        let t = (c + d * r.clone()).inv();
        let e = t.clone();
        let f = -r * t;
        (e, f)
    } else {
        let r = c.clone() / d.clone();
        let t = (c * r.clone() + d).inv();
        let e = r * t.clone();
        let f = -t;
        (e, f)
    }
}

pub trait Number:
    PartialOrd
    + Clone
    + Add<Output = Self>
    + Sub<Output = Self>
    + Div<Output = Self>
    + Mul<Output = Self>
    + Neg<Output = Self>
{
    /// Returns the absolute value of the input.
    /// `x.abs() := |a|`
    fn abs(&self) -> Self;

    /// Returns the reciprocal of the input.
    /// `x.inv() := 1/x`
    fn inv(&self) -> Self;

    /// Returns `true` if the input is zero.
    /// `x.is_zero() := x=0`
    fn is_zero(&self) -> bool;
}

macro_rules! impl_trait {
    ($t:ty) => {
        impl Number for $t {
            fn abs(&self) -> Self {
                <$t>::abs(*self)
            }

            fn inv(&self) -> Self {
                <$t>::recip(*self)
            }

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

impl_trait!(f32);
impl_trait!(f64);

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

    use super::*;

    fn p2(n: i32) -> f64 {
        2.0_f64.powf(f64::from(n))
    }

    #[test]
    fn complex_division_a() {
        let a = compdiv(1., 1., 1., 1e307);
        assert_eq!(
            (1.000_000_000_000_000_1e-307, -1.000_000_000_000_000_1e-307),
            a
        );
    }

    #[test]
    fn complex_division_b() {
        let b = compdiv(1., 1., 1e-307, 1e-307);
        assert_eq!((1.000_000_000_000_000_1e307, 0.), b);
    }

    #[test]
    fn complex_division_c() {
        let c = compdiv(1e307, 1e-307, 1e204, 1e-204);
        assert_eq!((1e103, -1e-305), c);
    }

    #[test]
    fn complex_division_1() {
        let d1 = compdiv(1., 1., 1., p2(1023));
        assert_eq!((p2(-1023), -p2(-1023)), d1);
    }

    #[test]
    fn complex_division_2() {
        let d2 = compdiv(1., 1., p2(-1023), p2(-1023));
        assert_eq!((p2(1023), 0.), d2);
    }

    #[test]
    fn complex_division_3() {
        let d3 = compdiv(p2(1023), p2(-1023), p2(677), p2(-677));
        assert_eq!((p2(346), -p2(-1008)), d3);
    }

    #[test]
    fn complex_division_5() {
        let d5 = compdiv(p2(1020), p2(-844), p2(656), p2(-780));
        assert_eq!((p2(364), -p2(-1072)), d5);
    }

    #[test]
    fn complex_division_6() {
        let d6 = compdiv(p2(-71), p2(1021), p2(1001), p2(-323));
        assert_eq!((p2(-1072), p2(20)), d6);
    }

    /// Test if both elements of the tuple are NaN.
    fn is_nan((a, b): (f32, f32)) -> bool {
        a.is_nan() && b.is_nan()
    }

    #[test]
    fn complex_division_limits() {
        let c1 = (1., 1.);
        let c2 = (1., -1.);
        let c3 = (-1., 1.);
        let c4 = (-1., -1.);

        let zero = (0.0_f32, 0.0_f32);

        assert!(is_nan(compdiv(c1.0, c1.1, zero.0, zero.1)));
        assert!(is_nan(compdiv(c2.0, c2.1, zero.0, zero.1)));
        assert!(is_nan(compdiv(c3.0, c3.1, zero.0, zero.1)));
        assert!(is_nan(compdiv(c4.0, c4.1, zero.0, zero.1)));
    }

    #[test]
    fn complex_division_and_inversion() {
        // Number are chosen to keep precision during operations
        let a = compdiv(1., 8., 1., -16.);
        let b = compdiv(1., -16., 1., 8.);
        assert_eq!(a, compinv(b.0, b.1));
    }

    proptest! {
        #[test]
        fn qc_complex_inversion_algorithm(a: f64, b: f64) {
            prop_assume!(a != 0. && b != 0.);
            assert_eq!(compdiv(1., 0., a, b), compinv(a, b));
        }
    }

    #[test]
    fn complex_inversion_limits() {
        let c = compinv(0.0_f32, 0.0_f32);
        assert!(c.0.is_nan());
        assert!(c.1.is_nan());
    }

    #[test]
    fn complex_inversion_real() {
        let n = 1234.45;
        let inv = compinv(n, 0.);
        assert_eq!((n.inv(), 0.), inv);
    }

    #[test]
    fn complex_inversion_imaginary() {
        let n = 1234.45;
        let inv = compinv(0., n);
        assert_eq!((0., -n.inv()), inv);
    }

    #[test]
    fn complex_double_inversion() {
        let expected = (1., -13.);
        let a = compinv(expected.0, expected.1);
        let orig = compinv(a.0, a.1);
        assert_eq!(expected, orig);
    }
}