Sadly, you are limited in that no matter how hard you try, you will not be able to compute it in double precision when it is that small. It will underflow. That leaves you with several choices.
- Use a tool that can compute in higher precision, that is either my HPF or syms.
- Use logs, and an approximation that will be valid for small x and large a.
If you cannot store the number in double precision, you really have no choice beyond those two cases. The number 1e-358 simply does not exist in the world of double precision. It is just zero. But let me give a few examples of things I might try.
Also, we need to worry about how MATLAB defines the incomplete gamma. That is, MATLAB normalizes by gamma(a). From the help for gammainc, we see:
The incomplete gamma function is defined as:
gammainc(x,a) = 1 ./ gamma(a) .*
integral from 0 to x of t^(a-1) exp(-t) dt
So, if you want to compute the same result as say Wolfram Alpha, which probably does not have that same normalization built in, we need to be careful.
Looking at the wiki site, I see a recurrence relation for the lower incomplete gamme. I'll call it g(a,x). Thus...
g(a+1,x) = a*g(a,x) - x^a*exp(-x)
How does this help? Say we want a 16 digit approximation to g(a,x). I only want 17 digits, because I'll use gammaln at the very end. I'll do it first using HPF. The trick is to run the recurrence BACKWARDS. My computations suggest we will need to take roughly 4 terms for 16 digits of precision, but I'll take back up the reverse recurrence by 10 terms to be safe.
a = 100;
x = hpf('0.01');
a0 = a + 10;
gax = hpf(0);
for ai = a0:-1:a
gax = (gax + x^ai*exp(-x))/ai;
gax = gax/exp(hpf(gammaln(100)));
Of course, I could have done the same using the symbolic toolbox. The trick is always to start out at zero, and then run backwards. This is a common trick for many problems. Personally, I've always loved the idea of running a recurrence backwards. It often works remarkably well, at least if you do so on the right problem.
Next, I could have used a simple series approximation. In fact, this should converge pretty well, as long as x is small. We would have:
For S as the sum from 0 to infinity of the terms
x^k/gamma(a + k + 1)
g(a,x) = x^a*gamma(a)*exp(-x)*S
And of course, if we use the MATLAB incomplete gamma form, where we divide by gamma(a), it gets even simpler. I'll do this one using syms and vpa, just for fun, and to show you how you might approach it with syms.
a = sym(100);
x = sym('0.01');
S = 1/gamma(sym(a)+1);
for k = 1:50
S = S + x^k/gamma(a + k + 1);
gax = S*x^a*exp(-x);
I should be able to do this another way, if I want to compute a result in a log form, so entirely in double precision.
Here it is easiest to use gammainc with the 'scaledlower' option, which takes all the nasty stuff out.
x = 0.01;
a = 100;
gaxln = log(gammainc(x,a,'scaledlower'));
gaxln = gaxln - gammaln(a + 1) - x + a*log(x);
Done entirely in double precision, the natural log of your result is:
You cannot exponentiate it of couse as a double, because it will underflow. But we can use vpa to verify the result.
Which looks reasonable. If you need to work in double precision, this may be your best choice, as long as the log is sufficient.
For a very simple approximation, we can use the asymptotic behavior of the lower incomplete gamma. Again, looking at the wiki page I cited, I see the limit as:
g(a,x)/x^a --> 1/a
as x --> 0.
Therefore, we can compute an approximation to the natural log of that lower incomplete gamma as
x = 0.01;
a = 100;
gaxlnapprox = a*log(x) - log(a) - gammaln(a)
which agrees to about 5 decimal digits with the value computed above. So not terribly bad as an approximation.
So you have some alternatives. With some thought, I'd bet I can come up with a few more.