From c0a6df55b17b767340e2bd65e98df136b356ff2c Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Wed, 1 Apr 2020 22:36:05 +0200
Subject: [PATCH 01/10] Addition of matrix exponent for static size matrices.

---
 src/linalg/exp.rs   | 445 ++++++++++++++++++++++++++++++++++++++++++++
 src/linalg/mod.rs   |   1 +
 tests/linalg/exp.rs |  97 ++++++++++
 tests/linalg/mod.rs |   2 +
 4 files changed, 545 insertions(+)
 create mode 100644 src/linalg/exp.rs
 create mode 100644 tests/linalg/exp.rs

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
new file mode 100644
index 00000000..8e8227d5
--- /dev/null
+++ b/src/linalg/exp.rs
@@ -0,0 +1,445 @@
+use crate::{
+    base::{
+        allocator::Allocator,
+        dimension::{DimMin, DimMinimum, DimName},
+        DefaultAllocator,
+    },
+    try_convert, ComplexField, MatrixN, RealField,
+};
+
+// https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/matfuncs.py
+struct ExpmPadeHelper<N, R>
+where
+    N: RealField,
+    R: DimName + DimMin<R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    use_exact_norm: bool,
+    ident: MatrixN<N, R>,
+
+    a: MatrixN<N, R>,
+    a2: Option<MatrixN<N, R>>,
+    a4: Option<MatrixN<N, R>>,
+    a6: Option<MatrixN<N, R>>,
+    a8: Option<MatrixN<N, R>>,
+    a10: Option<MatrixN<N, R>>,
+
+    d4_exact: Option<N>,
+    d6_exact: Option<N>,
+    d8_exact: Option<N>,
+    d10_exact: Option<N>,
+
+    d4_approx: Option<N>,
+    d6_approx: Option<N>,
+    d8_approx: Option<N>,
+    d10_approx: Option<N>,
+}
+
+impl<N, R> ExpmPadeHelper<N, R>
+where
+    N: RealField,
+    R: DimName + DimMin<R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    fn new(a: MatrixN<N, R>, use_exact_norm: bool) -> Self {
+        ExpmPadeHelper {
+            use_exact_norm,
+            ident: MatrixN::<N, R>::identity(),
+            a,
+            a2: None,
+            a4: None,
+            a6: None,
+            a8: None,
+            a10: None,
+            d4_exact: None,
+            d6_exact: None,
+            d8_exact: None,
+            d10_exact: None,
+            d4_approx: None,
+            d6_approx: None,
+            d8_approx: None,
+            d10_approx: None,
+        }
+    }
+
+    fn a2(&self) -> &MatrixN<N, R> {
+        if self.a2.is_none() {
+            let ap = &self.a2 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(&self.a * &self.a);
+            }
+        }
+        self.a2.as_ref().unwrap()
+    }
+
+    fn a4(&self) -> &MatrixN<N, R> {
+        if self.a4.is_none() {
+            let ap = &self.a4 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let a2 = self.a2();
+            unsafe {
+                *ap = Some(a2 * a2);
+            }
+        }
+        self.a4.as_ref().unwrap()
+    }
+
+    fn a6(&self) -> &MatrixN<N, R> {
+        if self.a6.is_none() {
+            let a2 = self.a2();
+            let a4 = self.a4();
+            let ap = &self.a6 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(a4 * a2);
+            }
+        }
+        self.a6.as_ref().unwrap()
+    }
+
+    fn a8(&self) -> &MatrixN<N, R> {
+        if self.a8.is_none() {
+            let a2 = self.a2();
+            let a6 = self.a6();
+            let ap = &self.a8 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(a6 * a2);
+            }
+        }
+        self.a8.as_ref().unwrap()
+    }
+
+    fn a10(&mut self) -> &MatrixN<N, R> {
+        if self.a10.is_none() {
+            let a4 = self.a4();
+            let a6 = self.a6();
+            let ap = &self.a10 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            unsafe {
+                *ap = Some(a6 * a4);
+            }
+        }
+        self.a10.as_ref().unwrap()
+    }
+
+    fn d4_tight(&mut self) -> N {
+        if self.d4_exact.is_none() {
+            self.d4_exact = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+        }
+        self.d4_exact.unwrap()
+    }
+
+    fn d6_tight(&mut self) -> N {
+        if self.d6_exact.is_none() {
+            self.d6_exact = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+        }
+        self.d6_exact.unwrap()
+    }
+
+    fn d8_tight(&mut self) -> N {
+        if self.d8_exact.is_none() {
+            self.d8_exact = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+        }
+        self.d8_exact.unwrap()
+    }
+
+    fn d10_tight(&mut self) -> N {
+        if self.d10_exact.is_none() {
+            self.d10_exact = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+        }
+        self.d10_exact.unwrap()
+    }
+
+    fn d4_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d4_tight();
+        }
+
+        if self.d4_exact.is_some() {
+            return self.d4_exact.unwrap();
+        }
+
+        if self.d4_approx.is_none() {
+            self.d4_approx = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+        }
+
+        self.d4_approx.unwrap()
+    }
+
+    fn d6_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d6_tight();
+        }
+
+        if self.d6_exact.is_some() {
+            return self.d6_exact.unwrap();
+        }
+
+        if self.d6_approx.is_none() {
+            self.d6_approx = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+        }
+
+        self.d6_approx.unwrap()
+    }
+
+    fn d8_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d8_tight();
+        }
+
+        if self.d8_exact.is_some() {
+            return self.d8_exact.unwrap();
+        }
+
+        if self.d8_approx.is_none() {
+            self.d8_approx = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+        }
+
+        self.d8_approx.unwrap()
+    }
+
+    fn d10_loose(&mut self) -> N {
+        if self.use_exact_norm {
+            return self.d10_tight();
+        }
+
+        if self.d10_exact.is_some() {
+            return self.d10_exact.unwrap();
+        }
+
+        if self.d10_approx.is_none() {
+            self.d10_approx = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+        }
+
+        self.d10_approx.unwrap()
+    }
+
+    fn pade3(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(120.0).unwrap(),
+            N::from_f64(60.0).unwrap(),
+            N::from_f64(12.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u = &self.a * (self.a2() * b[3] + &self.ident * b[1]);
+        let v = self.a2() * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade5(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(30240.0).unwrap(),
+            N::from_f64(15120.0).unwrap(),
+            N::from_f64(3360.0).unwrap(),
+            N::from_f64(420.0).unwrap(),
+            N::from_f64(30.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u = &self.a * (self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
+        let v = self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade7(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(17297280.0).unwrap(),
+            N::from_f64(8648640.0).unwrap(),
+            N::from_f64(1995840.0).unwrap(),
+            N::from_f64(277200.0).unwrap(),
+            N::from_f64(25200.0).unwrap(),
+            N::from_f64(1512.0).unwrap(),
+            N::from_f64(56.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u =
+            &self.a * (self.a6() * b[7] + self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
+        let v = self.a6() * b[6] + self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade9(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(17643225600.0).unwrap(),
+            N::from_f64(8821612800.0).unwrap(),
+            N::from_f64(2075673600.0).unwrap(),
+            N::from_f64(302702400.0).unwrap(),
+            N::from_f64(30270240.0).unwrap(),
+            N::from_f64(2162160.0).unwrap(),
+            N::from_f64(110880.0).unwrap(),
+            N::from_f64(3960.0).unwrap(),
+            N::from_f64(90.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let u = &self.a
+            * (self.a8() * b[9]
+                + self.a6() * b[7]
+                + self.a4() * b[5]
+                + self.a2() * b[3]
+                + &self.ident * b[1]);
+        let v = self.a8() * b[8]
+            + self.a6() * b[6]
+            + self.a4() * b[4]
+            + self.a2() * b[2]
+            + &self.ident * b[0];
+        (u, v)
+    }
+
+    fn pade13_scaled(&mut self, s: u64) -> (MatrixN<N, R>, MatrixN<N, R>) {
+        let b = [
+            N::from_f64(64764752532480000.0).unwrap(),
+            N::from_f64(32382376266240000.0).unwrap(),
+            N::from_f64(7771770303897600.0).unwrap(),
+            N::from_f64(1187353796428800.0).unwrap(),
+            N::from_f64(129060195264000.0).unwrap(),
+            N::from_f64(10559470521600.0).unwrap(),
+            N::from_f64(670442572800.0).unwrap(),
+            N::from_f64(33522128640.0).unwrap(),
+            N::from_f64(1323241920.0).unwrap(),
+            N::from_f64(40840800.0).unwrap(),
+            N::from_f64(960960.0).unwrap(),
+            N::from_f64(16380.0).unwrap(),
+            N::from_f64(182.0).unwrap(),
+            N::from_f64(1.0).unwrap(),
+        ];
+        let s = s as f64;
+
+        let mb = &self.a * N::from_f64(2.0.powf(-s)).unwrap();
+        let mb2 = self.a2() * N::from_f64(2.0.powf(-2.0 * s)).unwrap();
+        let mb4 = self.a4() * N::from_f64(2.0.powf(-4.0 * s)).unwrap();
+        let mb6 = self.a6() * N::from_f64(2.0.powf(-6.0 * s)).unwrap();
+
+        let u2 = &mb6 * (&mb6 * b[13] + &mb4 * b[11] + &mb2 * b[9]);
+        let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]);
+        let v2 = &mb6 * (&mb6 * b[12] + &mb4 * b[10] + &mb2 * b[8]);
+        let v = v2 + &mb6 * b[6] + &mb4 * b[4] + &mb2 * b[2] + &self.ident * b[0];
+        (u, v)
+    }
+}
+
+fn factorial(n: u128) -> u128 {
+    if n == 1 {
+        return 1;
+    }
+    n * factorial(n - 1)
+}
+
+fn onenorm_matrix_power_nnm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
+where
+    N: RealField,
+    R: DimName,
+    DefaultAllocator: Allocator<N, R, R>,
+{
+    let mut v = MatrixN::<N, R>::repeat(N::from_f64(1.0).unwrap());
+    let m = a.transpose();
+
+    for _ in 0..p {
+        v = &m * v;
+    }
+
+    v.amax()
+}
+
+fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
+where
+    N: RealField,
+    R: DimName,
+    DefaultAllocator: Allocator<N, R, R>,
+{
+    // 2m choose m = (2m)!/(m! * (2m-m)!)
+
+    let a_abs_onenorm = onenorm_matrix_power_nnm(&a.abs(), 2 * m + 1);
+
+    if a_abs_onenorm == N::zero() {
+        return 0;
+    }
+
+    let choose_2m_m =
+        factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128));
+    let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1);
+    let alpha = a_abs_onenorm / a.amax();
+    let alpha = alpha / N::from_u128(abs_c_recip).unwrap();
+
+    let u = N::from_f64(2_f64.powf(-53.0)).unwrap();
+    let log2_alpha_div_u = try_convert((alpha / u).log2()).unwrap();
+    let value = (log2_alpha_div_u / (2.0 * m as f64)).ceil();
+    if value > 0.0 {
+        value as u64
+    } else {
+        0
+    }
+}
+
+fn solve_p_q<N, R>(u: MatrixN<N, R>, v: MatrixN<N, R>) -> MatrixN<N, R>
+where
+    N: ComplexField,
+    R: DimMin<R, Output = R> + DimName,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    let p = &u + &v;
+    let q = &v - &u;
+
+    q.lu().solve(&p).unwrap()
+}
+
+impl<N: RealField, R: DimMin<R, Output = R> + DimName> MatrixN<N, R>
+where
+    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+{
+    /// Computes exp of this matrix
+    pub fn exp(&self) -> Self {
+        // Simple case
+        if self.nrows() == 1 {
+            return self.clone().map(|v| v.exp());
+        }
+
+        let mut h = ExpmPadeHelper::new(self.clone(), true);
+
+        let eta_1 = N::max(h.d4_loose(), h.d6_loose());
+        if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
+            let (u, v) = h.pade3();
+            return solve_p_q(u, v);
+        }
+
+        let eta_2 = N::max(h.d4_tight(), h.d6_loose());
+        if eta_2 < N::from_f64(2.539398330063230e-001).unwrap() && ell(&h.a, 5) == 0 {
+            let (u, v) = h.pade5();
+            return solve_p_q(u, v);
+        }
+
+        let eta_3 = N::max(h.d6_tight(), h.d8_loose());
+        if eta_3 < N::from_f64(9.504178996162932e-001).unwrap() && ell(&h.a, 7) == 0 {
+            let (u, v) = h.pade7();
+            return solve_p_q(u, v);
+        }
+        if eta_3 < N::from_f64(2.097847961257068e+000).unwrap() && ell(&h.a, 9) == 0 {
+            let (u, v) = h.pade9();
+            return solve_p_q(u, v);
+        }
+
+        let eta_4 = N::max(h.d8_loose(), h.d10_loose());
+        let eta_5 = N::min(eta_3, eta_4);
+        let theta_13 = N::from_f64(4.25).unwrap();
+
+        let mut s = if eta_5 == N::zero() {
+            0
+        } else {
+            let l2 = try_convert((eta_5 / theta_13).log2().ceil()).unwrap();
+
+            if l2 < 0.0 {
+                0
+            } else {
+                l2 as u64
+            }
+        };
+
+        s += ell(
+            &(&h.a * N::from_f64(2.0_f64.powf(-(s as f64))).unwrap()),
+            13,
+        );
+
+        let (u, v) = h.pade13_scaled(s);
+        let mut x = solve_p_q(u, v);
+
+        for _ in 0..s {
+            x = &x * &x;
+        }
+        x
+    }
+}
diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index de1108f7..42bfde63 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -5,6 +5,7 @@ mod bidiagonal;
 mod cholesky;
 mod convolution;
 mod determinant;
+mod exp;
 mod full_piv_lu;
 pub mod givens;
 mod hessenberg;
diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs
new file mode 100644
index 00000000..ac3f8dd7
--- /dev/null
+++ b/tests/linalg/exp.rs
@@ -0,0 +1,97 @@
+#[cfg(test)]
+mod tests {
+    //https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py
+    #[test]
+    fn exp() {
+        use nalgebra::{Matrix1, Matrix2, Matrix3};
+
+        {
+            let m = Matrix1::new(1.0);
+
+            let f = m.exp();
+
+            assert_eq!(f, Matrix1::new(1_f64.exp()));
+        }
+
+        {
+            let m = Matrix2::new(0.0, 1.0, 0.0, 0.0);
+
+            assert_eq!(m.exp(), Matrix2::new(1.0, 1.0, 0.0, 1.0));
+        }
+
+        {
+            let a: f64 = 1.0;
+            let b: f64 = 2.0;
+            let c: f64 = 3.0;
+            let d: f64 = 4.0;
+            let m = Matrix2::new(a, b, c, d);
+
+            let delta = ((a - d).powf(2.0) + 4.0 * b * c).sqrt();
+            let delta_2 = delta / 2.0;
+            let ad_2 = (a + d) / 2.0;
+            let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh());
+            let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh();
+            let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh();
+            let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
+
+            let f = Matrix2::new(m11, m12, m21, m22) / delta;
+            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+        }
+
+        {
+            // https://mathworld.wolfram.com/MatrixExponential.html
+            use rand::{
+                distributions::{Distribution, Uniform},
+                thread_rng,
+            };
+            let mut rng = thread_rng();
+            let dist = Uniform::new(-10.0, 10.0);
+            loop {
+                let a: f64 = dist.sample(&mut rng);
+                let b: f64 = dist.sample(&mut rng);
+                let c: f64 = dist.sample(&mut rng);
+                let d: f64 = dist.sample(&mut rng);
+                let m = Matrix2::new(a, b, c, d);
+
+                let delta_sq = (a - d).powf(2.0) + 4.0 * b * c;
+                if delta_sq < 0.0 {
+                    continue;
+                }
+
+                let delta = delta_sq.sqrt();
+                let delta_2 = delta / 2.0;
+                let ad_2 = (a + d) / 2.0;
+                let m11 = ad_2.exp() * (delta * delta_2.cosh() + (a - d) * delta_2.sinh());
+                let m12 = 2.0 * b * ad_2.exp() * delta_2.sinh();
+                let m21 = 2.0 * c * ad_2.exp() * delta_2.sinh();
+                let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
+
+                let f = Matrix2::new(m11, m12, m21, m22) / delta;
+                println!("a: {}", m);
+                assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+                break;
+            }
+        }
+
+        {
+            let m = Matrix3::new(1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0);
+
+            let e1 = 1.0_f64.exp();
+            let e2 = 2.0_f64.exp();
+
+            let f = Matrix3::new(
+                e1,
+                3.0 * e1,
+                15.0 * (e2 - 2.0 * e1),
+                0.0,
+                e1,
+                5.0 * (e2 - e1),
+                0.0,
+                0.0,
+                e2,
+            );
+
+            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+        }
+    }
+}
diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs
index 234cac39..c810df1f 100644
--- a/tests/linalg/mod.rs
+++ b/tests/linalg/mod.rs
@@ -1,7 +1,9 @@
 mod balancing;
 mod bidiagonal;
 mod cholesky;
