use crate::{
algebra::abstr::{
cast,
cast::{FromPrimitive, NumCast, ToPrimitive},
AbelianGroup, AbelianGroupAdd, AbelianGroupMul, Addition, CommutativeRing,
Complex as ComplexT, Field, Group, GroupAdd, GroupMul, Identity, Loop, Magma, MagmaAdd,
MagmaMul, Monoid, MonoidAdd, MonoidMul, Multiplication, One, Quasigroup, Real, Ring,
Scalar, Semigroup, SemigroupAdd, SemigroupMul, Sign, Zero,
},
elementary::{Exponential, Hyperbolic, Power, Trigonometry},
algebra::abstr::AbsDiffEq,
};
use std::{
cmp::Ordering,
fmt,
fmt::Display,
ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign},
};
#[cfg(feature = "blaslapack")]
use crate::algebra::abstr::{Blas, Lapack};
#[macro_export]
macro_rules! Complex {
($r:expr, $i:expr) => {
Complex::new($r, $i);
};
}
#[macro_export]
macro_rules! Complex32 {
($r:expr, $i:expr) => {
(Complex::new($r, $i));
};
}
#[macro_export]
macro_rules! Complex64 {
($r:expr, $i:expr) => {
(Complex::new($r, $i));
};
}
pub type Complex32 = Complex<f32>;
pub type Complex64 = Complex<f64>;
#[derive(Debug, Clone, Copy)]
pub struct Complex<T>
where T: Real
{
re: T,
im: T,
}
impl<T> Complex<T> where T: Real
{
pub fn new(re: T, im: T) -> Complex<T>
{
Complex { re, im }
}
}
impl<T> Neg for Complex<T> where T: Real
{
type Output = Complex<T>;
fn neg(self: Self) -> Complex<T>
{
Complex { re: -self.re,
im: -self.im }
}
}
impl<T> PartialOrd for Complex<T> where T: Real
{
fn partial_cmp(self: &Self, other: &Self) -> Option<Ordering>
{
self.re.partial_cmp(&other.re)
}
}
impl<T> ComplexT for Complex<T> where T: Real
{
fn conj(self: Self) -> Complex<T>
{
Complex { re: self.re,
im: -self.im }
}
fn arg(self: Self) -> Self
{
Complex { re: self.im.arctan2(self.re),
im: T::zero() }
}
}
#[cfg(feature = "blaslapack")]
impl<T> Lapack for Complex<T> where T: Real
{
fn xgehrd(_n: i32,
_ilo: i32,
_ihi: i32,
_a: &mut [Self],
_lda: i32,
_tau: &mut [Self],
_work: &mut [Self],
_lwork: i32,
_info: &mut i32)
{
unimplemented!();
}
fn xgehrd_work_size(_n: i32,
_ilo: i32,
_ihi: i32,
_a: &mut [Self],
_lda: i32,
_tau: &mut [Self],
_info: &mut i32)
-> i32
{
unimplemented!();
}
fn xorghr(_n: i32,
_ilo: i32,
_ihi: i32,
_a: &mut [Self],
_lda: i32,
_tau: &[Self],
_work: &mut [Self],
_lwork: i32,
_info: &mut i32)
{
unimplemented!();
}
fn xorghr_work_size(_n: i32,
_ilo: i32,
_ihi: i32,
_a: &mut [Self],
_lda: i32,
_tau: &[Self],
_info: &mut i32)
-> i32
{
unimplemented!();
}
fn xgeev(_jobvl: u8,
_jobvr: u8,
_n: i32,
_a: &mut [Self],
_lda: i32,
_wr: &mut [Self],
_wi: &mut [Self],
_vl: &mut [Self],
_ldvl: i32,
_vr: &mut [Self],
_ldvr: i32,
_work: &mut [Self],
_lwork: i32,
_info: &mut i32)
{
unimplemented!();
}
fn xgeev_work_size(_jobvl: u8,
_jobvr: u8,
_n: i32,
_a: &mut [Self],
_lda: i32,
_wr: &mut [Self],
_wi: &mut [Self],
_vl: &mut [Self],
_ldvl: i32,
_vr: &mut [Self],
_ldvr: i32,
_info: &mut i32)
-> i32
{
unimplemented!();
}
fn xgetrf(_m: i32, _n: i32, _a: &mut [Self], _lda: i32, _ipiv: &mut [i32], _info: &mut i32)
{
unimplemented!();
}
fn xgeqrf(_m: i32,
_n: i32,
_a: &mut [Self],
_lda: i32,
_tau: &mut [Self],
_work: &mut [Self],
_lwork: i32,
_info: &mut i32)
{
unimplemented!();
}
fn xgeqrf_work_size(_m: i32,
_n: i32,
_a: &mut [Self],
_lda: i32,
_tau: &mut [Self],
_info: &mut i32)
-> i32
{
unimplemented!();
}
fn xorgqr(_m: i32,
_n: i32,
_k: i32,
_a: &mut [Self],
_lda: i32,
_tau: &mut [Self],
_work: &mut [Self],
_lwork: i32,
_info: &mut i32)
{
unimplemented!();
}
fn xorgqr_work_size(_m: i32,
_n: i32,
_k: i32,
_a: &mut [Self],
_lda: i32,
_tau: &mut [Self],
_info: &mut i32)
-> i32
{
unimplemented!();
}
fn xgetri(_n: i32,
_a: &mut [Self],
_lda: i32,
_ipiv: &mut [i32],
_work: &mut [Self],
_lwork: i32,
_info: &mut i32)
{
unimplemented!();
}
fn xgetri_work_size(_n: i32,
_a: &mut [Self],
_lda: i32,
_ipiv: &mut [i32],
_info: &mut i32)
-> i32
{
unimplemented!();
}
fn xpotrf(_uplo: char, _n: i32, _a: &mut [Self], _lda: i32, _info: &mut i32)
{
unimplemented!();
}
fn xgetrs(_n: i32,
_nrhs: i32,
_a: &mut [Self],
_lda: i32,
_ipiv: &mut [i32],
_b: &mut [Self],
_ldb: i32,
_info: &mut i32)
{
unimplemented!();
}
}
#[cfg(feature = "blaslapack")]
impl<T> Blas for Complex<T> where T: Real
{
fn xgemm(_transa: u8,
_transb: u8,
_m: i32,
_n: i32,
_k: i32,
_alpha: Self,
_a: &[Self],
_lda: i32,
_b: &[Self],
_ldb: i32,
_beta: Self,
_c: &mut [Self],
_ldc: i32)
{
unimplemented!();
}
fn xtrsm(_side: char,
_uplo: char,
_transa: char,
_diag: char,
_m: i32,
_n: i32,
_alpha: Self,
_a: &[Self],
_lda: i32,
_b: &mut [Self],
_ldb: i32)
{
unimplemented!();
}
fn xscal(_n: i32, _a: Self, _x: &mut [Self], _inc: i32)
{
unimplemented!();
}
fn xaxpy(_n: i32, _a: Self, _x: &[Self], _incx: i32, _y: &mut [Self], _incy: i32)
{
unimplemented!();
}
}
impl<T> Sign for Complex<T> where T: Real
{
fn sign(self: &Self) -> Self
{
unimplemented!()
}
fn abs(self: &Self) -> Self
{
let root: T = T::from_f64(2.0);
Complex { re: (self.re * self.re + self.im * self.im).root(root),
im: T::zero() }
}
fn is_positive(self: &Self) -> bool
{
unimplemented!();
}
fn is_negative(&self) -> bool
{
unimplemented!();
}
}
impl<T> Scalar for Complex<T> where T: Real
{
}
impl<T> PartialEq for Complex<T> where T: Real
{
fn eq<'a, 'b>(self: &'a Self, rhs: &Self) -> bool
{
if self.re == rhs.re && self.im == rhs.im
{
return true;
}
false
}
}
impl<T> Display for Complex<T> where T: Real
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result
{
write!(f, "{} + {}i", self.re, self.im)
}
}
impl<T> Zero for Complex<T> where T: Real
{
fn zero() -> Self
{
return Complex::new(T::zero(), T::zero());
}
}
impl<T> Add for Complex<T> where T: Real
{
type Output = Complex<T>;
fn add(self, rhs: Complex<T>) -> Complex<T>
{
Complex { re: self.re + rhs.re,
im: self.im + rhs.im }
}
}
impl<'a, 'b, T> Add<&'b Complex<T>> for &'a Complex<T>
where T: Real
{
type Output = Complex<T>;
fn add(self: Self, rhs: &'b Complex<T>) -> Self::Output
{
Complex { re: self.re + rhs.re,
im: self.im + rhs.im }
}
}
impl<T> AddAssign for Complex<T> where T: Real
{
fn add_assign<'a>(self: &'a mut Self, other: Self)
{
self.re += other.re;
self.im += other.im;
}
}
impl<T> One for Complex<T> where T: Real
{
fn one() -> Self
{
return Complex::new(T::one(), T::zero());
}
}
impl<T> Mul for Complex<T>
where T: Real
{
type Output = Complex<T>;
fn mul(self: Self, other: Self) -> Self::Output
{
Complex { re: self.re * other.re - self.im * other.im,
im: self.im * other.re + self.re * other.im }
}
}
impl<'a, 'b, T> Mul<&'b Complex<T>> for &'a Complex<T>
where T: Real
{
type Output = Complex<T>;
fn mul(self: Self, rhs: &'b Complex<T>) -> Self::Output
{
Complex { re: self.re * rhs.re - self.im * rhs.im,
im: self.im * rhs.re + self.re * rhs.im }
}
}
impl<T> MulAssign for Complex<T>
where T: Real
{
fn mul_assign(self: &mut Self, other: Self)
{
let temp: Complex<T> = *self;
self.re = temp.re * other.re - temp.im * other.im;
self.im = temp.im * other.re + temp.re * other.im;
}
}
impl<T> Sub for Complex<T> where T: Real
{
type Output = Complex<T>;
fn sub(self: Self, rhs: Self) -> Self::Output
{
&self - &rhs
}
}
impl<'a, 'b, T> Sub<&'b Complex<T>> for &'a Complex<T>
where T: Real
{
type Output = Complex<T>;
fn sub(self: Self, rhs: &'b Complex<T>) -> Self::Output
{
Complex { re: self.re - rhs.re,
im: self.im - rhs.im }
}
}
impl<T> SubAssign for Complex<T> where T: Real
{
fn sub_assign<'a>(self: &'a mut Self, other: Self)
{
self.re -= other.re;
self.im -= other.im;
}
}
impl<T> Div for Complex<T> where T: Real
{
type Output = Complex<T>;
fn div(self: Self, rhs: Self) -> Self::Output
{
&self / &rhs
}
}
impl<'a, 'b, T> Div<&'b Complex<T>> for &'a Complex<T> where T: Real
{
type Output = Complex<T>;
fn div(self: Self, rhs: &'b Complex<T>) -> Self::Output
{
let divisor: T = rhs.re * rhs.re + rhs.im * rhs.im;
if divisor == T::zero()
{
panic!("rhs * rhs.conj() == 0")
}
let quot: Complex<T> = self * &Complex::new(rhs.re, -rhs.im);
return Complex::new(quot.re / divisor, quot.im / divisor);
}
}
impl<T> DivAssign for Complex<T> where T: Real
{
fn div_assign<'a>(self: &'a mut Self, other: Self)
{
*self = *self / other;
}
}
macro_rules! impl_to_primitive {
($ty:ty, $to:ident) => {
fn $to(&self) -> $ty
{
return self.re.$to();
}
};
}
impl<T: ToPrimitive> ToPrimitive for Complex<T> where T: Real
{
impl_to_primitive!(u8, to_u8);
impl_to_primitive!(u16, to_u16);
impl_to_primitive!(u32, to_u32);
impl_to_primitive!(u64, to_u64);
impl_to_primitive!(u128, to_u128);
impl_to_primitive!(i8, to_i8);
impl_to_primitive!(i16, to_i16);
impl_to_primitive!(i32, to_i32);
impl_to_primitive!(i64, to_i64);
impl_to_primitive!(i128, to_i128);
impl_to_primitive!(f32, to_f32);
impl_to_primitive!(f64, to_f64);
}
impl<T> Exponential for Complex<T> where T: Real
{
fn e() -> Self
{
Complex { re: T::e(),
im: T::zero() }
}
fn exp(self: Self) -> Self
{
let k: T = self.re.exp();
Complex { re: k * self.im.cos(),
im: k * self.im.sin() }
}
fn ln(self: Self) -> Self
{
Complex { re: self.abs().re.ln(),
im: self.arg().re }
}
}
impl<T> Trigonometry for Complex<T> where T: Real
{
fn pi() -> Self
{
Complex { re: T::pi(),
im: T::zero() }
}
fn sin(self: Self) -> Self
{
let a: Complex<T> = Complex!(-self.im, self.re);
let b: Complex<T> = Complex!(self.im, -self.re);
let c: Complex<T> = Complex!(T::zero(), T::one() + T::one());
(a.exp() - b.exp()) / c
}
fn cos(self: Self) -> Self
{
let a: Complex<T> = Complex!(-self.im, self.re);
let b: Complex<T> = Complex!(self.im, -self.re);
let c: Complex<T> = Complex!(T::one() + T::one(), T::zero());
(a.exp() + b.exp()) / c
}
fn tan(self: Self) -> Self
{
return self.sin() / self.cos();
}
fn cot(self: Self) -> Self
{
Complex::one() / self.tan()
}
fn sec(self: Self) -> Self
{
Complex::one() / self.cos()
}
fn csc(self: Self) -> Self
{
Complex::one() / self.sin()
}
fn arcsin(self: Self) -> Self
{
let mi: Complex<T> = Complex!(T::zero(), -T::one());
let iz: Complex<T> = Complex!(-self.im, self.re);
let exp: Complex<T> = Complex!(T::one() / (T::one() + T::one()), T::zero());
mi * (iz + (Complex::one() - self * self).pow(exp)).ln()
}
fn arccos(self: Self) -> Self
{
Complex!(T::pi() / (T::one() + T::one()), T::zero()) - self.arcsin()
}
fn arctan(self: Self) -> Self
{
let two: T = T::one() + T::one();
let re: T;
if self.re == T::zero()
{
if self.im.abs() <= T::one()
{
re = T::zero()
}
else
{
if self.im > T::zero()
{
re = T::pi() / two;
}
else
{
re = -T::pi() / two;
}
}
}
else
{
if self.re > T::zero()
{
re =
(((self.re * self.re + self.im * self.im - T::one()) / (two * self.re)).arctan()
+ T::pi() / two)
/ two
}
else
{
re =
(((self.re * self.re + self.im * self.im - T::one()) / (two * self.re)).arctan()
- T::pi() / two)
/ two
}
}
let im: T =
((two * self.im) / (self.re * self.re + self.im * self.im + T::one())).artanh() / two;
Complex!(re, im)
}
fn arctan2(self: Self, _other: Self) -> Self
{
unimplemented!()
}
fn arccot(self: Self) -> Self
{
if self.re == T::zero()
{
if self.im == T::one() || self.im == -T::one()
{
panic!()
}
}
(Complex::one() / self).arctan()
}
fn arcsec(self: Self) -> Self
{
if self.im == T::zero()
{
if self.re == -T::one() || self.re == T::zero() || self.re == T::one()
{
panic!()
}
}
(Complex::one() / self).arccos()
}
fn arccsc(self: Self) -> Self
{
(Complex::one() / self).arcsin()
}
}
impl<T> Power for Complex<T> where T: Real
{
fn pow(self: Self, exp: Self) -> Self
{
let r: T = self.abs().re;
let phi: T = self.arg().re;
let k: T = r.pow(exp.re) * (-exp.im * phi).exp();
let theta: T = r.ln() * exp.im + exp.re * phi;
let re: T = k * theta.cos();
let im: T = k * theta.sin();
Complex!(re, im)
}
fn root(self: Self, _root: Self) -> Self
{
unimplemented!();
}
fn sqrt(self: Self) -> Self
{
unimplemented!();
}
}
impl<T> Hyperbolic for Complex<T> where T: Real
{
fn sinh(self: Self) -> Self
{
Complex!(T::zero(), -T::one()) * Complex!(-self.im, self.re).sin()
}
fn cosh(self: Self) -> Self
{
Complex!(-self.im, self.re).cos()
}
fn tanh(self: Self) -> Self
{
self.sinh() / self.cosh()
}
fn coth(self: Self) -> Self
{
self.cosh() / self.sinh()
}
fn sech(self: Self) -> Self
{
Complex!(-self.im, self.re).sec()
}
fn csch(self: Self) -> Self
{
Complex!(T::zero(), -T::one()) * Complex!(-self.im, self.re).csc()
}
fn arsinh(self: Self) -> Self
{
let p: Complex<T> = Complex!(T::one() / (T::one() + T::one()), T::zero());
(self + (self * self + Complex::one()).pow(p)).ln()
}
fn arcosh(self: Self) -> Self
{
let p: Complex<T> = Complex!(T::one() / (T::one() + T::one()), T::zero());
(self + (self * self - Complex::one()).pow(p)).ln()
}
fn artanh(self: Self) -> Self
{
let f: Complex<T> = Complex!(T::one() / (T::one() + T::one()), T::zero());
((Complex::one() + self) / (Complex::one() - self)).ln() * f
}
fn arcoth(self: Self) -> Self
{
let f: Complex<T> = Complex!(T::one() / (T::one() + T::one()), T::zero());
((self + Complex::one()) / (self - Complex::one())).ln() * f
}
fn arsech(self: Self) -> Self
{
(Complex::one() / self).arcosh()
}
fn arcsch(self: Self) -> Self
{
(Complex::one() / self).arsinh()
}
}
impl<T> FromPrimitive for Complex<T> where T: Real
{
fn from_i64(_n: i64) -> Self
{
unimplemented!();
}
fn from_i128(_n: i128) -> Self
{
unimplemented!();
}
fn from_u64(n: u64) -> Self
{
Complex { re: cast::cast(n),
im: T::zero() }
}
fn from_u128(n: u128) -> Self
{
Complex { re: cast::cast(n),
im: T::zero() }
}
fn from_f64(n: f64) -> Self
{
Complex { re: cast::cast(n),
im: T::zero() }
}
}
impl<T> NumCast for Complex<T> where T: Real
{
fn from<K: ToPrimitive>(n: K) -> Self
{
Complex { re: cast::cast(n.to_f64()),
im: T::zero() }
}
}
impl<T> Identity<Addition> for Complex<T> where T: Identity<Addition> + Real
{
fn id() -> Self
{
return Complex::new(Identity::<Addition>::id(), Identity::<Addition>::id());
}
}
impl<T> Identity<Multiplication> for Complex<T>
where T: Identity<Multiplication> + Identity<Addition> + Real
{
fn id() -> Self
{
return Complex::new(Identity::<Multiplication>::id(), Identity::<Addition>::id());
}
}
impl<T> Magma<Addition> for Complex<T> where T: Real
{
fn operate(self, rhs: Self) -> Self
{
return Complex::new(self.re + rhs.re, self.im + rhs.im);
}
}
impl<T> AbsDiffEq for Complex<T>
where T: Real
{
type Epsilon = Self;
fn default_epsilon() -> Self
{
unimplemented!()
}
fn abs_diff_eq(&self, _other: &Self, _epsilon: Self) -> bool
{
unimplemented!();
}
}
impl<T> MagmaAdd for Complex<T> where T: Real
{
}
impl<T> Magma<Multiplication> for Complex<T> where T: Real
{
fn operate(self, rhs: Self) -> Self
{
return Complex::new(self.re * rhs.re - self.im * rhs.im,
self.re * rhs.im + self.im * rhs.re);
}
}
impl<T> MagmaMul for Complex<T> where T: Real
{
}
impl<T> Quasigroup<Addition> for Complex<T> where T: Real
{
}
impl<T> Quasigroup<Multiplication> for Complex<T> where T: Real
{
}
impl<T> AbelianGroup<Addition> for Complex<T> where T: Real
{
}
impl<T> AbelianGroupAdd for Complex<T> where T: Real
{
}
impl<T> AbelianGroup<Multiplication> for Complex<T> where T: Real
{
}
impl<T> AbelianGroupMul for Complex<T> where T: Real
{
}
impl<T> Loop<Addition> for Complex<T> where T: Real
{
}
impl<T> Loop<Multiplication> for Complex<T> where T: Real
{
}
impl<T> CommutativeRing for Complex<T> where T: Real
{
}
impl<T> Ring for Complex<T> where T: Real
{
}
impl<T> Monoid<Addition> for Complex<T> where T: Real
{
}
impl<T> MonoidAdd for Complex<T> where T: Real
{
}
impl<T> Monoid<Multiplication> for Complex<T> where T: Real
{
}
impl<T> MonoidMul for Complex<T> where T: Real
{
}
impl<T> Semigroup<Addition> for Complex<T> where T: Real
{
}
impl<T> SemigroupAdd for Complex<T> where T: Real
{
}
impl<T> Semigroup<Multiplication> for Complex<T> where T: Real
{
}
impl<T> SemigroupMul for Complex<T> where T: Real
{
}
impl<T> Field for Complex<T> where T: Real
{
}
impl<T> Group<Addition> for Complex<T> where T: Real
{
}
impl<T> GroupAdd for Complex<T> where T: Real
{
}
impl<T> Group<Multiplication> for Complex<T> where T: Real
{
}
impl<T> GroupMul for Complex<T> where T: Real
{
}