Today I learned about Spouge's formula to approximate the factorial.

Spouge's formula allows one to approximate the value of the gamma function. In case you don't know, the gamma function is like a generalisation of the factorial.

In fact, the following equality is true:

\[ \Gamma(z + 1) = z!\]

where \(\Gamma\) is the gamma function.

What Spouge's formula tells us is that

\[ \Gamma(z + 1) = (z + a)^{z + \frac12}e^{-z-a}\left( c_0 + \sum_{k=1}^{a-1} \frac{c_k}{z+k} + \epsilon_a(z) \right)\]

In the equality above, \(a\) is an arbitrary positive integer and \(\epsilon_a(z)\) is the error term. Thus, if we drop \(\epsilon_a(z)\), we get

\[ \Gamma(z + 1) = z! \approx (z + a)^{z + \frac12}e^{-z-a}\left( c_0 + \sum_{k=1}^{a-1} \frac{c_k}{z+k} \right)\]

The coefficients \(c_k\) are given by:

\[ \begin{cases} c_0 = \sqrt{2\pi} \\ c_k = \frac{(-1)^{k-1}}{(k - 1)!}(-k + a)^{k - \frac12}e^{-k+a}, ~ k \in \{1, 2, \cdots, a-1\} \end{cases}\]

By picking a suitable value of \(a\), one can approximate the value of \(z!\) up to a desired number of decimal places. Although we need the factorial function to compute the coefficients \(c_k\), those coefficients only need the factorial of numbers up to \(a - 2\). If we are approximating \(z!\), where \(a << z\), then this approximation saves us some work.

In order to determine the number of correct decimal places of the result, one needs to control the error term \(\epsilon_a(z)\). If \(a > 2\) and the \(Re(z) > 0\) (which is always true if \(z\) is a positive integer), then

\[ \epsilon_a(z) \leq a^{-\frac12}(2\pi)^{-a-\frac12}\]

By determining the value of \(a^{-\frac12}(2\pi)^{-a-\frac12}\), we can tell how many digits of the result will be correct. For example, with \(a = 10\), we get

\[ a^{-\frac12}(2\pi)^{-a-\frac12} \approx 1.31556 \times 10^{-9} ~ ,\]

meaning we will get 8 correct digits.

However, notice that the approximating formula must, itself, be computed with enough precision for the final result to hold as many correct digits as expected. In other words, if a higher value of \(a\) is picked so that the final result is more accurate, then we need to control the accuracy used when computing the coefficients \(c_k\) and the formula itself.

I'll leave it as an exercise for you, the reader, to implement this approximation in your favourite programming language.

In APL, (and disregarding the accuracy issues) it can look something like this:

```
β Computes the `c_k` coefficients:
Cks β {(.5*β¨β2),((!ks-1)Γ·β¨Β―1*ks-1)Γ((β΅-ks)*ks-.5)Γ*β΅-ksβ1+β³β΅-1}
β Computes the approximation of the gamma function:
GammaApprox β {((β΅+βΊ)*β΅+.5)Γ(*-β΅+βΊ)Γ(β’Γ·1,β΅+1ββ³ββ’)Cks βΊ}
β Computes an upper bound for the error term:
Err β {(β΅*Β―.5)Γ(β2)*-β΅+.5}
a β 10
Err a
1.315562187EΒ―9 β Thus, we expect 8 decimal places to be correct.
z β 100
a GammaApprox z
9.332621544E157
!z
9.332621544E157
```

That's it for now! Stay tuned and I'll see you around!

I write about Python every week. Join +16.000 others who are taking their Python π skills to the next level π, one email at a time.

- Spouge's approximation, Wikipedia, https://en.wikipedia.org/wiki/Spouge%27s_approximation [last accessed 10-01-2022];