logicaffeine_base/
matrix_exp.rs

1//! The exact matrix exponential `e^{t·M}` of a rational matrix.
2//!
3//! The characteristic polynomial `χ` of `M` annihilates it, so `e^{t·M} = r(M)` for the polynomial
4//! `r` of degree below `n` that agrees with `e^{t·x}` on the roots of `χ` to their multiplicities:
5//! `r^{(j)}(λ) = t^j·e^{t·λ}` for `j < m(λ)` (Hermite interpolation). Writing `r = Σ c_k·x^k`, those
6//! conditions are a confluent Vandermonde system `V·c = (t^j·e^{t·λ})`, so `c = V⁻¹·(t^j·e^{t·λ})`
7//! with `V⁻¹` exact over the complex algebraic numbers, and each entry of `e^{t·M}` is a sum
8//! `Σ_λ β_λ·e^{t·λ}`. A conjugate pair of eigenvalues `a ± b·i` contributes
9//! `2·e^{t·a}·(Re β·cos(t·b) − Im β·sin(t·b))`, so every entry is a real number in the exponential
10//! exact form: `e^{[[0, 1], [−1, 0]]}` is `[[cos(1), sin(1)], [−sin(1), cos(1)]]` exactly.
11
12use crate::algebraic::Algebraic;
13use crate::complex_algebraic::{complex_roots, ComplexAlgebraic};
14use crate::exact_real::ExactReal;
15use crate::numeric::Rational;
16use crate::poly::Poly;
17
18/// `det(x·I − M)` by Faddeev–LeVerrier; `None` for a matrix that is not square.
19pub fn characteristic_polynomial(m: &[Vec<Rational>]) -> Option<Poly> {
20    let n = m.len();
21    if n == 0 || m.iter().any(|r| r.len() != n) {
22        return None;
23    }
24    let mul = |a: &[Vec<Rational>], b: &[Vec<Rational>]| -> Vec<Vec<Rational>> {
25        (0..n)
26            .map(|i| {
27                (0..n)
28                    .map(|j| (0..n).fold(Rational::zero(), |acc, k| acc.add(&a[i][k].mul(&b[k][j]))))
29                    .collect()
30            })
31            .collect()
32    };
33    let mut c = vec![Rational::zero(); n + 1];
34    c[n] = Rational::one();
35    let mut acc = vec![vec![Rational::zero(); n]; n];
36    for k in 1..=n {
37        acc = mul(m, &acc);
38        for (i, row) in acc.iter_mut().enumerate() {
39            row[i] = row[i].add(&c[n - k + 1]);
40        }
41        let am = mul(m, &acc);
42        let trace = (0..n).fold(Rational::zero(), |t, i| t.add(&am[i][i]));
43        c[n - k] = trace
44            .negated()
45            .div(&Rational::from_i64(k as i64))
46            .expect("k ≥ 1");
47    }
48    Some(Poly::from_coeffs(c))
49}
50
51/// `k!/(k − j)!`, the `j`-th derivative's coefficient of `x^k`; `0` for `j > k`.
52fn falling(k: usize, j: usize) -> i64 {
53    if j > k {
54        return 0;
55    }
56    ((k - j + 1)..=k).map(|v| v as i64).product::<i64>().max(1)
57}
58
59/// The inverse of a square matrix over the complex algebraic numbers; `None` when singular.
60fn invert(mut a: Vec<Vec<ComplexAlgebraic>>) -> Option<Vec<Vec<ComplexAlgebraic>>> {
61    let n = a.len();
62    let mut inv: Vec<Vec<ComplexAlgebraic>> = (0..n)
63        .map(|i| {
64            (0..n)
65                .map(|j| {
66                    if i == j {
67                        ComplexAlgebraic::one()
68                    } else {
69                        ComplexAlgebraic::zero()
70                    }
71                })
72                .collect()
73        })
74        .collect();
75    for col in 0..n {
76        let pivot = (col..n).find(|&r| !a[r][col].is_zero())?;
77        a.swap(col, pivot);
78        inv.swap(col, pivot);
79        let scale = a[col][col].recip()?;
80        for j in 0..n {
81            a[col][j] = a[col][j].mul(&scale);
82            inv[col][j] = inv[col][j].mul(&scale);
83        }
84        for r in 0..n {
85            if r == col || a[r][col].is_zero() {
86                continue;
87            }
88            let f = a[r][col].clone();
89            for j in 0..n {
90                let (pa, pi) = (a[col][j].mul(&f), inv[col][j].mul(&f));
91                a[r][j] = a[r][j].sub(&pa);
92                inv[r][j] = inv[r][j].sub(&pi);
93            }
94        }
95    }
96    Some(inv)
97}
98
99/// `e^{t·M}`, every entry exact; `None` for a matrix that is not square.
100pub fn exp_matrix(m: &[Vec<Rational>], t: &Rational) -> Option<Vec<Vec<ExactReal>>> {
101    let n = m.len();
102    let chi = characteristic_polynomial(m)?;
103    let mut nodes: Vec<(ComplexAlgebraic, usize)> = Vec::new();
104    for (factor, multiplicity) in chi.factor_over_q() {
105        for lambda in complex_roots(&factor) {
106            nodes.push((lambda, multiplicity as usize));
107        }
108    }
109    let mut rows: Vec<Vec<ComplexAlgebraic>> = Vec::with_capacity(n);
110    let mut labels: Vec<(usize, usize)> = Vec::with_capacity(n);
111    for (i, (lambda, multiplicity)) in nodes.iter().enumerate() {
112        for j in 0..*multiplicity {
113            rows.push(
114                (0..n)
115                    .map(|k| {
116                        let f = falling(k, j);
117                        if f == 0 {
118                            ComplexAlgebraic::zero()
119                        } else {
120                            lambda.pow((k - j) as u32).scale(&Rational::from_i64(f))
121                        }
122                    })
123                    .collect(),
124            );
125            labels.push((i, j));
126        }
127    }
128    let w = invert(rows)?;
129    let identity: Vec<Vec<Rational>> = (0..n)
130        .map(|i| {
131            (0..n)
132                .map(|j| if i == j { Rational::one() } else { Rational::zero() })
133                .collect()
134        })
135        .collect();
136    let mut powers = vec![identity];
137    for _ in 1..n {
138        let last = powers.last().expect("the identity is there");
139        let next: Vec<Vec<Rational>> = (0..n)
140            .map(|i| {
141                (0..n)
142                    .map(|j| (0..n).fold(Rational::zero(), |acc, k| acc.add(&last[i][k].mul(&m[k][j]))))
143                    .collect()
144            })
145            .collect();
146        powers.push(next);
147    }
148    let mut out = vec![vec![ExactReal::from_i64(0); n]; n];
149    for (r, row) in out.iter_mut().enumerate() {
150        for (s, entry) in row.iter_mut().enumerate() {
151            let mut weights = vec![ComplexAlgebraic::zero(); nodes.len()];
152            for (col, (node, j)) in labels.iter().enumerate() {
153                let alpha = (0..n).fold(ComplexAlgebraic::zero(), |acc, k| {
154                    acc.add(&w[k][col].scale(&powers[k][r][s]))
155                });
156                weights[*node] = weights[*node].add(&alpha.scale(&t.pow(*j as i32).expect("t^j")));
157            }
158            *entry = assemble(&nodes, &weights, t);
159        }
160    }
161    Some(out)
162}
163
164/// `Σ_λ β_λ·e^{t·λ}`, a conjugate pair taken once by its upper root.
165fn assemble(nodes: &[(ComplexAlgebraic, usize)], weights: &[ComplexAlgebraic], t: &Rational) -> ExactReal {
166    let t_alg = Algebraic::from_rational(t.clone());
167    let mut total = ExactReal::from_i64(0);
168    for ((lambda, _), beta) in nodes.iter().zip(weights) {
169        if beta.is_zero() {
170            continue;
171        }
172        let growth = ExactReal::from_algebraic(lambda.re.mul(&t_alg)).exp();
173        if lambda.is_real() {
174            debug_assert!(beta.is_real(), "a real eigenvalue's weight is real");
175            total = total.add(&ExactReal::from_algebraic(beta.re.clone()).mul(&growth));
176        } else if lambda.im.signum() > 0 {
177            let angle = ExactReal::from_algebraic(lambda.im.mul(&t_alg));
178            let two = Algebraic::from_i64(2);
179            let cosine = ExactReal::from_algebraic(beta.re.mul(&two)).mul(&angle.cos());
180            let sine = ExactReal::from_algebraic(beta.im.mul(&two)).mul(&angle.sin());
181            total = total.add(&growth.mul(&cosine.sub(&sine)));
182        }
183    }
184    total
185}
186
187#[cfg(test)]
188#[path = "../../logicaffeine_tests/unit/logicaffeine_base/matrix_exp/tests.rs"]
189mod tests;
logicaffeine_base/matrix_exp.rs

logicaffeine_base/
matrix_exp.rs
