[Largely edited version]

I was looking at this issue about the derivative of the incomplete gamma function w.r.t the parameter `a`

.

Looks like from this page of wolfram alpha docs that the right formula is :

which depends on the regularized hypergeometric function. After a certain time of translation, namely using this page to match the wolfram alpha defintion to the one from the HypegeometricFunctions.jl package (also, the formula is about non-normalized upper incomplete gamma function, while `gamma_inc`

gives the normalized version), I have found a working implementation :

```
using SpecialFunctions, HypergeometricFunctions
function Γ(a,x)
return gamma_inc(a,x)[2]
end
function ∂Γ_∂a(a,x)
g = gamma(a)
dg = digamma(a)
G = Γ(a,x)
lx = log(x)
r = pFq([a,a],[a+1,a+1],-x)
return a^(-2) * x^a * r/g + (1 - G)*(dg - lx)
end
function emp(a,x)
ϵ = 1e-5
return (Γ(a+ϵ,x) - Γ(a-ϵ,x))/2ϵ
end
a = 1.9
x = 1.8
Γ(a,x)
(∂Γ_∂a(a,x), emp(a,x))
#a simple test:
a = 0.1:0.1:5
x = 0.1:0.1:5
diff(a,x) = abs(∂Γ_∂a(a,x) - emp(a,x))
maximum(diff.(a,x')) # 2e-10 on my machine.
```

**How could I:**

- Check in a better way that this is yielding the right derivative ?
- Tell ForwardDiff to use this instead of erroring out on
`ForwardDiff.gradient(a -> gamma_inc(a,1)[2],1.2)`

.