[vera-rubin]

Differentiating gamma_inc(a, x, 0) with respect to x at a = x = 150.0 or 200.0 returns Inf or NaN:

using ForwardDiff, SpecialFunctions
for a in (150.0, 200.0)
    @show ForwardDiff.derivative(x -> gamma_inc(a, x, 0)[1], a)
end

The documented normalization and derivative identity give finite derivatives of about 0.03255540946 and 0.02819772769, respectively.

We evaluate the derivative in log space using loggamma to avoid overflowing intermediate values while retaining the power rule at zero for higher derivatives.

The new regression tests and the existing dual-number tests pass.

Tested versions: ForwardDiff 1.4.6, Julia 1.10.12, SpecialFunctions 2.9.0.

— ChipLab

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.25%. Comparing base (a3c0f4f) to head (a90c217).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #862      +/-   ##
==========================================
+ Coverage   91.23%   91.25%   +0.02%     
==========================================
  Files          11       11              
  Lines        1072     1075       +3     
==========================================
+ Hits          978      981       +3     
  Misses         94       94              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[devmotion]

Thanks for spotting the overflow and for the PR!

Unfortunately, computing the derivative in log space suffers from cancellation in the exponent for large a. Against BigFloat references, the relative error is ~7e-10 at a = x = 1e6 and ~1e-5 at a = x = 1e10 in Float64, and in Float32 ~7e-4 at a = 1e4 and ~16% at a = 1e6.

SpecialFunctions already computes x^a e^(-x) / Γ(a) accurately internally, so I opened JuliaMath/SpecialFunctions.jl#568 to expose the derivative as a public gamma_inc_dx(a, x) (along with helpers for its own derivatives), and #871 to use it here. With it, first and second derivatives are accurate to a few ulps over this range, including at x = 0.

I'd suggest continuing in #871, which also includes a fallback for older SpecialFunctions versions.