`\(::SMatrix{3,3}, ::SVector{3})` is prone to under/overflow

Julia v1.11.1 with StaticArrays v1.9.8.

Discovered while attempting

julia> using StaticArrays

julia> M = @SMatrix [0 0 1f1; 0 0 0; 0 0 0]

3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
0.0  0.0  10.0
0.0  0.0   0.0
0.0  0.0   0.0

julia> exp(M)
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
NaN  NaN  NaN
NaN  NaN  NaN
NaN  NaN  NaN

This is a duplicate of #785 up to this point, but some investigation reveals the underlying issue, which I discuss here.

The problematic calculation is in _exp where

julia> U = @SMatrix Float32[0.0 0.0 1.6191188f17; 0.0 0.0 0.0; 0.0 0.0 0.0]
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
0.0  0.0  1.61912f17
0.0  0.0  0.0
0.0  0.0  0.0

julia> V = @SMatrix Float32[6.4764752f16 0.0 0.0; 0.0 6.4764752f16 0.0; 0.0 0.0 6.4764752f16]
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
6.47648f16  0.0         0.0
0.0         6.47648f16  0.0
0.0         0.0         6.47648f16

julia> VmU = V - U; VpU = V + U;

julia> VmU \ VpU
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
   NaN  0.0  NaN
   0.0  NaN  0.0
   0.0  0.0  NaN

This can be resolved by

julia> lu(VmU) \ (VpU)
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
1.0  0.0  5.0
0.0  1.0  0.0
0.0  0.0  1.0

The issue is with _solve(::Size{(3,3)}, ::Size{(3,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}. The dependence of _solve on det makes it vulnerable to under/overflow. We can replicate the issue in Float64

julia> @SMatrix([1e200 0 0; 0 1e200 0; 0 0 1e200]) \ @SVector([1e0,0,0])
3-element SVector{3, Float64} with indices SOneTo(3):
NaN
NaN
NaN

This issue is closely related to #959, but that issue has to do with numerical instability and not overflow. Although it is likely that finding a resolution to that would improve the situation here. This issue is most acute in the 3x3 case, but also affects the 2x2 specialization.

I see that _solve is much faster than lu

julia> using BenchmarkTools

julia> @btime \($VmU, $VpU);
  5.466 ns (0 allocations: 0 bytes)

julia> @btime \(lu($VmU), $VpU);
  48.222 ns (0 allocations: 0 bytes)

but to what extent is this worth it? Is there a way we can keep the performance of the closed-form inverse without overflow? For example, can we scale the values in an attempt to keep det from overflowing? Can lu be made faster? At the very least, it seems that using lu in _exp would close #785.

Start by inspecting the _solve(::Size{(3,3)}, ::Size{(3,)}, ...) specialization and how _exp calls it. Reproduce the Float32 and Float64 examples, then compare the closed-form path with lu. Done means the affected solves no longer produce NaN from determinant overflow while the relevant behavior and performance tradeoff are addressed.