Difficulty calculating the lower incomplete gamma function for very large and small values

Hi,

I've been using the expression gammainc(z,a)*gamma(a) to calculate the lower incomplete gamma function along with function vpa (to handle large numbers).

However, even when using vpa this fails at a certain point:

gamma(vpa(1e7))
 
ans =
 
1.2024234005159034561401534879443e+65657052
 
>> gamma(vpa(1e8))
 
ans =
 
Inf

I have tried taking logs, but have been unable to combine this with how I was calculating the lower incomplete gamma function (see previous question https://www.mathworks.com/matlabcentral/answers/540815-issue-with-gammainc-x-a-for-small-x-and-larger-a?s_tid=prof_contriblnk).

Is there a way to calculate both these large and small (very close to zero) values with necessary precision?

Thank you!

[Steven Lord]

Why are you trying to compute the gamma function for inputs on the order of tens or hundreds of millions, the results of which would have tens of millions (or more) digits? What's your use case? If you're computing ratios of factorials, for example, there may be ways to compute the ratios without having to compute the (extremely) large values separately.

[David Goodmanson]

Hi John,

in the expression gammainc(z,a)*gamma(a) are you interested in small x and large a, or large x and medium size a, or large x and large a with x>>a, or ---?

[John Fullerton]

Hi,

I need to evaluate the expression for all combinations of large and small 'a' and 'z' (Both varying roughly between ~1e-8 and ~1e8).

Thanks for your comments.

[David Goodmanson]

The function gammainc(z,a)*gamma(a) is simply

Integral{0,x} x^(a-1) e^-x dx.

For small x you can just expand e^-x in a Taylor series and integrate term-by-term to obtain

result = x^a ( 1/a - x/(a+1) + x^2/(2!(a+2)) - x^3/(3!(a+3)) ...  ) 

which is easy to calculate, although the factor in front takes special handling and in many cases can be expressed only in terms of its exponent.

Large x is harder, especially if 'a' is also large.

The range of 1e-8 to 1e8 for x is pretty reasonable, but since 'a' appears in the exponent in the integrand, large values of 'a' can lead to crazy large or small values for the result. For example if x = 1e-6 and a = 1e8 then the result is 10(-8) * 10^(-6*10^8). And if 'a' is 10^7 and x is larger than that, the result is (10^7)! which is approximately 10^(6.57*10^7). These numbers are so large (or small) that it is hard to see what their purpose might be. Could you explain the context for these values?

[John Fullerton]

Hi,

We're just trying to plot an equation that happens to contain the lower incomplete gamma function, and for the parameters we are interested in that's how a and x vary. It's based off previous work in Mathematica which I guess has some way of handling these crazy large and small numbers.

Is there no way of keeping precision for large values of a?

Thank you!