From a90c2172a090b6dab6926eaeb0bec273d09921c1 Mon Sep 17 00:00:00 2001
From: vera-rubin <322284620+vera-rubin@users.noreply.github.com>
Date: Thu, 15 Oct 2026 06:32:38 +0000
Subject: [PATCH] Avoid overflow in gamma_inc derivatives

---
 src/dual.jl      |  8 +++++++-
 test/DualTest.jl | 12 ++++++++++++
 2 files changed, 19 insertions(+), 1 deletion(-)

diff --git a/src/dual.jl b/src/dual.jl
index 6d13dec3..fc4a249f 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -878,7 +878,13 @@ end
 function SpecialFunctions.gamma_inc(a::Real, d::Dual{T,<:Real}, ind::Integer) where {T}
     x = value(d)
     p, q = SpecialFunctions.gamma_inc(a, x, ind)
-    ∂p = exp(-x) * x^(a - 1) / SpecialFunctions.gamma(a) * partials(d)
+    logx = log(x)
+    if isfinite(logx)
+        ∂p = exp(-x + (a - 1) * logx - SpecialFunctions.loggamma(a)) * partials(d)
+    else
+        # Preserve the power rule at zero for higher-order derivatives.
+        ∂p = exp(-x) * x^(a - 1) / SpecialFunctions.gamma(a) * partials(d)
+    end
     return (Dual{T}(p, ∂p), Dual{T}(q, -∂p))
 end
 
diff --git a/test/DualTest.jl b/test/DualTest.jl
index 67d9c9f7..3e932976 100644
--- a/test/DualTest.jl
+++ b/test/DualTest.jl
@@ -667,6 +667,18 @@ end
 # bug fixes #
 #############
 
+@testset "gamma_inc derivatives" begin
+    for V in (Float32, Float64), (a, x) in ((150, 150), (200, 200), (1, 0), (2, 0)), ind in 0:2
+        # Compute the reference in arbitrary precision to avoid overflowing intermediates.
+        expected = V(exp(-big(x)) * big(x)^(a - 1) / factorial(big(a - 1)))
+        p, q = gamma_inc(V(a), Dual{TestTag}(V(x), one(V)), ind)
+        @test partials(p, 1) ≈ expected
+        @test partials(q, 1) ≈ -expected
+    end
+    @test ForwardDiff.derivative(x -> ForwardDiff.derivative(y -> gamma_inc(2.0, y, 0)[1], x), 0.0) ≈ 1.0
+    @test ForwardDiff.derivative(x -> ForwardDiff.derivative(y -> gamma_inc(2.0, y, 0)[2], x), 0.0) ≈ -1.0
+end
+
 @testset "Exponentiation of zero" begin
     x0 = 0.0
     x1 = Dual{:t1}(x0, 1.0)