Could anyone help/guide me towards properly analyzing the overflow as from the value of n = 10?

Question

The small program below calculates the absolute and relative errors in the Stirling's approximation $$n! = \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$$ For n = 1,2,...10.

Program StirlingTable implicit none

integer :: n, nf
real :: stirling, abs_error, pi, rel_error

pi = acos(-1.0) ! definition of the variable pi
nf = 1       ! starting value for n!

do n = 1, 10
    nf = nf*n
    stirling = sqrt(2.0*pi*n)*exp(-real(n))*((n)**n)
    abs_error = abs(nf - stirling)
    rel_error = abs_error/nf
    write(*,*) n, nf, stirling, &
                             abs_error, rel_error


end do
         read(*,*)

End Program StirlingTable

When compiling, the results give

Could anyone help me analyze the reasons behind why the absolute error increases and the relative error decreases with the increase in n, especially within the context of the code.

That is, which part of the code would be responsible for the overflow error.

Thank you

Answer 1

The problem is due to integer overflow in the computation of $n^n$. In the case of $n=10$ the computed result is $10^{10} \equiv 1410065408 \mod 2^{31}$, where $2^{31}$ corresponds to the largest positive number of type "integer" in Fortran. The problem is fixed by forcing the compiler to issue floating point instructions as demonstrated below

Program StirlingTable

  implicit none

  integer :: n, nf
  real    :: stirling, abs_error, pi, rel_error
  real    :: realn

  pi = acos(-1.0) ! definition of the variable pi
  nf = 1          ! starting value for n!

  do n = 1, 13
     nf = nf*n;
     realn=real(n);
     stirling = sqrt(2.0*pi*n)*exp(-real(n))*(realn**n)
     abs_error = abs(nf - stirling)
     rel_error = abs_error/nf
     write(*,"(I3,I12,3E14.7)") n, nf, stirling, abs_error, rel_error
  end do

End Program StirlingTable

A new problem emerges as the calculation of $n!$ overflows at $n=13$ and the computed value of nf equals $13! \equiv 1932053504 \mod 2^{31}$.

  1           1 0.9221370E+00 0.7786298E-01 0.7786298E-01
  2           2 0.1919004E+01 0.8099556E-01 0.4049778E-01
  3           6 0.5836209E+01 0.1637907E+00 0.2729845E-01
  4          24 0.2350618E+02 0.4938240E+00 0.2057600E-01
  5         120 0.1180192E+03 0.1980827E+01 0.1650689E-01
  6         720 0.7100782E+03 0.9921753E+01 0.1378021E-01
  7        5040 0.4980396E+04 0.5960400E+02 0.1182619E-01
  8       40320 0.3990239E+05 0.4176055E+03 0.1035728E-01
  9      362880 0.3595368E+06 0.3343156E+04 0.9212842E-02
 10     3628800 0.3598696E+07 0.3010400E+05 0.8295855E-02
 11    39916800 0.3961563E+08 0.3011720E+06 0.7544993E-02
 12   479001600 0.4756875E+09 0.3314080E+07 0.6918724E-02
 13  1932053504 0.6187240E+10 0.4255186E+10 0.2202416E+01

It could be a worthwhile addition to the exercise to also compute the logarithm of $n!$ and the logarithm of Stirling's approximation using the addition formula $\log(ab) = \log(a) + \log(b)$ as this would circumvent the current range issues.

Answer 2

Stirling approximation states that $$ n! \sim \sqrt{2\pi n} \left( \frac{n}{e} \right)^n, $$ meaning $$ \lim_{n \to \infty} \frac{n!}{\sqrt{2\pi n} \left( \frac{n}{e} \right)^n} = 1. $$

This does not mean that the following $$ \lim_{n \to \infty} \left( n! - \sqrt{2\pi n} \left( \frac{n}{e} \right)^n\right) = 0 $$ is true.

To feel how this can happen consider two simpler equivalent sequences, for instance $$ n^2 + n\sim n^2. $$ The absolute error is $$ (n^2 + n) - n^2 = n \to \infty, $$ but the relative error is $$ \frac{(n^2 + n) - n^2}{n^2} = \frac{1}{n} \to 0. $$

Answer 3

Let $$A=n! \qquad \text{and} \qquad B=\sqrt{2\pi n} \left( \frac{n}{e} \right)^n$$

Using Sirling approximation, we have $$\log \left(\frac{A}{B}\right) = \frac{1}{12 n}+O\left(\frac{1}{n^3}\right)\implies \frac{A}{B}=1+ \frac{1}{12 n}+O\left(\frac{1}{n^2}\right) $$ Then $$A-B\sim\frac{B}{12 n}=\frac 16 \sqrt{\frac{\pi }{2n}} \left( \frac{n}{e} \right)^n$$ increases very fast (just as $n \log(n)$) but $$\frac{A-B}{B}\sim \frac{1}{12 n}$$ decreases.