+mod convolution;
 mod eigen;
+mod exp;
 mod full_piv_lu;
 mod hessenberg;
 mod inverse;

From 2696ffd7a1837b73756b9512b32781532b4c66ca Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Thu, 2 Apr 2020 09:38:18 +0200
Subject: [PATCH 02/10] Fixed bug by introducing one norm

---
 src/linalg/exp.rs | 51 +++++++++++++++++++++++++++++++++++++----------
 1 file changed, 40 insertions(+), 11 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 8e8227d5..877ce0d3 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -121,28 +121,28 @@ where
 
     fn d4_tight(&mut self) -> N {
         if self.d4_exact.is_none() {
-            self.d4_exact = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+            self.d4_exact = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
         }
         self.d4_exact.unwrap()
     }
 
     fn d6_tight(&mut self) -> N {
         if self.d6_exact.is_none() {
-            self.d6_exact = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_exact = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
         }
         self.d6_exact.unwrap()
     }
 
     fn d8_tight(&mut self) -> N {
         if self.d8_exact.is_none() {
-            self.d8_exact = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_exact = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
         }
         self.d8_exact.unwrap()
     }
 
     fn d10_tight(&mut self) -> N {
         if self.d10_exact.is_none() {
-            self.d10_exact = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_exact = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
         }
         self.d10_exact.unwrap()
     }
