Incomplete gamma function continued fraction converges initially to wrong value

Question

I'm trying to develop an algorithm to calculate the incomplete gamma function for complex input values to arbitrary precision.

For some ranges of input values, I use the continued fraction formula:

$$\Gamma(a, z)= \cfrac{z^a e^{-z}}{1+z-a+ \cfrac{a-1}{3+z-a+ \cfrac{2(a-2)}{5+z-a+ \cfrac{3(a-3)} {7+z-a+ \cfrac{4(a-4)}{9+z-a+ \ddots}}}}}$$

I'm calculating the continued fraction with the (modified) Lentz's algorithm.

For a wide range of input values this algorithm works fine (if $a$ and $z$ aren't too large in absolute value and $z$ is not too close to the negative real axis and $|a|$ is roughly the same or less than $|z|$).

However there are some relatively narrow ranges of $a$ and $z$ where the continued fraction does not converge well. It initially converges to a wrong value - possibly to a precision of 100 digits or more - until it diverges away from the wrong value and then finally converges to the correct (?) value.

For example:

$\Gamma(-1 \cdot 10^{-2} + 4 \cdot 10^2i, 1 \cdot 10^{-2} + 1 \cdot 10^2i)$

When I run the algorithm with my arbitrary precision software using 150 digits of precision the continued fraction (excluding the factor of $z^a e^{-z}$) converges as follows.

In the first 240 iterations it converges to approximately this value (to a precision of 101 digits): $-3.4812455650395559724083225442044685679172703381183511740233299928770701153447585747686976555814151526 \cdot 10^{-6} + 3.3333097415343713576386150215385070073982379130010751824025858800036447026117907322447199186033463044446\cdot 10^{-3}i$

It circles around that point for maybe 120 iterations more and then starts to diverge away from it.

At 1240 iterations (from the beginning) it has diverged fully away and starts to move to a completely different point. The value changes substantially.

At around 1300 iterations it starts to converge to a different point (which I guess is the correct one):

$-8.88 \cdot 10^{-2} + 8.39 \cdot 10^{-2}i$

It keeps converging and after 2647 iterations the algorithm stops as it thinks it has the value correctly to a precision of 150 digits.

Assuming that Mathematica returns the correct value, it is:

Gamma[-1*10^-2 + 4*10^2 I, 1*10^-2 + 1*10^2 I] // N
-1.05428*10^-274 - 1.20736*10^-274 I

I can also verify the convergence of the continued fraction with Mathematica, using exact values, so I'm fairly certain that the problem isn't just error propagation in the Lentz's algorithm:

z^a Exp[-z] ContinuedFractionK[If[n == 1, 1, (n - 1) (a - n + 1)], (2 n - 1 + z - a), {n, 1, 240}] /. a -> -1*10^-2 + 4*10^2 I /. z -> 1*10^-2 + 1*10^2 I // N
-4.36799*10^-276 - 1.7607*10^-277 I

z^a Exp[-z] ContinuedFractionK[If[n == 1, 1, (n - 1) (a - n + 1)], (2 n - 1 + z - a), {n, 1, 1240}] /. a -> -1*10^-2 + 4*10^2 I /. z -> 1*10^-2 + 1*10^2 I // N
-4.36917*10^-276 - 1.75839*10^-277 I

z^a Exp[-z] ContinuedFractionK[If[n == 1, 1, (n - 1) (a - n + 1)], (2 n - 1 + z - a), {n, 1, 1300}] /. a -> -1*10^-2 + 4*10^2 I /. z -> 1*10^-2 + 1*10^2 I // N
-1.05267*10^-274 - 1.21171*10^-274 I

z^a Exp[-z] ContinuedFractionK[If[n == 1, 1, (n - 1) (a - n + 1)], (2 n - 1 + z - a), {n, 1, 2647}] /. a -> -1*10^-2 + 4*10^2 I /. z -> 1*10^-2 + 1*10^2 I // N
-1.05428*10^-274 - 1.20736*10^-274 I

Problem: If I want to calculate the result e.g. to just 50 digits of precision then there's really no way to know that it initially converges to a wrong value (because it converges to the wrong value to a precision of 101 digits).

Questions:

• How can I know, when the algorithm initially converges to an incorrect value and will only converge to the correct value later with more iterations and higher precision?
• Is there any way to estimate what is the minimum precision needed to detect that it initially converges to an incorrect value?
• Is there a better algorithm (or set of algorithms) for the complex incomplete gamma function that converges for all possible input values (for $\gamma(a,z)$ I already have a different continued fraction and the sum formulae for small $|z|$ and for negative integer $a$ yet a different algorithm)?

Answer 1

To answer my own question, this is the ad hoc algorithm that I came up with.

The basic idea is to observe if the continued fraction converges (close enough to the initial working precision, e.g. 50 digits) and then diverges (close enough to zero precision). Both of these can be easily seen from the delta in Lentz's algorithm; how close delta is to one.

Empirical evidence suggests that this approach works even if the initial working precision is far less than what is the precision where the continued fraction starts to diverge after the initial convergence. I have no mathematical proof of this whatsoever.

The problem with this approach is still that you don't know how many iterations of the continued fraction should be run, to see if it will diverge or not, after the initial convergence. Heuristically I'm using the point where the $|a_n|$ or $|b_n|$ start increasing again (whichever $n$ is larger), after initially decreasing. However in practice this seems to often give a much larger value than what is actually needed to observe the situation. Again I have no mathematical proof of this so I have no idea if it is actually true or not.

Thus, depending on if the continued fraction for the upper gamma function or lower gamma function is used, we can compute a minimum number of iterations that should indicate if the algorithm ultimately will converge or not.

In pseudo-code the algorithm is then:

• Compute the minimum number of iterations
• Observe the maximum precision that the algorithm has reached so far
• If the current precision is significantly lower than the maximum precision that has been reached so far, and the maximum precision is close enough to the working precision, abort the algorithm, increase the working precision (e.g. double it) and start from the beginning
• Repeat until the minimum number of iterations has been reached and the desired precision has been reached

This way the working precision is increased, until it is sufficient to hold the needed significant digits for the initial mis-convergence, divergence and then actual convergence to the correct value.

It should be noted that this approach does not actually work for the input values in the above example, because the minimum number of iterations there is zero. How this is solved is that if $|a|$ and $|z|$ are relatively close to each other (within a factor of 100 for example) then it's not certain, whether the upper gamma continued fraction or the lower gamma continued fraction works better. This can be tested by running both algorithms with a low fixed precision (e.g. 50 digits) to a maximum of e.g. 50 iterations and by comparing, which one converges faster. In the above example, the lower gamma continued fraction is chosen, and it does not seem to exhibit this convergence problem with the input values in the example. Once more I have no proof whatsoever that this approach would work in all cases.

With some other input values the approach works, e.g. with $\Gamma(600 + 1000i, -8 + 1000i)$

• The minimum number of iterations is 602
• On the first attempt, with working precision 56 it reaches precision 56 and then goes down to precision 9 when n=624 and the algorithm restarts with a higher working precision
• On the second attempt, with working precision 96 it reaches precision 96 and then goes down to precision 9 when n=764 and the algorithm restarts with a higher working precision
• On the third attempt, with working precision 176 it reaches precision 125 and then goes down to precision 0 and then reaches the target precision and this is accepted as the result should be now accurate, n=1360

Interestingly, even Mathematica does not always give consistent answers, for example I get with Mathematica 11.0:

In[1]:= N[Gamma[600 + 1000 I, -8 + 1000 I]]

Out[1]= -1.081029909234896*10^1115 - 1.16223826754520*10^1114 I

In[2]:= N[Gamma[600 + 1000 I, -8 + 1000 I], 150]

Out[2]= 7.\
1855853396209541562579560379061894335094757512323456565316919197701502\
8877260926174769914003649112663217993409550429804705379135275340692346\
731089907*10^1130 + 
 1.3260769500785825064134190817894642899524344648528512709130316162950\
2184076799914763130652533608907210383212359711735409448320421787750355\
051891174834*10^1130 I