| 845 | |
| 846 | # Continued fraction approximation of the incomplete gamma function. |
| 847 | def gf(a, x, epsilon=3.e-7, iterations=200): |
| 848 | ln = gammaln(a) |
| 849 | g0 = 0.0 |
| 850 | a0 = 1.0 |
| 851 | b0 = 0.0 |
| 852 | a1 = x |
| 853 | b1 = 1.0 |
| 854 | f = 1.0 |
| 855 | for i in range(1, iterations): |
| 856 | a0 = (a1 + a0 * (i - a)) * f |
| 857 | b0 = (b1 + b0 * (i - a)) * f |
| 858 | a1 = x * a0 + a1 * i * f |
| 859 | b1 = x * b0 + b1 * i * f |
| 860 | if a1 != 0.0: |
| 861 | f = 1.0 / a1 |
| 862 | g = b1 * f |
| 863 | if abs((g - g0) / g) < epsilon: |
| 864 | return (g * exp(-x + a * log(x) - ln), ln) |
| 865 | g0 = g |
| 866 | raise StopIteration, (abs((g-g0) / g)) |
| 867 | |
| 868 | if a <= 0.0: |
| 869 | return 1.0 |