@@ -157,7 +157,7 @@ where
         }
 
         if self.d4_approx.is_none() {
-            self.d4_approx = Some(self.a4().amax().powf(N::from_f64(0.25).unwrap()));
+            self.d4_approx = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
         }
 
         self.d4_approx.unwrap()
@@ -173,7 +173,7 @@ where
         }
 
         if self.d6_approx.is_none() {
-            self.d6_approx = Some(self.a6().amax().powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_approx = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
         }
 
         self.d6_approx.unwrap()
@@ -189,7 +189,7 @@ where
         }
 
         if self.d8_approx.is_none() {
-            self.d8_approx = Some(self.a8().amax().powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_approx = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
         }
 
         self.d8_approx.unwrap()
@@ -205,7 +205,7 @@ where
         }
 
         if self.d10_approx.is_none() {
-            self.d10_approx = Some(self.a10().amax().powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_approx = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
         }
 
         self.d10_approx.unwrap()
@@ -333,7 +333,7 @@ where
         v = &m * v;
     }
 
-    v.amax()
+    one_norm(&v)
 }
 
 fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
@@ -353,7 +353,7 @@ where
     let choose_2m_m =
         factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128));
     let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1);
-    let alpha = a_abs_onenorm / a.amax();
+    let alpha = a_abs_onenorm / one_norm(a);
     let alpha = alpha / N::from_u128(abs_c_recip).unwrap();
 
     let u = N::from_f64(2_f64.powf(-53.0)).unwrap();
@@ -378,6 +378,24 @@ where
     q.lu().solve(&p).unwrap()
 }
 
+pub fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
+where
+    N: RealField,
+    R: DimName,
+    DefaultAllocator: Allocator<N, R, R>,
+{
+    let mut col_sums = vec![N::zero(); m.ncols()];
+    for i in 0..m.ncols() {
+        let col = m.column(i);
+        col.iter().for_each(|v| col_sums[i] += v.abs());
+    }
+    let mut max = col_sums[0];
+    for i in 1..col_sums.len() {
+        max = N::max(max, col_sums[i]);
+    }
+    max
+}
+
 impl<N: RealField, R: DimMin<R, Output = R> + DimName> MatrixN<N, R>
 where
     DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
@@ -389,7 +407,7 @@ where
             return self.clone().map(|v| v.exp());
         }
 
-        let mut h = ExpmPadeHelper::new(self.clone(), true);
+        let mut h = ExpmPadeHelper::new(self.clone(), false);
 
         let eta_1 = N::max(h.d4_loose(), h.d6_loose());
         if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
@@ -443,3 +461,14 @@ where
         x
     }
 }
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn one_norm() {
+        use crate::Matrix3;
+        let m = Matrix3::new(-3.0, 5.0, 7.0, 2.0, 6.0, 4.0, 0.0, 2.0, 8.0);
+
+        assert_eq!(super::one_norm(&m), 19.0);
+    }
+}

From e1565616778f4635088d47f88bb46c1ff8bc4c4d Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Thu, 2 Apr 2020 09:39:50 +0200
Subject: [PATCH 03/10] one_norm is not a public function

---
 src/linalg/exp.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 877ce0d3..5b76dcaf 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -378,7 +378,7 @@ where
     q.lu().solve(&p).unwrap()
 }
 
-pub fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
+fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
 where
     N: RealField,
     R: DimName,

From b3ef66fd14b17ffb54f8cb560a5866c81a6aff4a Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Thu, 2 Apr 2020 10:55:00 +0200
Subject: [PATCH 04/10] Fixed bug in onenorm_matrix_power_norm

---
 src/linalg/exp.rs | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 5b76dcaf..7c087159 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -320,31 +320,32 @@ fn factorial(n: u128) -> u128 {
     n * factorial(n - 1)
 }
 
