[devmotion]

This PR adds gamma_inc_dx(a, x), a public (not exported) function for the derivative of the regularized incomplete gamma functions with respect to x:

$$\frac{\partial P(a, x)}{\partial x} = - \frac{\partial Q(a, x)}{\partial x} = \frac{x^{a - 1} e^{-x}}{\Gamma(a)}$$

For a > 0 it is the density of the gamma distribution with shape a and unit scale.

The PR also adds a ChainRules rule for gamma_inc_dx with derivatives in both arguments, and uses gamma_inc_dx to fix the ChainRules rules of gamma_inc and loggamma(a, x).

Motivation

The naive formula exp(-x) * x^(a - 1) / gamma(a) returns NaN or 0 once an intermediate result over- or underflows, e.g. a = x = 150, even though the result is perfectly representable. The log-space form avoids this but suffers from cancellation for large a: its relative error is 7e-10 at a = 1e6 and 1e-5 at a = 1e10 in Float64.

This came up in JuliaDiff/ForwardDiff.jl#862. ForwardDiff will call gamma_inc_dx when it is defined.

StatsFuns can reuse it too:

• StatsFuns.gammapdf is computed as exp(gammalogpdf(...)), the log form above. gamma_inc_dx(k, x / θ) / θ would be accurate for large shapes.
• The derivatives of the gamma tails in Differentiate the gamma tail functions under ForwardDiff JuliaStats/StatsFuns.jl#231 currently use the same log form.

Since gammapdf is commonly differentiated with respect to the shape, e.g. to fit k by maximum likelihood, gamma_inc_dx needs an accurate derivative in a as well as in x.

gamma_inc_dx

• The direct formula is used whenever all intermediate results are guaranteed to be normal and finite. Its accuracy is about 1 ulp. The guarantee is checked a priori from the arguments (x < 708, -170 < a < 171, |(a - 1) * log2(x)| < 1022), not by computing and discarding values.
• Otherwise it uses rgammax(a, x) / x (for a ≥ -0.5), which is accurate for large a and x, or the log form for very negative a.
• Float16 and Float32 are computed in Float64. BigFloat always uses the direct formula.
• The special values at x = 0, x = Inf, non-positive integer a, and NaN follow the limits of the formula.

Accuracy against 512-bit references:

• Float64: within 8 ulps outside the fallback regimes; the fallbacks are within the conditioning of the problem.
• Float32: within 0.5 ulps.
• BigFloat (150, 150): exact.

Timings:

(a, x)	before (naive formula)	gamma_inc_dx
(2.5, 1)	22.6 ns	27 ns
(30, 25)	58 ns	60 ns
(150, 150)	107 ns (wrong result)	30 ns
(1e6, 1e6)	wrong result	23.5 ns

ChainRules

• gamma_inc: the rule uses gamma_inc_dx.
• gamma_inc_dx (new rule): this gives ChainRules-based AD accurate second derivatives of gamma_inc in x, and makes the new function differentiable in both arguments. The derivatives are computed by two internal helpers that ForwardDiff can share:
  • ∂f/∂x = f · (a − 1 − x)/x. a - 1 - x is exact near the peak, where (a - 1)/x - 1 would cancel. Only when f underflows does it switch to gamma_inc_dx(a - 1, x) - f, which is also exact at x = 0 and x = ∞.
  • ∂f/∂a = f · (log x − ψ(a)). For a ≥ 10 the factor is computed as log1p((x - a)/a) + (log a − ψ(a)), with the asymptotic expansion of log a − ψ(a). This avoids the cancellation of log(x) - digamma(a), which loses 6e6 eps at a = x = 1e6 and 8e10 eps at a = x = 1e10; the new form is within 1 eps there.
  • Poles of ψ(a) (a = 0, −1, −2, …, where f = 0): the reflection formula ψ(a) = ψ(1 − a) − π cot(πa) gives a form that is smooth in a, so nested derivatives are correct there too.
  • Accuracy: both partials are within 8 eps of 512-bit references in Float32 and Float64, at large a (150, 1e6, 1e10), near the peak, and at the poles.
  • Timings: ∂x takes 3.5 ns and ∂a 8 to 17 ns, against 24 to 62 ns for gamma_inc_dx itself.
  • Types: Float16 and Float32 are computed in Float64.
• loggamma(a, x):
  • The old rule -exp(-(x + Ω)) * x^(a - 1) cancels. It has a relative error of 6e-7 at (2, 1e10) and 5e-11 at (30, 1e6).
  • Since the primal is computed as a*log(x) + log(expintx(1 - a, x)) - x, the derivative is -1 / (x * expintx(1 - a, x)) with no cancellation. Its relative error is ≤ 4 eps in the tested cases.
  • The frule/rrule are now hand-written so that the primal and the derivative share expintx.
  • The plain loggamma(a, x) primal is unchanged and no slower.

Other changes

loggamma(a, x) now promotes with map(float, promotereal(a, x)).

• loggamma(2, 3f0) stays Float32.
• This makes the separate loggamma(a::Integer, x) method redundant, so it is removed.
• The odd-integer branch for negative x uses isinteger(a) && isodd(a).

The asymptotic expansion at the end of digamma is factored out into _log_minus_digamma_asymptotic, so that ∂f/∂a can reuse it. The final step is now evaluated as log(z) - muladd(t², P(t²), t/2). This changes about 25% of results for z ≥ 8 in the last bit, all toward the correctly rounded value. Over 200 000 points, the maximum error drops from 1.59 to 1.32 ulps.

The version is bumped to 2.10.0.

Notes

#550 touches the same loggamma(a, x) code (removing the Integer method and the negative-x branch). The two PRs will conflict and should be coordinated.

ForwardDiff's Dual rules will call the internal helpers _gamma_inc_dx_dx and _gamma_inc_dx_da, behind an isdefined check. Nested Duals pass a Dual f to _gamma_inc_dx_dx, which calls floatmin(f), so that relies on JuliaDiff/ForwardDiff.jl#855.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 97.43590% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.66%. Comparing base (d8ffcc1) to head (971dd1a).

Files with missing lines	Patch %	Lines
src/gamma_inc.jl	96.66%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #568      +/-   ##
==========================================
+ Coverage   94.61%   94.66%   +0.05%     
==========================================
  Files          14       14              
  Lines        3025     3094      +69     
==========================================
+ Hits         2862     2929      +67     
- Misses        163      165       +2     

Flag	Coverage Δ
unittests	94.66% <97.43%> (+0.05%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ 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.