Numerical calculation of incomplete gamma functions by the trapezoidal rule (Q1076477)

The trapezoidal rule is applied to the numerical calculation of the known integral representation of the complementary incomplete gamma function \[ \Gamma (a,x)=(2e^{-x}x^ a/\Gamma (1-a))\int^{+\infty}_{0}e^{- u^ 2}u^{-2a+1}/x+\quad u^ 2\quad du \] in the region \(a<-1\) and \(x>0\). This application of the trapezoidal rule is not standard since the function \(u^{-2a+1}(x+u^ 2)^{-1}\) in the integrand is not even, but numerical tests show that using the rule may still be convenient. The explanation of this surprising fact requires a careful investigation of the related Euler-Maclaurin formula and expecially of its remainder term. The outcome is a simple numerical procedure for obtaining values of incomplete gamma functions with good accuracy in the stated region. The explained method can be applied to the numerical evaluation of other important special functions.