If you are not careful, you will need to use a tool with literally thusands of digits of precision, as those terms will vary by thousands of digits in magnitude. Worse, if you even try to compute any of those terms, you will find that double precision either underflows or overflows.
vpa(factorial(sym(1247)))
8.458643820551953682887859707083e+3320
As such, in order to compute that series, you will need to use symbolic computations, or use a tool like my HPF. I'll do it using HPF, although I could have as easily used the symbolc toolbox too.
DefaultNumberOfDigits 100
0.8494300028975462166085239594512340095256328874564113849900813900015574782646329472378321829350836040
That is the lazy way to solve the problem. It involves no real effort, no thought about how to compute the result.
Can you solve the problem using double precision? Well yes. But again, you MUST be careful. And you must understand how to work with double precision in a way that does not underflow or overflow before you have gotten the result you need. The trick is often to use logs on these problems.
T = T + log(1211.5) - log(n+1);
In this sum, you should realize that almost every single term is effectively zero. So by using natural logs, I computed the log of every term, then exponentiate only after I have computed the log of the complete term.
As you should see, the two computations agree reasonably well. The double precision term is wrong in the last few digits, but you also need ot consider that every term in that series varies by as much as roughly 3 decimal digits. So that is probably as good as I can easily do using doubles.
Best Answer