add matrix exponentials

src/Quaternion.jl

@@ -293,3 +293,25 @@ function slerp(qa::Quaternion{T}, qb::Quaternion{T}, t::T) where {T}
         qa.v3 * ratio_a + qm.v3 * ratio_b,
     )
 end
+
+function _complex_representation(A::Matrix{Quaternion{T}}) where {T}
+    # convert a quatenion matrix to corresponding complex matrix
+    a = map(t->t.s, A)
+    b = map(t->t.v1, A)
+    c = map(t->t.v2, A)
+    d = map(t->t.v3, A)
+
+    return [a+im*b c+im*d;-c+im*d a-im*b]
+end
+
+function _quaternion_representation(A::Matrix{Complex{T}}) where {T}
+    n = round(Int, size(A, 1)/2)
+
+    X = @view A[1:n, 1:n]
+    Y = @view A[1:n, n+1:2n]
+
+    return [Quaternion(X[i,j].re, X[i,j].im, Y[i,j].re, Y[i,j].im) for i in 1:n, j in 1:n]
+end
+
+# A quick way of doing the quaternionic matrix exponential
+exp(A::Matrix{Quaternion{T}}) where {T} = _quaternion_representation(exp(_complex_representation(A)))

test/test_Quaternion.jl

@@ -138,4 +138,13 @@ for _ in 1:100
         @test q ⊗ slerp(q1, q2, t) ≈ slerp(q ⊗ q1, q ⊗ q2, t)
         @test q ⊗ linpol(q1, q2, t) ≈ linpol(q ⊗ q1, q ⊗ q2, t)
     end
+    let # test matrix exponential
+        z = zeros(Quaternion{Float64}, 3, 3)
+        id = Matrix{Quaternion{Float64}}(I, 3, 3)
+        d = [Quaternion(randn(4)...) for i in 1:3]
+        expd = exp.(d)
+        D = exp(Array(Diagonal(d)))
+        @test exp(z) ≈ id
+        @test D ≈ Array(Diagonal(expd))
+    end
 end