Thanks…but it seems that this one only sticks to Float64, right? it might use polynomial of Float64 internal?
julia> A = SMatrix{3, 3, Float32}(rand(Float32, 3, 3))
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
0.276386 0.954482 0.714192
0.574739 0.3517 0.418869
0.555151 0.3695 0.359785
julia> exponential!(A, ExpMethodGeneric())
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
2.11407 1.78511 1.49504
1.13127 2.05394 1.03571
1.08678 1.04209 1.95143
julia> ForwardDiff.jacobian(x -> exponential!(x, ExpMethodGeneric()), A)
9×9 SMatrix{9, 9, Float64, 81} with indices SOneTo(9)×SOneTo(9):
1.8236 0.838051 0.680035 0.523368 … 0.503557 0.158246 0.126138
0.523368 1.82291 0.451652 0.101791 0.0980227 0.506292 0.0816883
0.503557 0.441118 1.79091 0.0980227 0.0943939 0.0782926 0.501925
0.838051 0.265264 0.211498 1.82291 0.441118 0.131362 0.104497
0.164332 0.842565 0.137019 0.526213 0.0812944 0.443656 0.0674744
0.158246 0.131362 0.835204 0.506292 … 0.0782926 0.0645242 0.440187
0.680035 0.211498 0.16852 0.451652 1.79091 0.835204 0.677983
0.130985 0.683771 0.109075 0.0848232 0.521681 1.7899 0.450526
0.126138 0.104497 0.677983 0.0816883 0.501925 0.440187 1.75782