Stirling's Formula

> Compare n! to (n/e)^n * sqrt(2*Pi*n)

> [seq([100*n,evalf((100*n)!,10),evalf((100*n/exp(1))^(100*n)*sqrt(2*Pi*100*n),10)],n=1..10)];

> Compare ln(n!) to n*ln(n) -n+(1/2) * ln(2*Pi*n)

> S := n -> n*ln(n) -n+(1/2) * ln(2*Pi*n);

> [seq([100*n,evalf(ln((100*n)!),10),evalf(S(100*n),10)],n=1..10)];

> er := (n,k) -> sum(evalf(bernoulli(2*j)/((2*j-1)*(2*j)*n^(2*j-1)),50),j=1..k);

> [seq([10*k,evalf(er(10,10*k),50)],k=1..10)];

> The actual error is given by

> evalf(ln(10!)-S(10),50);

> This next program finds the smallest summand in the asymptotic series for the error function.

> summand := j -> abs(bernoulli(2*j)/((2*j-1)*(2*j)*10^(2*j-1)));

> k=1; while summand(k) > summand(k+1) do k:=k+1 end do; print(k);

> evalf(summand(32));

> This will be the greatest possible accuracy in the error function. We set the accuracy of our calculations to 50 digits to ensure that machine rounding does not affect accuracy of our results. We first compute the approximation of Log[100!] given by Stirling's formula.

> evalf(S(10),50);

> We next computee the error term given by the asymptotic series with 32 summands.

> sum(evalf(bernoulli(2*j)/((2*j-1)*(2*j)*10^(2*j-1)),50),j=1..32);

> We now compare the exponential of Stirling's formula, the exponential of Stirling's formula with the error term added in, and the actual value of 10!.

> round(evalf(exp(S(10)),50));

> round(evalf(exp(S(10)+sum(evalf(bernoulli(2*j)/((2*j-1)*(2*j)*10^(2*j-1)),50),j=1..32)),50));

> 10!;