I am looking to replace the generalised incomplete gamma function (which appears in a solution to a problem I've posted about here) with a numerically stable formula involving other functions. This is because I want to move over to other languages which do not have this function implemented.
To illustrate my point, the generalised incomplete gamma function when expanded is,
Gamma[a, b, c] = Gamma[a, b] - Gamma[a, c]
So naively you might think that the LHS and RHS should return the answer. The RHS however is numerically unclever and so has issues when $a$ is large and $b$ and $c$ are small. For example,
Gamma[100, 0.1, 0.01]
Gamma[100, 0.1] - Gamma[100, 0.01]
My question is what does the generalised incomplete gamma function actually do in the background? And so, how can I replicate this behaviour using simpler functions?
Note 1: I do not want answers explicitly involving calls to NIntegrate or Integrate since these are inefficient to implement. Rather I would like a solution involving other functions.
Note 2: I know the definition of the incomplete Gamma in Mathematica is the upper incomplete gamma. If there were a way to write a numerically stable form using the lower incomplete gamma, that may avoid the issue since the Gamma[a] terms would cancel out.