///|
/// The matrix exponential of a square matrix, by scaling and squaring followed
/// by a truncated Taylor series.
///
/// This is what turns a continuous state matrix into a discrete one; the
/// scaling keeps the series inside the range where twenty terms reach machine
/// precision.
pub fn expm(a : Mat, terms~ : Int32 = 20) -> Mat raise Refused {
let n = a.rows
if a.rows != a.cols {
raise Refused::Dimension(got=a.cols, want=a.rows)
}
if n == 0 {
raise Refused::Empty
}
let mut norm = 0.0
for i in 0.. norm {
norm = row_sum
}
}
let mut scale = 1.0
let mut squarings = 0
while norm / scale > 0.5 {
scale *= 2.0
squarings += 1
}
let scaled = a.scale(1.0 / scale)
let mut result = Mat::identity(n)
let mut term = Mat::identity(n)
for k in 1..<=terms {
term = scaled.mul(term).scale(1.0 / k.to_double())
result = result.add(term)
}
let mut out = result
for _ in 0..