-fn onenorm_matrix_power_nnm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
+/// Compute the 1-norm of a non-negative integer power of a non-negative matrix.
+fn onenorm_matrix_power_nonm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
 where
     N: RealField,
     R: DimName,
-    DefaultAllocator: Allocator<N, R, R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
 {
-    let mut v = MatrixN::<N, R>::repeat(N::from_f64(1.0).unwrap());
+    let mut v = crate::VectorN::<N, R>::repeat(N::from_f64(1.0).unwrap());
     let m = a.transpose();
 
     for _ in 0..p {
         v = &m * v;
     }
 
-    one_norm(&v)
+    v.max()
 }
 
 fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
 where
     N: RealField,
     R: DimName,
-    DefaultAllocator: Allocator<N, R, R>,
+    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
 {
     // 2m choose m = (2m)!/(m! * (2m-m)!)
 
-    let a_abs_onenorm = onenorm_matrix_power_nnm(&a.abs(), 2 * m + 1);
+    let a_abs_onenorm = onenorm_matrix_power_nonm(&a.abs(), 2 * m + 1);
 
     if a_abs_onenorm == N::zero() {
         return 0;
@@ -396,9 +397,11 @@ where
     max
 }
 
-impl<N: RealField, R: DimMin<R, Output = R> + DimName> MatrixN<N, R>
+impl<N: RealField, R> MatrixN<N, R>
 where
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    R: DimMin<R, Output = R> + DimName,
+    DefaultAllocator:
+        Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>> + Allocator<N, R>,
 {
     /// Computes exp of this matrix
     pub fn exp(&self) -> Self {

From 4ec94408b5652e12b779df462479e15df8bc7f33 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Fri, 3 Apr 2020 11:11:14 +0200
Subject: [PATCH 05/10] Pub use of exp in the linalg module

---
 src/linalg/exp.rs | 6 ++++--
 src/linalg/mod.rs | 1 +
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 7c087159..8822027b 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -1,3 +1,5 @@
+//! This module provides the matrix exponent (exp) function to square matrices.
+//!
 use crate::{
     base::{
         allocator::Allocator,
@@ -403,14 +405,14 @@ where
     DefaultAllocator:
         Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>> + Allocator<N, R>,
 {
-    /// Computes exp of this matrix
+    /// Computes exponential of this matrix
     pub fn exp(&self) -> Self {
         // Simple case
         if self.nrows() == 1 {
             return self.clone().map(|v| v.exp());
         }
 
-        let mut h = ExpmPadeHelper::new(self.clone(), false);
+        let mut h = ExpmPadeHelper::new(self.clone(), true);
 
         let eta_1 = N::max(h.d4_loose(), h.d6_loose());
         if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index 42bfde63..f96cef0c 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -27,6 +27,7 @@ mod symmetric_tridiagonal;
 pub use self::bidiagonal::*;
 pub use self::cholesky::*;
 pub use self::convolution::*;
+pub use self::exp::*;
 pub use self::full_piv_lu::*;
 pub use self::hessenberg::*;
 pub use self::lu::*;

From dbbf87a3dd475c454a5f69cfb355f9fd6b4da5c3 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Tue, 7 Apr 2020 09:44:06 +0200
Subject: [PATCH 06/10] Rebased against dev

---
 tests/linalg/mod.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs
index c810df1f..7fc01396 100644
--- a/tests/linalg/mod.rs
+++ b/tests/linalg/mod.rs
@@ -12,5 +12,4 @@ mod qr;
 mod schur;
 mod solve;
 mod svd;
-mod convolution;
 mod tridiagonal;

From c7d9e415ce3531b2ffb8399c2c4e3682c54068bf Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Tue, 7 Apr 2020 09:55:58 +0200
Subject: [PATCH 07/10] Converted tests to use relative_eq

---
 tests/linalg/exp.rs | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs
index ac3f8dd7..37abc4d2 100644
--- a/tests/linalg/exp.rs
+++ b/tests/linalg/exp.rs
@@ -10,13 +10,17 @@ mod tests {
 
             let f = m.exp();
 
-            assert_eq!(f, Matrix1::new(1_f64.exp()));
+            assert!(relative_eq!(f, Matrix1::new(1_f64.exp()), epsilon = 1.0e-7));
         }
 
         {
             let m = Matrix2::new(0.0, 1.0, 0.0, 0.0);
 
-            assert_eq!(m.exp(), Matrix2::new(1.0, 1.0, 0.0, 1.0));
+            assert!(relative_eq!(
+                m.exp(),
+                Matrix2::new(1.0, 1.0, 0.0, 1.0),
+                epsilon = 1.0e-7
+            ));
         }
 
         {
@@ -35,7 +39,7 @@ mod tests {
             let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
 
             let f = Matrix2::new(m11, m12, m21, m22) / delta;
-            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+            assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
         }
 
         {
@@ -68,7 +72,7 @@ mod tests {
 
                 let f = Matrix2::new(m11, m12, m21, m22) / delta;
                 println!("a: {}", m);
-                assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+                assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
                 break;
             }
         }
@@ -91,7 +95,7 @@ mod tests {
                 e2,
             );
 
-            assert!((f - m.exp()).iter().all(|v| v.abs() <= 0.00005));
+            assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
         }
     }
 }

From 0a3ee99cdb0dfa9b9202ced71ae5c277dc18b0f4 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Sun, 12 Apr 2020 11:46:00 +0200
Subject: [PATCH 08/10] Changed dimension name R to D Changed N::from_x to
 crate::convert

---
 src/linalg/exp.rs | 230 ++++++++++++++++++++++------------------------
 1 file changed, 111 insertions(+), 119 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 8822027b..f68cbe6c 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -6,25 +6,25 @@ use crate::{
         dimension::{DimMin, DimMinimum, DimName},
         DefaultAllocator,
     },
-    try_convert, ComplexField, MatrixN, RealField,
+    convert, try_convert, ComplexField, MatrixN, RealField,
 };
 
 // https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/matfuncs.py
-struct ExpmPadeHelper<N, R>
+struct ExpmPadeHelper<N, D>
 where
     N: RealField,
-    R: DimName + DimMin<R>,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    D: DimName + DimMin<D>,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     use_exact_norm: bool,
-    ident: MatrixN<N, R>,
+    ident: MatrixN<N, D>,
 
-    a: MatrixN<N, R>,
-    a2: Option<MatrixN<N, R>>,
-    a4: Option<MatrixN<N, R>>,
-    a6: Option<MatrixN<N, R>>,
-    a8: Option<MatrixN<N, R>>,
-    a10: Option<MatrixN<N, R>>,
+    a: MatrixN<N, D>,
+    a2: Option<MatrixN<N, D>>,
+    a4: Option<MatrixN<N, D>>,
+    a6: Option<MatrixN<N, D>>,
+    a8: Option<MatrixN<N, D>>,
+    a10: Option<MatrixN<N, D>>,
 
     d4_exact: Option<N>,
     d6_exact: Option<N>,
@@ -37,16 +37,16 @@ where
     d10_approx: Option<N>,
 }
 
-impl<N, R> ExpmPadeHelper<N, R>
+impl<N, D> ExpmPadeHelper<N, D>
 where
     N: RealField,
-    R: DimName + DimMin<R>,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    D: DimName + DimMin<D>,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
-    fn new(a: MatrixN<N, R>, use_exact_norm: bool) -> Self {
+    fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
         ExpmPadeHelper {
             use_exact_norm,
-            ident: MatrixN::<N, R>::identity(),
+            ident: MatrixN::<N, D>::identity(),
             a,
             a2: None,
             a4: None,
@@ -64,9 +64,9 @@ where
         }
     }
 
-    fn a2(&self) -> &MatrixN<N, R> {
+    fn a2(&self) -> &MatrixN<N, D> {
         if self.a2.is_none() {
-            let ap = &self.a2 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a2 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(&self.a * &self.a);
             }
@@ -74,9 +74,9 @@ where
         self.a2.as_ref().unwrap()
     }
 
-    fn a4(&self) -> &MatrixN<N, R> {
+    fn a4(&self) -> &MatrixN<N, D> {
         if self.a4.is_none() {
-            let ap = &self.a4 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a4 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             let a2 = self.a2();
             unsafe {
                 *ap = Some(a2 * a2);
@@ -85,11 +85,11 @@ where
         self.a4.as_ref().unwrap()
     }
 
-    fn a6(&self) -> &MatrixN<N, R> {
+    fn a6(&self) -> &MatrixN<N, D> {
         if self.a6.is_none() {
             let a2 = self.a2();
             let a4 = self.a4();
-            let ap = &self.a6 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a6 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(a4 * a2);
             }
@@ -97,11 +97,11 @@ where
         self.a6.as_ref().unwrap()
     }
 
-    fn a8(&self) -> &MatrixN<N, R> {
+    fn a8(&self) -> &MatrixN<N, D> {
         if self.a8.is_none() {
             let a2 = self.a2();
             let a6 = self.a6();
-            let ap = &self.a8 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a8 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(a6 * a2);
             }
@@ -109,11 +109,11 @@ where
         self.a8.as_ref().unwrap()
     }
 
-    fn a10(&mut self) -> &MatrixN<N, R> {
+    fn a10(&mut self) -> &MatrixN<N, D> {
         if self.a10.is_none() {
             let a4 = self.a4();
             let a6 = self.a6();
-            let ap = &self.a10 as *const Option<MatrixN<N, R>> as *mut Option<MatrixN<N, R>>;
+            let ap = &self.a10 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
             unsafe {
                 *ap = Some(a6 * a4);
             }
@@ -123,28 +123,28 @@ where
 
     fn d4_tight(&mut self) -> N {
         if self.d4_exact.is_none() {
-            self.d4_exact = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
+            self.d4_exact = Some(one_norm(self.a4()).powf(convert(0.25)));
         }
         self.d4_exact.unwrap()
     }
 
     fn d6_tight(&mut self) -> N {
         if self.d6_exact.is_none() {
-            self.d6_exact = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_exact = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
         }
         self.d6_exact.unwrap()
     }
 
     fn d8_tight(&mut self) -> N {
         if self.d8_exact.is_none() {
-            self.d8_exact = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_exact = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
         }
         self.d8_exact.unwrap()
     }
 
     fn d10_tight(&mut self) -> N {
         if self.d10_exact.is_none() {
-            self.d10_exact = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_exact = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
         }
         self.d10_exact.unwrap()
     }
@@ -159,7 +159,7 @@ where
         }
 
         if self.d4_approx.is_none() {
-            self.d4_approx = Some(one_norm(self.a4()).powf(N::from_f64(0.25).unwrap()));
+            self.d4_approx = Some(one_norm(self.a4()).powf(convert(0.25)));
         }
 
         self.d4_approx.unwrap()
@@ -175,7 +175,7 @@ where
         }
 
         if self.d6_approx.is_none() {
-            self.d6_approx = Some(one_norm(self.a6()).powf(N::from_f64(1.0 / 6.0).unwrap()));
+            self.d6_approx = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
         }
 
         self.d6_approx.unwrap()
@@ -191,7 +191,7 @@ where
         }
 
         if self.d8_approx.is_none() {
-            self.d8_approx = Some(one_norm(self.a8()).powf(N::from_f64(1.0 / 8.0).unwrap()));
+            self.d8_approx = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
         }
 
         self.d8_approx.unwrap()
@@ -207,48 +207,43 @@ where
         }
 
         if self.d10_approx.is_none() {
-            self.d10_approx = Some(one_norm(self.a10()).powf(N::from_f64(1.0 / 10.0).unwrap()));
+            self.d10_approx = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
         }
 
         self.d10_approx.unwrap()
     }
 
-    fn pade3(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(120.0).unwrap(),
-            N::from_f64(60.0).unwrap(),
-            N::from_f64(12.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
-        ];
+    fn pade3(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 4] = [convert(120.0), convert(60.0), convert(12.0), convert(1.0)];
         let u = &self.a * (self.a2() * b[3] + &self.ident * b[1]);
         let v = self.a2() * b[2] + &self.ident * b[0];
         (u, v)
     }
 
-    fn pade5(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(30240.0).unwrap(),
-            N::from_f64(15120.0).unwrap(),
-            N::from_f64(3360.0).unwrap(),
-            N::from_f64(420.0).unwrap(),
-            N::from_f64(30.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade5(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 6] = [
+            convert(30240.0),
+            convert(15120.0),
+            convert(3360.0),
+            convert(420.0),
+            convert(30.0),
+            convert(1.0),
         ];
         let u = &self.a * (self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
         let v = self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
         (u, v)
     }
 
-    fn pade7(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(17297280.0).unwrap(),
-            N::from_f64(8648640.0).unwrap(),
-            N::from_f64(1995840.0).unwrap(),
-            N::from_f64(277200.0).unwrap(),
-            N::from_f64(25200.0).unwrap(),
-            N::from_f64(1512.0).unwrap(),
-            N::from_f64(56.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade7(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 8] = [
+            convert(17297280.0),
+            convert(8648640.0),
+            convert(1995840.0),
+            convert(277200.0),
+            convert(25200.0),
+            convert(1512.0),
+            convert(56.0),
+            convert(1.0),
         ];
         let u =
             &self.a * (self.a6() * b[7] + self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
@@ -256,18 +251,18 @@ where
         (u, v)
     }
 
-    fn pade9(&mut self) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(17643225600.0).unwrap(),
-            N::from_f64(8821612800.0).unwrap(),
-            N::from_f64(2075673600.0).unwrap(),
-            N::from_f64(302702400.0).unwrap(),
-            N::from_f64(30270240.0).unwrap(),
-            N::from_f64(2162160.0).unwrap(),
-            N::from_f64(110880.0).unwrap(),
-            N::from_f64(3960.0).unwrap(),
-            N::from_f64(90.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade9(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 10] = [
+            convert(17643225600.0),
+            convert(8821612800.0),
+            convert(2075673600.0),
+            convert(302702400.0),
+            convert(30270240.0),
+            convert(2162160.0),
+            convert(110880.0),
+            convert(3960.0),
+            convert(90.0),
+            convert(1.0),
         ];
         let u = &self.a
             * (self.a8() * b[9]
@@ -283,29 +278,29 @@ where
         (u, v)
     }
 
-    fn pade13_scaled(&mut self, s: u64) -> (MatrixN<N, R>, MatrixN<N, R>) {
-        let b = [
-            N::from_f64(64764752532480000.0).unwrap(),
-            N::from_f64(32382376266240000.0).unwrap(),
-            N::from_f64(7771770303897600.0).unwrap(),
-            N::from_f64(1187353796428800.0).unwrap(),
-            N::from_f64(129060195264000.0).unwrap(),
-            N::from_f64(10559470521600.0).unwrap(),
-            N::from_f64(670442572800.0).unwrap(),
-            N::from_f64(33522128640.0).unwrap(),
-            N::from_f64(1323241920.0).unwrap(),
-            N::from_f64(40840800.0).unwrap(),
-            N::from_f64(960960.0).unwrap(),
-            N::from_f64(16380.0).unwrap(),
-            N::from_f64(182.0).unwrap(),
-            N::from_f64(1.0).unwrap(),
+    fn pade13_scaled(&mut self, s: u64) -> (MatrixN<N, D>, MatrixN<N, D>) {
+        let b: [N; 14] = [
+            convert(64764752532480000.0),
+            convert(32382376266240000.0),
+            convert(7771770303897600.0),
+            convert(1187353796428800.0),
+            convert(129060195264000.0),
+            convert(10559470521600.0),
+            convert(670442572800.0),
+            convert(33522128640.0),
+            convert(1323241920.0),
+            convert(40840800.0),
+            convert(960960.0),
+            convert(16380.0),
+            convert(182.0),
+            convert(1.0),
         ];
         let s = s as f64;
 
-        let mb = &self.a * N::from_f64(2.0.powf(-s)).unwrap();
-        let mb2 = self.a2() * N::from_f64(2.0.powf(-2.0 * s)).unwrap();
-        let mb4 = self.a4() * N::from_f64(2.0.powf(-4.0 * s)).unwrap();
-        let mb6 = self.a6() * N::from_f64(2.0.powf(-6.0 * s)).unwrap();
+        let mb = &self.a * convert::<f64, N>(2.0_f64.powf(-s));
+        let mb2 = self.a2() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
+        let mb4 = self.a4() * convert::<f64, N>(2.0.powf(-4.0 * s));
+        let mb6 = self.a6() * convert::<f64, N>(2.0.powf(-6.0 * s));
 
         let u2 = &mb6 * (&mb6 * b[13] + &mb4 * b[11] + &mb2 * b[9]);
         let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]);
@@ -323,13 +318,13 @@ fn factorial(n: u128) -> u128 {
 }
 
 /// Compute the 1-norm of a non-negative integer power of a non-negative matrix.
-fn onenorm_matrix_power_nonm<N, R>(a: &MatrixN<N, R>, p: u64) -> N
+fn onenorm_matrix_power_nonm<N, D>(a: &MatrixN<N, D>, p: u64) -> N
 where
     N: RealField,
-    R: DimName,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
+    D: DimName,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
-    let mut v = crate::VectorN::<N, R>::repeat(N::from_f64(1.0).unwrap());
+    let mut v = crate::VectorN::<N, D>::repeat(convert(1.0));
     let m = a.transpose();
 
     for _ in 0..p {
@@ -339,11 +334,11 @@ where
     v.max()
 }
 
-fn ell<N, R>(a: &MatrixN<N, R>, m: u64) -> u64
+fn ell<N, D>(a: &MatrixN<N, D>, m: u64) -> u64
 where
     N: RealField,
-    R: DimName,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<N, R>,
+    D: DimName,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
     // 2m choose m = (2m)!/(m! * (2m-m)!)
 
@@ -357,10 +352,10 @@ where
         factorial(2 * m as u128) / (factorial(m as u128) * factorial(2 * m as u128 - m as u128));
     let abs_c_recip = choose_2m_m * factorial(2 * m as u128 + 1);
     let alpha = a_abs_onenorm / one_norm(a);
-    let alpha = alpha / N::from_u128(abs_c_recip).unwrap();
+    let alpha: f64 = try_convert(alpha).unwrap() / abs_c_recip as f64;
 
-    let u = N::from_f64(2_f64.powf(-53.0)).unwrap();
-    let log2_alpha_div_u = try_convert((alpha / u).log2()).unwrap();
+    let u = 2_f64.powf(-53.0);
+    let log2_alpha_div_u = (alpha / u).log2();
     let value = (log2_alpha_div_u / (2.0 * m as f64)).ceil();
     if value > 0.0 {
         value as u64
@@ -369,11 +364,11 @@ where
     }
 }
 
-fn solve_p_q<N, R>(u: MatrixN<N, R>, v: MatrixN<N, R>) -> MatrixN<N, R>
+fn solve_p_q<N, D>(u: MatrixN<N, D>, v: MatrixN<N, D>) -> MatrixN<N, D>
 where
     N: ComplexField,
-    R: DimMin<R, Output = R> + DimName,
-    DefaultAllocator: Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>>,
+    D: DimMin<D, Output = D> + DimName,
+    DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     let p = &u + &v;
     let q = &v - &u;
@@ -381,11 +376,11 @@ where
     q.lu().solve(&p).unwrap()
 }
 
-fn one_norm<N, R>(m: &MatrixN<N, R>) -> N
+fn one_norm<N, D>(m: &MatrixN<N, D>) -> N
 where
     N: RealField,
-    R: DimName,
-    DefaultAllocator: Allocator<N, R, R>,
+    D: DimName,
+    DefaultAllocator: Allocator<N, D, D>,
 {
     let mut col_sums = vec![N::zero(); m.ncols()];
     for i in 0..m.ncols() {
@@ -399,11 +394,11 @@ where
     max
 }
 
-impl<N: RealField, R> MatrixN<N, R>
+impl<N: RealField, D> MatrixN<N, D>
 where
-    R: DimMin<R, Output = R> + DimName,
+    D: DimMin<D, Output = D> + DimName,
     DefaultAllocator:
-        Allocator<N, R, R> + Allocator<(usize, usize), DimMinimum<R, R>> + Allocator<N, R>,
+        Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>> + Allocator<N, D>,
 {
     /// Computes exponential of this matrix
     pub fn exp(&self) -> Self {
@@ -415,30 +410,30 @@ where
         let mut h = ExpmPadeHelper::new(self.clone(), true);
 
         let eta_1 = N::max(h.d4_loose(), h.d6_loose());
-        if eta_1 < N::from_f64(1.495585217958292e-002).unwrap() && ell(&h.a, 3) == 0 {
+        if eta_1 < convert(1.495585217958292e-002) && ell(&h.a, 3) == 0 {
             let (u, v) = h.pade3();
             return solve_p_q(u, v);
         }
 
         let eta_2 = N::max(h.d4_tight(), h.d6_loose());
-        if eta_2 < N::from_f64(2.539398330063230e-001).unwrap() && ell(&h.a, 5) == 0 {
+        if eta_2 < convert(2.539398330063230e-001) && ell(&h.a, 5) == 0 {
             let (u, v) = h.pade5();
             return solve_p_q(u, v);
         }
 
         let eta_3 = N::max(h.d6_tight(), h.d8_loose());
-        if eta_3 < N::from_f64(9.504178996162932e-001).unwrap() && ell(&h.a, 7) == 0 {
+        if eta_3 < convert(9.504178996162932e-001) && ell(&h.a, 7) == 0 {
             let (u, v) = h.pade7();
             return solve_p_q(u, v);
         }
-        if eta_3 < N::from_f64(2.097847961257068e+000).unwrap() && ell(&h.a, 9) == 0 {
+        if eta_3 < convert(2.097847961257068e+000) && ell(&h.a, 9) == 0 {
             let (u, v) = h.pade9();
             return solve_p_q(u, v);
         }
 
         let eta_4 = N::max(h.d8_loose(), h.d10_loose());
         let eta_5 = N::min(eta_3, eta_4);
-        let theta_13 = N::from_f64(4.25).unwrap();
+        let theta_13 = convert(4.25);
 
         let mut s = if eta_5 == N::zero() {
             0
@@ -452,10 +447,7 @@ where
             }
         };
 
-        s += ell(
-            &(&h.a * N::from_f64(2.0_f64.powf(-(s as f64))).unwrap()),
-            13,
-        );
+        s += ell(&(&h.a * convert::<f64, N>(2.0_f64.powf(-(s as f64)))), 13);
 
         let (u, v) = h.pade13_scaled(s);
         let mut x = solve_p_q(u, v);

From e914afe2afd243b568b2d8c41e09a621f465851b Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Sun, 12 Apr 2020 11:59:06 +0200
Subject: [PATCH 09/10] Added support for dynamic matrices

---
 src/linalg/exp.rs   | 23 +++++++++++++----------
 tests/linalg/exp.rs | 30 +++++++++++++++++++++++++++++-
 2 files changed, 42 insertions(+), 11 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index f68cbe6c..e07ddd2f 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -3,7 +3,8 @@
 use crate::{
     base::{
         allocator::Allocator,
-        dimension::{DimMin, DimMinimum, DimName},
+        dimension::{Dim, DimMin, DimMinimum, U1},
+        storage::Storage,
         DefaultAllocator,
     },
     convert, try_convert, ComplexField, MatrixN, RealField,
@@ -13,7 +14,7 @@ use crate::{
 struct ExpmPadeHelper<N, D>
 where
     N: RealField,
-    D: DimName + DimMin<D>,
+    D: DimMin<D>,
     DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     use_exact_norm: bool,
@@ -40,13 +41,14 @@ where
 impl<N, D> ExpmPadeHelper<N, D>
 where
     N: RealField,
-    D: DimName + DimMin<D>,
+    D: DimMin<D>,
     DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     fn new(a: MatrixN<N, D>, use_exact_norm: bool) -> Self {
+        let (nrows, ncols) = a.data.shape();
         ExpmPadeHelper {
             use_exact_norm,
-            ident: MatrixN::<N, D>::identity(),
+            ident: MatrixN::<N, D>::identity_generic(nrows, ncols),
             a,
             a2: None,
             a4: None,
@@ -321,10 +323,11 @@ fn factorial(n: u128) -> u128 {
 fn onenorm_matrix_power_nonm<N, D>(a: &MatrixN<N, D>, p: u64) -> N
 where
     N: RealField,
-    D: DimName,
+    D: Dim,
     DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
-    let mut v = crate::VectorN::<N, D>::repeat(convert(1.0));
+    let nrows = a.data.shape().0;
+    let mut v = crate::VectorN::<N, D>::repeat_generic(nrows, U1, convert(1.0));
     let m = a.transpose();
 
     for _ in 0..p {
@@ -337,7 +340,7 @@ where
 fn ell<N, D>(a: &MatrixN<N, D>, m: u64) -> u64
 where
     N: RealField,
-    D: DimName,
+    D: Dim,
     DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
 {
     // 2m choose m = (2m)!/(m! * (2m-m)!)
@@ -367,7 +370,7 @@ where
 fn solve_p_q<N, D>(u: MatrixN<N, D>, v: MatrixN<N, D>) -> MatrixN<N, D>
 where
     N: ComplexField,
-    D: DimMin<D, Output = D> + DimName,
+    D: DimMin<D, Output = D>,
     DefaultAllocator: Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>>,
 {
     let p = &u + &v;
@@ -379,7 +382,7 @@ where
 fn one_norm<N, D>(m: &MatrixN<N, D>) -> N
 where
     N: RealField,
-    D: DimName,
+    D: Dim,
     DefaultAllocator: Allocator<N, D, D>,
 {
     let mut col_sums = vec![N::zero(); m.ncols()];
@@ -396,7 +399,7 @@ where
 
 impl<N: RealField, D> MatrixN<N, D>
 where
-    D: DimMin<D, Output = D> + DimName,
+    D: DimMin<D, Output = D>,
     DefaultAllocator:
         Allocator<N, D, D> + Allocator<(usize, usize), DimMinimum<D, D>> + Allocator<N, D>,
 {
diff --git a/tests/linalg/exp.rs b/tests/linalg/exp.rs
index 37abc4d2..75122107 100644
--- a/tests/linalg/exp.rs
+++ b/tests/linalg/exp.rs
@@ -2,7 +2,7 @@
 mod tests {
     //https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/sparse/linalg/tests/test_matfuncs.py
     #[test]
-    fn exp() {
+    fn exp_static() {
         use nalgebra::{Matrix1, Matrix2, Matrix3};
 
         {
@@ -98,4 +98,32 @@ mod tests {
             assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
         }
     }
+
+    #[test]
+    fn exp_dynamic() {
+        use nalgebra::DMatrix;
+
+        let m = DMatrix::from_row_slice(3, 3, &[1.0, 3.0, 0.0, 0.0, 1.0, 5.0, 0.0, 0.0, 2.0]);
+
+        let e1 = 1.0_f64.exp();
+        let e2 = 2.0_f64.exp();
+
+        let f = DMatrix::from_row_slice(
+            3,
+            3,
+            &[
+                e1,
+                3.0 * e1,
+                15.0 * (e2 - 2.0 * e1),
+                0.0,
+                e1,
+                5.0 * (e2 - e1),
+                0.0,
+                0.0,
+                e2,
+            ],
+        );
+
+        assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
+    }
 }

From 583a8fb1102f075e4457208f3c1cf68b21707384 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fredrik.jansson.se@gmail.com>
Date: Sun, 12 Apr 2020 13:27:11 +0200
Subject: [PATCH 10/10] Removed unsafe immutable hack

---
 src/linalg/exp.rs | 146 ++++++++++++++++++++++++++--------------------
 1 file changed, 83 insertions(+), 63 deletions(-)

diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index e07ddd2f..86bd81a7 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -66,87 +66,78 @@ where
         }
     }
 
-    fn a2(&self) -> &MatrixN<N, D> {
+    fn calc_a2(&mut self) {
         if self.a2.is_none() {
-            let ap = &self.a2 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(&self.a * &self.a);
-            }
+            self.a2 = Some(&self.a * &self.a);
         }
-        self.a2.as_ref().unwrap()
     }
 
-    fn a4(&self) -> &MatrixN<N, D> {
+    fn calc_a4(&mut self) {
         if self.a4.is_none() {
-            let ap = &self.a4 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            let a2 = self.a2();
-            unsafe {
-                *ap = Some(a2 * a2);
-            }
+            self.calc_a2();
+            let a2 = self.a2.as_ref().unwrap();
+            self.a4 = Some(a2 * a2);
         }
-        self.a4.as_ref().unwrap()
     }
 
-    fn a6(&self) -> &MatrixN<N, D> {
+    fn calc_a6(&mut self) {
         if self.a6.is_none() {
-            let a2 = self.a2();
-            let a4 = self.a4();
-            let ap = &self.a6 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(a4 * a2);
-            }
+            self.calc_a2();
+            self.calc_a4();
+            let a2 = self.a2.as_ref().unwrap();
+            let a4 = self.a4.as_ref().unwrap();
+            self.a6 = Some(a4 * a2);
         }
-        self.a6.as_ref().unwrap()
     }
 
-    fn a8(&self) -> &MatrixN<N, D> {
+    fn calc_a8(&mut self) {
         if self.a8.is_none() {
-            let a2 = self.a2();
-            let a6 = self.a6();
-            let ap = &self.a8 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(a6 * a2);
-            }
+            self.calc_a2();
+            self.calc_a6();
+            let a2 = self.a2.as_ref().unwrap();
+            let a6 = self.a6.as_ref().unwrap();
+            self.a8 = Some(a6 * a2);
         }
-        self.a8.as_ref().unwrap()
     }
 
-    fn a10(&mut self) -> &MatrixN<N, D> {
+    fn calc_a10(&mut self) {
         if self.a10.is_none() {
-            let a4 = self.a4();
-            let a6 = self.a6();
-            let ap = &self.a10 as *const Option<MatrixN<N, D>> as *mut Option<MatrixN<N, D>>;
-            unsafe {
-                *ap = Some(a6 * a4);
-            }
+            self.calc_a4();
+            self.calc_a6();
+            let a4 = self.a4.as_ref().unwrap();
+            let a6 = self.a6.as_ref().unwrap();
+            self.a10 = Some(a6 * a4);
         }
-        self.a10.as_ref().unwrap()
     }
 
     fn d4_tight(&mut self) -> N {
         if self.d4_exact.is_none() {
-            self.d4_exact = Some(one_norm(self.a4()).powf(convert(0.25)));
+            self.calc_a4();
+            self.d4_exact = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25)));
         }
         self.d4_exact.unwrap()
     }
 
     fn d6_tight(&mut self) -> N {
         if self.d6_exact.is_none() {
-            self.d6_exact = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
+            self.calc_a6();
+            self.d6_exact = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0)));
         }
         self.d6_exact.unwrap()
     }
 
     fn d8_tight(&mut self) -> N {
         if self.d8_exact.is_none() {
-            self.d8_exact = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
+            self.calc_a8();
+            self.d8_exact = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0)));
         }
         self.d8_exact.unwrap()
     }
 
     fn d10_tight(&mut self) -> N {
         if self.d10_exact.is_none() {
-            self.d10_exact = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
+            self.calc_a10();
+            self.d10_exact = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0)));
         }
         self.d10_exact.unwrap()
     }
@@ -161,7 +152,8 @@ where
         }
 
         if self.d4_approx.is_none() {
-            self.d4_approx = Some(one_norm(self.a4()).powf(convert(0.25)));
+            self.calc_a4();
+            self.d4_approx = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25)));
         }
 
         self.d4_approx.unwrap()
@@ -177,7 +169,8 @@ where
         }
 
         if self.d6_approx.is_none() {
-            self.d6_approx = Some(one_norm(self.a6()).powf(convert(1.0 / 6.0)));
+            self.calc_a6();
+            self.d6_approx = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0)));
         }
 
         self.d6_approx.unwrap()
@@ -193,7 +186,8 @@ where
         }
 
         if self.d8_approx.is_none() {
-            self.d8_approx = Some(one_norm(self.a8()).powf(convert(1.0 / 8.0)));
+            self.calc_a8();
+            self.d8_approx = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0)));
         }
 
         self.d8_approx.unwrap()
@@ -209,7 +203,8 @@ where
         }
 
         if self.d10_approx.is_none() {
-            self.d10_approx = Some(one_norm(self.a10()).powf(convert(1.0 / 10.0)));
+            self.calc_a10();
+            self.d10_approx = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0)));
         }
 
         self.d10_approx.unwrap()
@@ -217,8 +212,10 @@ where
 
     fn pade3(&mut self) -> (MatrixN<N, D>, MatrixN<N, D>) {
         let b: [N; 4] = [convert(120.0), convert(60.0), convert(12.0), convert(1.0)];
-        let u = &self.a * (self.a2() * b[3] + &self.ident * b[1]);
-        let v = self.a2() * b[2] + &self.ident * b[0];
+        self.calc_a2();
+        let a2 = self.a2.as_ref().unwrap();
+        let u = &self.a * (a2 * b[3] + &self.ident * b[1]);
+        let v = a2 * b[2] + &self.ident * b[0];
         (u, v)
     }
 
@@ -231,8 +228,15 @@ where
             convert(30.0),
             convert(1.0),
         ];
