A fundamental, basic trick that anyone who is summing series should understand is WHY you benefit from the use of logs.
It is NOT that you can do exp(gammln(x)). That is a useless thing, becaue if gamma(x) would overflow, then so would the alternative form! Sorry, but it is just silly to think that would not happen. So how can you benefit? You can benefit from logs because the big terms in a product cancel out the small terms.
I'll just pick some random series, where I theoretically have no idea what the answer is, but I will expect that it will take multiple terms before convergence happens. So the nth term will be:
(x + 10)^n/gamma(n+1)
for n=0 to inf. (If you know what the sum is, SHHHSH! Keep it a secret!) So I expect that for small n, (x+10) will dominate things. But eventually gamma(n) will take over. Arbitrarily, I'll pick x = 1.
Ok, symsum will actually handle this.
syms n
x = 1;
symsum((x+10)^n/gamma(n+1),0,inf)
ans =
exp(11)
But suppose I did not have symsum? And suppose that I did not recognize that basic series? The terms of that series will have numerator and denominators that get really large.
n = (0:20)';
[(x+10).^n,gamma(n)]
ans =
1 Inf
11 1
121 1
1331 2
14641 6
161051 24
1771561 120
19487171 720
214358881 5040
2357947691 40320
25937424601 362880
285311670611 3628800
3138428376721 39916800
34522712143931 479001600
379749833583241 6227020800
4.17724816941565e+15 87178291200
4.59497298635722e+16 1307674368000
5.05447028499294e+17 20922789888000
5.55991731349223e+18 355687428096000
6.11590904484145e+19 6.402373705728e+15
6.7274999493256e+20 1.21645100408832e+17
But eventually, the gamma term will just take over, completely dominating the numerator. Depending on the value of x, conceivably, things might get so bad that we get overflows. Yet, I will assert the series itself converges for all real x. You still might need a high precision computation for large negative x, since there will be massive subtractive cancellation.
So, how might you sum this series? One obvious trick is to use gammln, the log of the gamma function. That is, compute the log of each term.
S = 0;
tol = eps;
term = inf;
n = 0;
while abs(term) > tol
term = exp(log(x+10)*n - gammaln(n+1));
S = S + term;
n = n + 1;
end
[n,S,exp(11)]
ans =
56 59874.1417151979 59874.1417151978
So we had to go all the way to n==56 to gain convergence. But there was no problem at all. In fact, for somewhat larger x, I'd bet I might see the parts of that term overflow, yet in theory, I could still have used the loop above and have gained convergence.
Now there are other trivial tricks I could have used. One is to simply compute one term directly from the one before. This can work as long as there is a simple way to get from term(n) to term(n+1). While you may still use logs here, you may be able to avoid it too. Then I might never even need to use gammaln. So in the following sum, I'll use the basic gamma identity
gamma(n+1) = gamma(n)*n
That allows me to write the sum as:
S = 0;
tol = eps;
term = 1;
n = 0;
while abs(term) > tol
S = S + term;
n = n + 1;
term = term*(x+10)/(n);
end
S = S + term;
[n,S,exp(11)]
ans =
55 59874.1417151978 59874.1417151978
See that I never needed to use logs here, or even gammaln, because each term was trivially computed from the term before.
There are other tricks for summing an infinite series. I'm actually a bit surprised you have not seen these old ones. If not, and you are going to do much of it, I'd suggest you do some serious reading. And, you also need to conside horw factorial(n) is related to the gamma function. Knowledge of some basic identities will help you, as well as understanding floating point arithmetic. So I would suggest you do some reading.
