4

N. M. Temme, "Special Functions" (Wiley 1996) gives the following expression that expresses the upper incomplete gamma function in terms of the ordinary gamma function, for integer orders:

$$ \Gamma(n,z) = \Gamma(n) e^{-z} \sum_{m=0}^{n-1} \frac{z^m}{m!} \\ n=1,2,... $$

My experiments indicate that this is a convenient way to compute the upper incomplete gamma function for small integer orders as the computation appears to be numerically stable. I tried orders up to n=50 and a wide range of real z.

Is there a similar expression that allows the straightforward computation of the upper incomplete gamma function in terms of the ordinary gamma function, for half-integer orders, that is, $ \Gamma(n+\frac{1}{2},z) $? I am aware that $ \Gamma(\frac{1}{2},z) = \sqrt{\pi} erfc(\sqrt{z}) $.

njuffa
  • 1,768
  • Apart from your question. How it csn be said that the above relation is numerically stable? As, the R.H.S also involves $\Gamma(n)$. – kaka Mar 24 '14 at 01:50
  • C++ provides a function tgamma() that computes $ \Gamma(x) $ accurately, if a good math library is provided. The remaining computation does not involve subtractive cancellation and excessive round-off error, so a straight rendering of the above into C++ code (for n <= 50) seems to be accurate to about 15 decimal digits if the entire computation is performed in double precision. I am not a numerical analyst, just a computer science guy, so I determined the stability empirically (which is why I stated that it "appears" to be stable). I took care to avoid intermediate underflow in exp(-z). – njuffa Mar 24 '14 at 03:15

1 Answers1

2

According to Maple, for nonnegative integers $n$ $$ \Gamma(1/2+n, t) = \text{pochhammer}(1/2,n) \sqrt{\pi}\; \text{erfc}(\sqrt{t}) + t^{n-1/2} e^{-t} \sum_{k=0}^{n-1} \text{pochhammer}(1/2-n,k) (-t)^{-k} $$

Note that $\text{pochhammer}(1/2,k) = 2^{1-2k} \dfrac{(2k-1)!}{(k-1)!}$.

Robert Israel
  • 448,999
  • Could you please clarify the definition of $pochhammer (\frac{1}{2}-n, k) $ ? I am not familiar with this and can't determine it from the above (there is n on the left-hand side but k on the right-hand side, and I have been unable to track it down in literature so far. Thanks. – njuffa Mar 24 '14 at 01:40
  • 1
    For any nonnegative integer $k$, $\text{pochhammer}(z,k) = \prod_{j=0}^{k-1} (z + j) = z (z+1) \ldots (z + k-1)$ (with $\text{pochhammer}(z,0) = 1$). This is $\Gamma(z+k)/\Gamma(z)$ if $z$ is not a nonpositive integer. – Robert Israel Mar 24 '14 at 05:53
  • Thanks for the clarification. With $ pochhammer(z,k) = \Gamma(z+k) / \Gamma(z) $ the straightforward translation into C++ works fine. – njuffa Mar 24 '14 at 07:16
  • In a now deleted answer, a poster pointed out this reference: M.A. Chaudhry, N.M. Temme, E.J.M. Veling, "Asymptotics and closed form of a generalized incomplete gamma function", Journal of Computational and Applied Mathematics 67 (1996) 371-379. On page 379, the following formula is given, which seems to be a more concise variant of the above: $$\Gamma(\alpha,x) = \Gamma(\alpha)\left{erfc(\sqrt{x})+e^{-x} \sum_{j=0}^{\alpha-\frac{3}{2}}\frac{x^{j+\frac{1}{2}}}{(j+\frac{1}{2})\Gamma({j+\frac{1}{2}})}\right}, \space \space \space \alpha = \frac{3}{2},\frac{5}{2}, ...,$$ – njuffa Aug 14 '15 at 16:29