matext4cgmath/
lib.rs

1//! Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition.
2
3#![cfg_attr(not(debug_assertions), deny(warnings))]
4#![deny(clippy::all, rust_2018_idioms)]
5#![warn(
6    missing_docs,
7    missing_debug_implementations,
8    trivial_casts,
9    trivial_numeric_casts,
10    unsafe_code,
11    unstable_features,
12    unused_import_braces,
13    unused_qualifications
14)]
15#![allow(clippy::many_single_char_names)]
16
17pub use cgmath;
18use cgmath::*;
19use num_complex::Complex;
20
21mod eigens;
22mod exp_decomp;
23/// solvers for low dimensional algebraic equations.
24pub mod solver;
25
26/// extension for eigen values
27pub trait EigenValues: VectorSpace {
28    /// the type of the array of eigen values
29    type EigenValues;
30    /// calculate eigen values.
31    ///
32    /// # Examples
33    ///
34    /// ```
35    /// use cgmath::*;
36    /// use matext4cgmath::*;
37    /// use num_complex::Complex;
38    /// const EPS: f64 = 1.0e-10;
39    ///
40    /// let mat = Matrix2::new(4.0, -2.0, 3.0, -1.0);
41    /// let mut eigens = mat.eigenvalues();
42    /// // Even in the case of real solutions, the order in the array is not guaranteed.
43    /// eigens.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
44    /// assert!(Complex::norm(eigens[0] - 1.0) < EPS);
45    /// assert!(Complex::norm(eigens[1] - 2.0) < EPS);
46    /// ```
47    fn eigenvalues(self) -> Self::EigenValues;
48}
49
50#[cfg_attr(doc, katexit::katexit)]
51/// [operator norms](https://en.wikipedia.org/wiki/Matrix_norm): $L^1$, $L^2$, and $L^\infty$.
52/// # Examples
53///
54/// ```
55/// use cgmath::*;
56/// use matext4cgmath::*;
57///
58/// let mat = Matrix2::new(1.0, 3.0, -2.0, 4.0);
59/// // L^1 operator norm is the maximum column absolute sumation.
60/// assert_eq!(mat.norm_l1(), 6.0);
61/// // L^2 operator norm is the maximum singular value.
62/// let ans_norm_l2 = (5.0 + f64::sqrt(5.0)) / f64::sqrt(2.0);
63/// assert!(f64::abs(mat.norm_l2() - ans_norm_l2) < 1.0e-10);
64/// // L^∞ operator norm is the maximum row absolute sumation.
65/// assert_eq!(mat.norm_linf(), 7.0);
66/// ```
67pub trait OperatorNorm: VectorSpace {
68    /// operator norm for absolute value sumation: $L^1$.
69    fn norm_l1(self) -> Self::Scalar;
70    /// operator norm for Euclidean norm: $L^2$.
71    fn norm_l2(self) -> Self::Scalar;
72    /// operator norm for absolute value maximum: $L^{\infty}$.
73    fn norm_linf(self) -> Self::Scalar;
74}
75
76/// calculate exponential value
77pub trait Exponential: OperatorNorm + std::ops::AddAssign<Self> + One
78where
79    Self::Scalar: BaseFloat,
80{
81    /// calculate exponential
82    ///
83    /// # Examples
84    ///
85    /// ```
86    /// use cgmath::*;
87    /// use matext4cgmath::*;
88    ///
89    /// let x = Matrix2::new(0.0, 1.0, -1.0, 0.0);
90    /// let res = x.exp();
91    /// let ans = Matrix2::from_angle(Rad(1.0));
92    /// assert!((res - ans).norm_l1() < 1.0e-10);
93    /// ```
94    fn exp(self) -> Self {
95        use num_traits::{Float, NumCast};
96        let eps = <Self::Scalar as Float>::epsilon();
97        let mut a = self;
98        let mut res = Self::one();
99        for i in 2_u16..=64 {
100            res += a;
101            a = a * self / <Self::Scalar as NumCast>::from(i).unwrap();
102            if a.norm_linf() < eps * self.norm_linf() {
103                break;
104            }
105        }
106        res
107    }
108}
109
110#[cfg_attr(doc, katexit::katexit)]
111/// some decompositions of matrix
112pub trait Decomposition: VectorSpace {
113    #[cfg_attr(doc, katexit::katexit)]
114    /// Returns $(K, A, N)$ of the Iwasawa decomposition: $M = KAN$.
115    ///
116    /// - $K$: orthonormal matrix
117    /// - $A$: diagonal matrix
118    /// - $N$: upper-half unipotent matrix
119    fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)>;
120}