My question is about incomplete Gamma (IncGamma) function,
$$ \Gamma(a,x)=\int_x^\infty t^{a-1}e^{-t} dt. $$ Using basic algebra, $$ \Gamma(a,x)-\Gamma(a,y) = \int_x^y t^{a-1}e^{-t} dt $$ clearly finding the integral between x and y is less complex and includes less numerical instability than two separate IncGamma functions.
With above introduction, my question is about directly computing the difference of two IncGamma functions with very high precision.
Are there any ways in particular in software?
If you can rewrite $y =x + \delta$, such that $\vert \delta \vert < \vert x \vert$, it can be shown that
$$\Gamma(a,x) - \Gamma(a,x+\delta) = e^{-x}x^{a-1}\sum_{n=0}^{\infty} \frac{(1-a)_n}{(-x)^n}(1-e^{-\delta} e_n(\delta))$$
where $(1-a)_n$ is a Pochhammer symbol and $e_n(\delta) = \sum_{k = 0}^n \frac{\delta^k}{k !}$
Take a look in here for reference.