NaN gradient occurring in `SpecialFunctions.gamma`, when using float arguments.

[vwiela]

I tried to use ForwardDiff on the upper incomplete gamma function as provided in SpecialFunctions.jl

Calling it with the first argument being an integer did work fine, but changing it to a float, caused the derivative being NaN. I tried to dig into it a bit and the following MwE does show the problem.

using ForwardDiff
using SpecialFunctions

function testing_gamma(a::Number, x::Number)
    return iszero(x) ? gamma(one(x)*a) : x^a * expint(1 - a, x)
end


promotereal(x::Real, y::Real) = promote(x, y)
test_gamma(a::Number,x::Number) = testing_gamma(promotereal(float(a), float(x))...)

If we want to get the gradient of testing_gamma wrt. to the second argument, this works fine

ForwardDiff.derivative(x -> testing_gamma(2.0, -x), 0.01)

However, on test_gamma it will return a NaN value.

ForwardDiff.derivative(x -> test_gamma(2.0, -x), 0.01)

Maybe someone has an idea why this behavior occurs and how it can best be fixed.

[vwiela]

@gdalle This is the issue.

[ChrisRackauckas]

See https://juliadiff.org/ForwardDiff.jl/dev/user/advanced/#Fixing-NaN/Inf-Issues

[devmotion]

As suggested in the previous comment, this problem can be fixed by enabling NaN-safe mode. With NaN-safe mode enabled I get:

julia> ForwardDiff.derivative(x -> testing_gamma(2.0, -x), 0.01)
-0.010100501670848416

julia> ForwardDiff.derivative(x -> test_gamma(2.0, -x), 0.01)
-0.010100501670848416