Preserve input element type in ExpMethodGeneric (fix Float32 → Float64 promotion)

[ChrisRackauckas-Claude]

Please ignore until reviewed by @ChrisRackauckas. Draft opened by an agent.

Problem

Reported on Discourse: How to use ForwardDiff.jl with inv and matrix exponential (post #15).

exponential!(x, ExpMethodGeneric()) silently promotes Float32 inputs to Float64:

julia> A = rand(SMatrix{3,3,Float32});

julia> exponential!(A, ExpMethodGeneric()) |> eltype
Float64            # expected Float32

julia> ForwardDiff.jacobian(x -> exponential!(x, ExpMethodGeneric()), A) |> eltype
Float64            # expected Float32

Cause

The specialized (13,13) Padé numerator methods hardcoded Float64 coefficients — UniformScaling{Float64}(...) for the immutable-matrix path and Float64 literals for the scalar path. The hardcoded Float64 coefficients force the whole @evalpoly Horner evaluation to promote to Float64, regardless of the input element type. (The immutable/SMatrix path is the one that hits these; mutable matrices go through exp_generic_mutable, which already promotes to Float64 by design and is untouched here.)

Fix

Take the coefficient type from float(eltype(x)) (matrix) / float(typeof(x)) (scalar) at runtime, using the exact Padé rationals. For Float64 inputs the resulting coefficients are bit-identical to the previous hardcoded literals (verified), so Float64 results are unchanged. Float32, ComplexF32, Double64, and ForwardDiff Dual inputs now retain their element type.

This is deliberately not @generated: deriving the coefficient type inside a generated-function body requires constructing the element type (e.g. a ForwardDiff Dual), which is world-age blocked and breaks differentiation. That is also why the existing generic @generated exp_pade_p(x, ::Val{k}, ::Val{m}) method (used for non-13 orders) could not simply be reused for the default k = 13 matrix/scalar case.

Tests

Added a testset asserting element-type preservation for ExpMethodGeneric across Float32/Float64 static matrices (result eltype and ForwardDiff jacobian eltype), ComplexF32, and the scalar path. Full suite passes locally on Julia 1.10 (364 pass, 1 pre-existing broken).

🤖 Generated with Claude Code