# Evaluation of the error term in Stirling's formula for factorial
**URL:** https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542
**Category:** Performance
**Created:** [March 6, 2018, 4:46pm UTC](https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542 "2018-03-06T16:46:06Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@dmbates](https://discourse.julialang.org/u/dmbates)
#### Post date: [March 6, 2018, 4:46pm UTC](https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542/1 "2018-03-06T16:46:06Z")
[Pull request #33](https://github.com/JuliaStats/StatsFuns.jl/pull/33/files) on the [`StatsFuns` package](https://github.com/JuliaStats/StatsFuns.jl) provides a link to Catherine Loader’s [`note`](https://lists.gnu.org/archive/html/octave-maintainers/2011-09/pdfK0uKOST642.pdf) on evaluating the binomial and Poisson probability mass function. One central part to this is evaluating the error in [Stirling’s approximation](https://en.wikipedia.org/wiki/Stirling%27s_approximation) to the log-factorial function. The version used is odd powers of `1/n` with coefficients that are inverses of integers. It starts as
\frac{1}{12n} -\frac{1}{320n^3}+\frac{1}{1260n^5}-\frac{1}{1680n^7}+\dots
Are there clever ways of evaluating such a polynomial so that, e.g., switching to BigInt would be transparent? The code used by Catherine Loader and translated more-or-less literally in the PR memoizes the values for 1 to 15 (where the error is greatest) then switches to Horner evaluation of polynomials with explicitly evaluated Float64 coefficients. I think it would be better to retain the integer representation of the inverse of the coefficients.
---
### Author:  [@dmbates](https://discourse.julialang.org/u/dmbates)
#### Post date: [March 6, 2018, 4:57pm UTC](https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542/2 "2018-03-06T16:57:27Z")
A bit more detail on the Stirling approximation. `log(factorial(n))` is approximated as `n * log(n) -n + log(2π*n)/2` For example
```julia
julia> n = big(35)
35
julia> log(factorial(n))
9.213617560368709248333303629689953216439487104597391356978356285232562502050595e+01
julia> n * log(n) - n + log(2π*n)/2
9.213379471607885835185644929557128163242134296609018476816483085946711233412839e+01
julia> n * log(n) - n + log(2π*n)/2 + inv(12n)
9.213617566845981073280883024795223401337372391847113714911721181184806471508052e+01
julia> n * log(n) - n + log(2π*n)/2 + inv(12n) - inv(320*n^3)
9.213617559557351335671553578731083459646410292721749283424840714712503264510946e+01
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 6, 2018, 5:03pm UTC](https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542/3 "2018-03-06T17:03:41Z")
Why not just use the `lgamma` function for the log of the factorial, since this already supports `BigFloat` with arbitrary precision?
You have to be careful with Stirling at arbitrary precision since it is an asymptotic series (i.e. the more precision you want, the bigger `n` needs to be).
---
### Author:  [@dmbates](https://discourse.julialang.org/u/dmbates)
#### Post date: [March 6, 2018, 5:07pm UTC](https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542/4 "2018-03-06T17:07:43Z")
As I was writing this I too was beginning to think that `lgamma` may be a better approach. Loader argues against it for the evaluation of the binomial pdf because of cancellation. I think it would be worthwhile checking out exactly how well it does relative to the simpler formula using `lgamma`.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 7, 2018, 2:28pm UTC](https://discourse.julialang.org/t/evaluation-of-the-error-term-in-stirlings-formula-for-factorial/9542/5 "2018-03-07T14:28:17Z")
> [@dmbates](#):
>
> Loader argues against it for the evaluation of the binomial pdf because of cancellation.
This is true. If we define `logp(x,n,p) = lgamma(n+1) - lgamma(x+1) - lgamma(n-x+1) + x*log(p) + (n-x)*log1p(-p)`, then `logp(1e6, 2e6, 0.5)` gives `-7.4801203440874815`, whereas the correct (exactly rounded) answer is `-7.480120346906837`, so in `Float64` precision you get 9 correct digits. (She claims that 7 digits are lost in this example, which is maybe too pessimistic by 1 digit.)
However, if you are using `BigFloat` anyway, then you can just throw more digits at it. This seems a lot easier than working out your own arbitrary-precision routine.
I should also mention that the combination of `lgamma` calls here is closely related to the `lbeta` function. Currently, `lbeta` in Base is computed via the combination of `lgamma` calls, but it might be nice to have a more sophisticated implementation (especially for `Float64`) that uses an asymptotic formula for large arguments.