-        let u = &self.a * (self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
-        let v = self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        self.calc_a2();
+        self.calc_a6();
+        let u = &self.a
+            * (self.a4.as_ref().unwrap() * b[5]
+                + self.a2.as_ref().unwrap() * b[3]
+                + &self.ident * b[1]);
+        let v = self.a4.as_ref().unwrap() * b[4]
+            + self.a2.as_ref().unwrap() * b[2]
+            + &self.ident * b[0];
         (u, v)
     }
 
@@ -247,9 +251,18 @@ where
             convert(56.0),
             convert(1.0),
         ];
-        let u =
-            &self.a * (self.a6() * b[7] + self.a4() * b[5] + self.a2() * b[3] + &self.ident * b[1]);
-        let v = self.a6() * b[6] + self.a4() * b[4] + self.a2() * b[2] + &self.ident * b[0];
+        self.calc_a2();
+        self.calc_a4();
+        self.calc_a6();
+        let u = &self.a
+            * (self.a6.as_ref().unwrap() * b[7]
+                + self.a4.as_ref().unwrap() * b[5]
+                + self.a2.as_ref().unwrap() * b[3]
+                + &self.ident * b[1]);
+        let v = self.a6.as_ref().unwrap() * b[6]
+            + self.a4.as_ref().unwrap() * b[4]
+            + self.a2.as_ref().unwrap() * b[2]
+            + &self.ident * b[0];
         (u, v)
     }
 
