Use expm1 for added precision

Author: sethaxen

src/expm.jl

@@ -50,15 +50,15 @@ end
     @inbounds d = A[4]
 
     z = sqrt((a - d)*(a - d) + 4*b*c )
-    e = exp((a + d - z)/2)
-    f = exp((a + d + z)/2)
+    e = expm1((a + d - z) / 2)
+    f = expm1((a + d + z) / 2)
     ϵ = eps()
     g = abs2(z) < ϵ^2 ? exp((a + d) / 2) * (1 + z^2 / 24) : (f - e) / z
 
-    m11 = (g * (a - d) + f + e) / 2
+    m11 = (g * (a - d) + f + e) / 2 + 1
     m12 = g * b
     m21 = g * c
-    m22 = (-g * (a - d) + f + e) / 2
+    m22 = (-g * (a - d) + f + e) / 2 + 1
 
     (newtype)((m11, m21, m12, m22))
 end

test/expm.jl

@@ -8,6 +8,8 @@ using StaticArrays, Test, LinearAlgebra
     @test exp(@SMatrix [ -3+1im  -1+2im;-4-5im   5-2im])::SMatrix ≈ exp(Complex{Float64}[ -3+1im  -1+2im;-4-5im   5-2im])
     # test for case `(a[1,1] - a[2,2])^2 + 4*a[2,1]*a[1,2] == 0`
     @test exp(@SMatrix [2+2im  1-1im; -2-2im  6+2im])::SMatrix ≈ exp(Complex{Float64}[2+2im  1-1im; -2-2im  6+2im])
+    # test for zeros matrix
+    @test exp(@SMatrix zeros(Complex{Float64}, 2, 2))::SMatrix ≈ Complex{Float64}[1 0; 0 1]
     @test exp(@SMatrix [1 2 0; 2 1 0; 0 0 1]) ≈ exp([1 2 0; 2 1 0; 0 0 1])
 
     for sz in (3,4), typ in (Float64, Complex{Float64})