Factorials grow very quickly, so quickly the results can easily exceed the limits of computer arithmetic. For example, 30! is too big to store in a typical integer and 200! is too big to store (even approximately) in a standard floating point number. It’s usually more convenient to work with the logarithms of factorials rather than factorials themselves.

So how might you compute log( *n*! )? You don’t want to compute n! first and then take logs because you’ll overflow for moderately large arguments. You want to compute the log factorial directly.

Since *n*! = 1 × 2 × 3 × … × *n*, log(*n*!) = log(1) + log(2) + log(3) + … + log(*n*). That gives a way to compute log(*n*!), but it’s slow for large arguments: the run time is proportional to the size of *n*.

If you only need to compute log(*n*!) for n within a moderate range, you could just tabulate the values. Calculate log(*n*!) for *n* = 1, 2, 3, …, *N* by any means, no matter how slow, and save the results in an array. Then at runtime, just look up the result.

But suppose you want to be able to compute log(*n*!) for any value of *n* such that log(*n*!) won’t overflow. That requires a very large value of *n*. Since log(*n*!) is on the order of *n* log (*n*) for large *n*, log(*n*!) won’t overflow even for some very large values of *n*. You’ll run out of memory to store your table before you’ll run out of values of *n* such that log(*n*!) doesn’t overflow.

So now the problem becomes how to evaluate log(*n*!) for large values of *n*. Say we tabulate log(*n*!) for *n* up to some size and then use a formula to calculate log(*n*!) for larger arguments. I’m going to switch notation now and work with the gamma function Γ(*x*) because most references state their results in terms of the gamma function rather than in terms of factorial. It’s easy to move between the two since Γ(*n*+1) = *n*!.

Stirling’s approximation says that

log Γ(*x*) ≈ (*x* – 1/2) log(*x*) – *x* + (1/2) log(2 π)

with an error on the order of 1/12x. So if *n* were around 1000, the approximation would be good to about four decimal places. What if we wanted more accuracy but not a bigger table? Stirling’s approximation above is part of an infinite series of approximations, an asymptotic series for log Γ(x):

log Γ(*x*) ≈ (*x* – 1/2) log(*x*) – *x* + (1/2) log(2 π) + 1/(12 *x*) – 1/(360 *x*^{3}) + 1/(1260 *x*^{5}) – …

This raises a couple questions.

- What is the form of a general term in the series?
- How accurate is the series when we truncate it at some point?

The answer to the first question is the term with *x*^{2m-1} in the denominator has coefficient *B*_{2m} / (2*m*(2*m*-1)) where the *B*‘s are Bernoulli numbers. Perhaps that’s not very satisfying, but that’s what it is.

Now to the question of accuracy. If you’ve never used asymptotic series before, you might be tempted to use one like you’d use a Taylor series: the more terms the better. But asymptotic series don’t work that way. Typically the error improves at first as you take more terms, reaches a minimum, then *increases* as you take more terms. Another difference is that while Taylor series approximations improve as arguments get smaller, asymptotic approximations improve as arguments get larger. That’s convenient for us since we’re looking for an approximation for n so large that it’s beyond our table of saved values.

For this particular series, the absolute value of the error in truncating the series is less than the absolute value of the first term that was left out. Suppose we make a look-up table for the values 1 through 256. If we truncate the series after 1/(12 *x*), the error will be less than 1/(360 *x*^{3}). If *x* > 256, log(*x*!) > 1419 and the error in the asymptotic approximation is less than 1/(360×2^{24}) = 1.65 × 10^{-10}. Since the number we’re computing has at least four digits and the result is good to 10 decimal places, we should have at least 14 significant figures, near the limits of floating point precision. (For more details, see Anatomy of a floating point number.)

In summary, one way to compute log factorial is to pre-compute log(*n*!) for *n* = 1, 2, 3, … 256 and store the results in an array. For values of *n* ≤ 256, look up the result from the table. For *n* > 256, return

(*x* – 1/2) log(*x*) – *x* + (1/2) log(2 π) + 1/(12 *x*)

with *x* = *n* + 1. This has been coded up here.

You could include the 1/(360 *x*^{3}) term or higher terms from the asymptotic series and use a smaller table. This would use less memory but would require more computation for arguments outside the range of the table.

The glibc implementation of the lgamma function (math.h) is similar, although it appears to use some other tricks to get good accuracy for small arguments.

Implementing log gamma is more complicated. I just considered factorial (i.e. integer arguments) here to simplify things a bit.

This post was motivated by this code for generating Poisson random values which only requires log factorial and not the full log gamma.