@@ -266,16 +279,20 @@ where
             convert(90.0),
             convert(1.0),
         ];
+        self.calc_a2();
+        self.calc_a4();
+        self.calc_a6();
+        self.calc_a8();
         let u = &self.a
-            * (self.a8() * b[9]
-                + self.a6() * b[7]
-                + self.a4() * b[5]
-                + self.a2() * b[3]
+            * (self.a8.as_ref().unwrap() * b[9]
+                + self.a6.as_ref().unwrap() * b[7]
+                + self.a4.as_ref().unwrap() * b[5]
+                + self.a2.as_ref().unwrap() * b[3]
                 + &self.ident * b[1]);
-        let v = self.a8() * b[8]
-            + self.a6() * b[6]
-            + self.a4() * b[4]
-            + self.a2() * b[2]
+        let v = self.a8.as_ref().unwrap() * b[8]
+            + self.a6.as_ref().unwrap() * b[6]
+            + self.a4.as_ref().unwrap() * b[4]
+            + self.a2.as_ref().unwrap() * b[2]
             + &self.ident * b[0];
         (u, v)
     }
@@ -300,9 +317,12 @@ where
         let s = s as f64;
 
         let mb = &self.a * convert::<f64, N>(2.0_f64.powf(-s));
-        let mb2 = self.a2() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
-        let mb4 = self.a4() * convert::<f64, N>(2.0.powf(-4.0 * s));
-        let mb6 = self.a6() * convert::<f64, N>(2.0.powf(-6.0 * s));
+        self.calc_a2();
+        self.calc_a4();
+        self.calc_a6();
+        let mb2 = self.a2.as_ref().unwrap() * convert::<f64, N>(2.0_f64.powf(-2.0 * s));
+        let mb4 = self.a4.as_ref().unwrap() * convert::<f64, N>(2.0.powf(-4.0 * s));
+        let mb6 = self.a6.as_ref().unwrap() * convert::<f64, N>(2.0.powf(-6.0 * s));
 
         let u2 = &mb6 * (&mb6 * b[13] + &mb4 * b[11] + &mb2 * b[9]);
         let u = &mb * (&u2 + &mb6 * b[7] + &mb4 * b[5] + &mb2 * b[3] + &self.ident * b[1]);