In a problem involving the (lower) incomplete gamma function $\gamma(a,x)$, with $a>0$ and $x<0$, I've decided to apply the following equivalence: $$ \gamma(a,x) = e^{-x}x^a \sum_{k=0}^{+\infty}{ \frac{ x^k }{\, a (a+1) \cdots (a+k) \,} } \, \mathrm{,} $$ which can be found here, for instance (and is valid for negative values of $x$).

The first two terms cancel with other terms in my equations, and I finally keep with this series: $$ \sum_{k=0}^{+\infty}{ \frac{ x^k }{\, a (a+1) \cdots (a+k) \,} } \, \mathrm{.} $$

The problem is that, when trying to estimate the value of that series, I get results that are not accurate enough.

How can I deal with this series? How can the result of this series be accurately estimated?

For instance, for $a=3/2$ and $x=-100/3$, the exact result of the series (according to Mathematica) is $0.02954292068693\dots$, while a computer routine I've implemented in another environment different from Mathematica returns an approximated value of $0.0295481649038644$ (an absolute difference of $0.000299250391986472$); another implementation that I've tried returns $0.0295299691304212$, for instance. Those small differences become bigger when the result is used in further calculations.

It is not a matter of the stopping criteria of the routines I've programmed. The returned values do not change if I add more iterations.

I cannot provide more information about how the routines were programmed. They simply try to implement the partial sum, the stopping criteria being no significant differences between two consecutive terms of the series.

The problem could be related with the fact that the terms of the series are small at first, then grow fast and then they decrease again (at least for this concrete example), but I am not really sure.


As I tried to explain above, it is not a matter of adding more iterations or calculating more terms of the summation. The problem seems to be that some of the terms of the summation are somehow inaccurate, and from that moment, the sum gets away from the actual result.

My question is mainly about other ways to cope with this summation —maybe by transforming it into some equivalent expression, maybe applying any result for this kind of series...


I realised that As said in a comment, the previous series is an alternating one for $x<0$. It means that each term in the series is more or less cancelled by the following one. So, why not try to re-write it by adding up each two consecutive terms of the series? I did it, and this is the result: $$ \sum_{k=0}^{+\infty}{ \frac{ x^k }{\, a (a+1) \cdots (a+k) \,} } = \sum_{j=0}^{+\infty}{ \frac{ x^{2j}(a + 2j + x + 1) }{\, a (a+1) \cdots (a+2j+1) \,} } \, \mathrm{.} $$

Besides, I also learned about Kahan summation algorithm (Kahan, 1965), which aims at minimising rounding errors in summations.

I tried both approaches (both together and separately) but the results I get are still unsatisfactory.


Kahan W (1965). Further remarks on reducing truncation errors. Communications of the ACM, 8(1):40. doi:10.1145/363707.363723

  • 1
    $\begingroup$ For negative values of $x$, you get an alternating series, and thus,$$x<0\implies\left|\sum_{k=n}^\infty\frac{x^k}{a(a+1)\dots(a+k)}\right|<\frac{|x|^n}{a(a+1)\dots(a+n)}$$ $\endgroup$ – Simply Beautiful Art Oct 2 '17 at 14:25
  • 1
    $\begingroup$ Since you are using fairly large values of $x$, you should expect the error to be greater, and thus, it will require more terms. Also note that you or Mathematica may be running into floating error type stuff. $\endgroup$ – Simply Beautiful Art Oct 2 '17 at 14:28
  • 1
    $\begingroup$ What are the specific steps taken to calculate the terms? Do you calculate each term from scratch or do you use the current term to generate the next? It is entirely possible that you can do Richardson extrapolation here. If $S_n$ is the $n$ partial sum, then I expect you to have an expansion of the form $$S - S_n = \alpha n^{-p}+ O(n^{-q})$$. This may be hard to prove, but a numerical test can be done. Compute $$\frac{S_{2n}-S_n}{S_{4n}-S_{2n}}$$ for $n=1,2,4,8,...$. Error estimates can be derived if the fractions converge to $2^p$ for some $p$. $\endgroup$ – Carl Christian Oct 3 '17 at 18:11
  • 1
    $\begingroup$ Please don't waste any time on my previous suggestion. I have written a few words to the effect of why it will not work. $\endgroup$ – Carl Christian Oct 3 '17 at 20:56

This is not an answer, but too long for a comment.

For your particular case of $a=3/2$ and $x=-100/3$ the largest absolute value of the individual term is nearly $10^{11}$ before decay sets in. Regardless of how cleverly we sum the terms, i.e. Kahan summation or treewise summation, we carry with us an error which is at least $10^{11}u \approx 10^{-5}$ in double precision even if the floating point representation of the terms is obtained. Here $u$ is the unit roundoff. The error of $10^{-5}$ has the same order of magnitude as the error you report.

It is clear that you are pushing against the limit of what brute force can do and a clever rewrite is needed as you have already concluded.

It may be that


is a good place to start.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.