The factorial function is rarely a good idea to use heavily. Large factorials tend to exceed the limits of double precision.
You generally have choices in such a problem. You can use loops, computing each term from the previous coefficient.
So in a Taylor series for exp(x), with terms of the form x^n/factorial(n), for each term you just mutiply the previous term by x, then divide by n. This takes advantage of the rule for factorial(n) as the product of the positive integers not exceeding n.
A ratio of factorials is simple, since most of those integers in a product just cancel out. So we have the simple identity for integers M<=N
factorial(N)/factorial(M) = prod(M+1:N)
Another way to deal with factorials is as you stated, to use logs. That is if your terms might look like this:
x^N * factorial(N)/factorial(M) = exp(N*log(X) + sum(log(1:N)) - sum(log(1:M))
Of course, again you could see that many of those terms will drop out of the sums.
Finally, use the gamma function, which is also available in log form. That is, remember that for integer N, we know
factorial(N) = gamma(N+1)
Often though, if we are working with logs, we can use the gammaln function, which computes the natural log of gamma(N), but it never computes the gamma function directly. And since factorial(N) will quickly overwhelm double precision computations, gammaln is a big improvement. That is, what is the ratio
factorial(2000)/factorial(1900)
ans =
NaN
We cannot directly compute this using factorials, even if we use logs.
log(factorial(2000)/factorial(1900))
ans =
NaN
log(factorial(2000)) - log(factorial(1900))
ans =
NaN
Nor can we compute it using other means, because the ratio itself is simply too large a number. I suppose we could use symbolic computations, thus
double(log(factorial(sym(2000))/factorial(sym(1900))))
ans =
757.57
But symbolic computations are highly computationally intensive. Things will get very slow.
However a simple solution in double precision arithmetic is to use gammaln.
gammaln(2000 + 1) - gammaln(1900 + 1)
ans =
757.57
So your problem is trivially solved, as we might start it like this:
q = 1:100;
N = 100;
u = 0.02;
x = sum(exp((q - 1)*log(u) + gammaln(N-1 + 1) - gammaln(N - q + 1)))
x =
3.0669e+09
