Fix exp of SMatrix with large entries #1295 (#1296)

Fix #1295 

Appreciate if you can review this @mikmoore.

Author: juliohm

Project.toml

@@ -1,6 +1,6 @@
 name = "StaticArrays"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.12"
+version = "1.9.13"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/expm.jl

@@ -14,6 +14,7 @@ end
     (newtype)((exp(A[1]), ))
 end
 
+# Bernstein, D. S. & So, W. 1993. "Some Explicit Formulas for the Matrix Exponential"
 @inline function _exp(::Size{(2,2)}, A::StaticMatrix{<:Any,<:Any,<:Real})
     T = typeof(exp(zero(eltype(A))))
     newtype = similar_type(A,T)
@@ -25,25 +26,50 @@ end
 
     v = (a-d)^2 + 4*b*c
 
+    m = (a + d) / 2
+
     if v > 0
-        z = sqrt(v)
-        z1 = cosh(z / 2)
-        z2 = sinh(z / 2) / z
+        # In this case the formulas in the entries of the matrix
+        # are a function of cosh and sinh, and could (in theory)
+        # follow the same code pattern of the other branches (v ≤ 0).
+        # However, cosh and sinh explode with large arguments and
+        # we use the following identity to avoid numerical issues:
+        #
+        # exp(m) * [c₁ * cohs(δ) + c₂ * sinh(δ)] =
+        #   c₁ * (e₊ + e₋) / 2 + c₂ * (e₊ - e₋) / 2
+        #
+        # where e₊ = exp(m + δ) and e₋ = exp(m - δ).
+        #
+        # See https://github.com/JuliaArrays/StaticArrays.jl/issues/1295
+        δ = sqrt(v) / 2
+        e₊ = exp(m + δ)
+        e₋ = exp(m - δ)
+        e₁ = (e₊ + e₋) / 2
+        e₂ = (e₊ - e₋) / 2
+        c₂ = (a - d) / 2δ
+        m11 = (e₁ + c₂ * e₂)
+        m12 = (b / δ) * e₂
+        m21 = (c / δ) * e₂
+        m22 = (e₁ - c₂ * e₂)
     elseif v < 0
         z = sqrt(-v)
+        r = exp(m)
         z1 = cos(z / 2)
         z2 = sin(z / 2) / z
+        m11 = r * (z1 + (a - d) * z2)
+        m12 = r * 2b * z2
+        m21 = r * 2c * z2
+        m22 = r * (z1 - (a - d) * z2)
     else # if v == 0
+        r = exp(m)
         z1 = T(1.0)
         z2 = T(0.5)
+        m11 = r * (z1 + (a - d) * z2)
+        m12 = r * 2b * z2
+        m21 = r * 2c * z2
+        m22 = r * (z1 - (a - d) * z2)
     end
 
-    r = exp((a + d) / 2)
-    m11 = r * (z1 + (a - d) * z2)
-    m12 = r * 2b * z2
-    m21 = r * 2c * z2
-    m22 = r * (z1 - (a - d) * z2)
-
     (newtype)((m11, m21, m12, m22))
 end
 

test/expm.jl

@@ -13,6 +13,10 @@ using StaticArrays, Test, LinearAlgebra
     @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])
 
+    # https://github.com/JuliaArrays/StaticArrays.jl/issues/1295
+    @test exp(@SMatrix [-800.0 800.0; 800.0 -800.0])::SMatrix ≈ [0.5 0.5; 0.5 0.5]
+    @test exp(@SMatrix [-800.0 800.0; 800.0 -800.0]) ≈ exp([-800.0 800.0; 800.0 -800.0])
+
     for sz in (3,4), typ in (Float64, Complex{Float64})
         A = rand(typ, sz, sz)
         nA = norm(A, 1)