We may need to calculate factorials \(n!\) from time to time. But we will easily get into a very large number as \(n\) increases. We all know Stirling’s formula for calculating factorials:
str[n_]:=n^n E^(-n) Sqrt[2 Pi n]
or \(n! \approx n^n e^{-n} \sqrt{2\pi n}\)
Error in Stirling’s formula
But how close is the Stirling’s formula? I did this in mathematica:
For[n=1, n<100, n=n+1, Print[n," ", N[n!/str[n]]]]
and the result is:
1 1.08444
2 1.04221
3 1.02806
4 1.02101
5 1.01678
6 1.01397
7 1.01197
8 1.01047
9 1.0093
10 1.00837
11 1.0076
12 1.00697
13 1.00643
14 1.00597
15 1.00557
16 1.00522
17 1.00491
18 1.00464
19 1.0044
20 1.00418
21 1.00398
22 1.00379
23 1.00363
24 1.00348
25 1.00334
26 1.00321
27 1.00309
28 1.00298
29 1.00288
30 1.00278
31 1.00269
32 1.00261
33 1.00253
34 1.00245
35 1.00238
36 1.00232
37 1.00225
38 1.0022
39 1.00214
40 1.00209
41 1.00203
42 1.00199
43 1.00194
44 1.0019
45 1.00185
46 1.00181
47 1.00177
48 1.00174
49 1.0017
50 1.00167
51 1.00164
52 1.0016
53 1.00157
54 1.00154
55 1.00152
56 1.00149
57 1.00146
58 1.00144
59 1.00141
60 1.00139
61 1.00137
62 1.00134
63 1.00132
64 1.0013
65 1.00128
66 1.00126
67 1.00124
68 1.00123
69 1.00121
70 1.00119
71 1.00117
72 1.00116
73 1.00114
74 1.00113
75 1.00111
76 1.0011
77 1.00108
78 1.00107
79 1.00106
80 1.00104
81 1.00103
82 1.00102
83 1.001
84 1.00099
85 1.00098
86 1.00097
87 1.00096
88 1.00095
89 1.00094
90 1.00093
91 1.00092
92 1.00091
93 1.0009
94 1.00089
95 1.00088
96 1.00087
97 1.00086
98 1.00085
99 1.00084
So you can see,
When | n!/str[n]-1 |
---|---|
n >= 9 |
\(< 1\%\) |
n >= 17 |
\(< 0.5\%\) |
n >= 28 |
\(< 0.3\%\) |
n >= 42 |
\(< 0.2\%\) |
n >= 84 |
\(< 0.1\%\) |
C-implementation
But actually, we can do more precise. According to MathWorld, we can have the inequality form for Stirling’s approximation:
\[n^n e^{-n} \sqrt{2\pi n} e^{1/(12n+1)} < n! < n^n e^{-n} \sqrt{2\pi n} e^{1/12n}\]also, we have the Gosper’s inequality which better approximate: \(n! \approx n^n e^{-n} \sqrt{\frac{2n+1}{3}\pi}\)
To compare, we can write a program as the one below:
/* Factorial approximations
* This program runs for n=1 to 1550, calculating approximations of n!
* Approximations used here are:
* 1. Stirling's formula
* 2. Lower bound of Stirling's inequality
* 3. Upper bound of Stirling's inequality
* 4. Gosper's approximation
* This program will show that, long double can only hold up to 1546!
*/
#include <stdio.h>
#include <math.h>
# define M_El 2.7182818284590452353602874713526625L /* e */
# define M_PIl 3.1415926535897932384626433832795029L /* pi */
long double factorial(int n) {
// Real factorial
long double val=1;
for(;n>1;n--) {
val *= n;
};
return val;
};
long double stirling(int n) {
// Stirling's approximation
return powl(n,n)*expl(-n)*sqrtl(2*M_PIl*n);
};
long double stirling_l(int n) {
// Stirling's lower bound
return powl(n,n)*expl(-n+1/(long double)(12*n+1))*sqrtl(2*M_PIl*n);
};
long double stirling_u(int n) {
// Stirling's lower bound
return powl(n,n)*expl(-n+1/(long double)(12*n))*sqrtl(2*M_PIl*n);
};
long double gosper(int n) {
// Gosper's approximation
return powl(n,n)*expl(-n)*sqrtl(M_PIl*(2*n+1/(long double)3));
};
int main(int argc, char**argv)
{
// Print a table showing the derivation between different factorial approximation
int n;
printf("%5s %15s %15s %15s %15s %15s\n",
"n","n!","Stirling","Stirling LB","Stirling UB","Gosper");
for(n=1;n<1550;n++) {
printf("%5d %15Lg %15Lg %15Lg %15Lg %15Lg\n",n,factorial(n),
stirling(n)/factorial(n)-1,
stirling_l(n)/factorial(n)-1,
stirling_u(n)/factorial(n)-1,
gosper(n)/factorial(n)-1);
};
return 0;
};
This program generates different approximations and compare it to the correct value of factorial. We can also shown by this program that, the long double precision in C can only hold up to factorial of 1546.
Below are the output of the C program up to factorial of 100. The number is the relative difference between approximation and actual value, i.e. approx/actual\(-1\):
n Stirling Stirling LB Stirling UB Gosper
1 -0.077863 -0.00412984 0.00227445 -0.00397819
2 -0.0404978 -0.0013398 0.000326024 -0.00131847
3 -0.0272984 -0.000650687 9.98571e-05 -0.000644108
4 -0.020576 -0.000382436 4.26621e-05 -0.000379603
5 -0.0165069 -0.000251217 2.19757e-05 -0.000249749
6 -0.0137803 -0.000177483 1.27601e-05 -0.000176626
7 -0.0118262 -0.000131995 8.05199e-06 -0.000131453
8 -0.0103573 -0.000101982 5.40142e-06 -0.000101616
9 -0.00921276 -8.1147e-05 3.79708e-06 -8.08896e-05
10 -0.00829596 -6.60984e-05 2.7699e-06 -6.59103e-05
11 -0.00754507 -5.4877e-05 2.08209e-06 -5.47353e-05
12 -0.00691879 -4.62873e-05 1.60434e-06 -4.6178e-05
13 -0.0063885 -3.95667e-05 1.26222e-06 -3.94805e-05
14 -0.0059337 -3.42098e-05 1.01084e-06 -3.41407e-05
15 -0.00553933 -2.98712e-05 8.22004e-07 -2.9815e-05
16 -0.00519412 -2.63084e-05 6.77414e-07 -2.62621e-05
17 -0.0048894 -2.33469e-05 5.64836e-07 -2.33082e-05
18 -0.00461846 -2.08586e-05 4.7588e-07 -2.0826e-05
19 -0.00437596 -1.87478e-05 4.04663e-07 -1.87201e-05
20 -0.00415765 -1.6942e-05 3.46975e-07 -1.69181e-05
21 -0.00396009 -1.53849e-05 2.9975e-07 -1.53643e-05
22 -0.00378045 -1.40331e-05 2.60719e-07 -1.40152e-05
23 -0.00361641 -1.28518e-05 2.28181e-07 -1.28362e-05
24 -0.003466 -1.18137e-05 2.00839e-07 -1.17999e-05
25 -0.00332761 -1.08964e-05 1.77697e-07 -1.08842e-05
26 -0.00319984 -1.0082e-05 1.57977e-07 -1.00711e-05
27 -0.00308152 -9.35556e-06 1.4107e-07 -9.34585e-06
28 -0.00297164 -8.7049e-06 1.26493e-07 -8.69619e-06
29 -0.00286932 -8.11982e-06 1.13856e-07 -8.11197e-06
30 -0.00277382 -7.5918e-06 1.02848e-07 -7.58471e-06
31 -0.00268447 -7.11366e-06 9.32145e-08 -7.10723e-06
32 -0.00260069 -6.6793e-06 8.47474e-08 -6.67346e-06
33 -0.00252199 -6.28354e-06 7.72755e-08 -6.27822e-06
34 -0.00244791 -5.92194e-06 7.06567e-08 -5.91707e-06
35 -0.00237806 -5.59068e-06 6.47727e-08 -5.58621e-06
36 -0.00231208 -5.28645e-06 5.95243e-08 -5.28235e-06
37 -0.00224966 -5.0064e-06 5.48279e-08 -5.00262e-06
38 -0.00219053 -4.74802e-06 5.06129e-08 -4.74453e-06
39 -0.00213442 -4.50915e-06 4.6819e-08 -4.50591e-06
40 -0.00208112 -4.28785e-06 4.3395e-08 -4.28485e-06
41 -0.00203042 -4.08245e-06 4.02969e-08 -4.07967e-06
42 -0.00198212 -3.89147e-06 3.74869e-08 -3.88888e-06
43 -0.00193607 -3.71358e-06 3.49321e-08 -3.71117e-06
44 -0.00189211 -3.54762e-06 3.26043e-08 -3.54536e-06
45 -0.00185011 -3.39253e-06 3.04789e-08 -3.39043e-06
46 -0.00180993 -3.2474e-06 2.85342e-08 -3.24543e-06
47 -0.00177145 -3.11138e-06 2.67515e-08 -3.10954e-06
48 -0.00173458 -2.98374e-06 2.51142e-08 -2.982e-06
49 -0.00169921 -2.86379e-06 2.36079e-08 -2.86216e-06
50 -0.00166526 -2.75093e-06 2.22197e-08 -2.7494e-06
51 -0.00163263 -2.64462e-06 2.09382e-08 -2.64317e-06
52 -0.00160126 -2.54435e-06 1.97534e-08 -2.54298e-06
53 -0.00157107 -2.44967e-06 1.86563e-08 -2.44838e-06
54 -0.001542 -2.36019e-06 1.7639e-08 -2.35897e-06
55 -0.00151399 -2.27551e-06 1.66943e-08 -2.27436e-06
56 -0.00148697 -2.19532e-06 1.58159e-08 -2.19423e-06
57 -0.00146091 -2.11929e-06 1.4998e-08 -2.11825e-06
58 -0.00143574 -2.04714e-06 1.42356e-08 -2.04616e-06
59 -0.00141142 -1.97862e-06 1.3524e-08 -1.97768e-06
60 -0.00138791 -1.91348e-06 1.28591e-08 -1.91259e-06
61 -0.00136518 -1.8515e-06 1.2237e-08 -1.85065e-06
62 -0.00134317 -1.79249e-06 1.16544e-08 -1.79168e-06
63 -0.00132187 -1.73625e-06 1.11082e-08 -1.73548e-06
64 -0.00130123 -1.68262e-06 1.05956e-08 -1.68189e-06
65 -0.00128122 -1.63144e-06 1.01141e-08 -1.63074e-06
66 -0.00126182 -1.58255e-06 9.66134e-09 -1.58188e-06
67 -0.001243 -1.53583e-06 9.23518e-09 -1.53519e-06
68 -0.00122473 -1.49115e-06 8.83373e-09 -1.49054e-06
69 -0.00120699 -1.44839e-06 8.45521e-09 -1.44781e-06
70 -0.00118976 -1.40745e-06 8.09801e-09 -1.40689e-06
71 -0.00117301 -1.36822e-06 7.76065e-09 -1.36768e-06
72 -0.00115673 -1.3306e-06 7.44177e-09 -1.33009e-06
73 -0.00114089 -1.29452e-06 7.14012e-09 -1.29402e-06
74 -0.00112549 -1.25988e-06 6.85456e-09 -1.2594e-06
75 -0.00111049 -1.22661e-06 6.58403e-09 -1.22616e-06
76 -0.00109588 -1.19465e-06 6.32754e-09 -1.19421e-06
77 -0.00108166 -1.16392e-06 6.08421e-09 -1.1635e-06
78 -0.0010678 -1.13436e-06 5.8532e-09 -1.13395e-06
79 -0.00105429 -1.10591e-06 5.63373e-09 -1.10552e-06
80 -0.00104112 -1.07851e-06 5.42511e-09 -1.07814e-06
81 -0.00102827 -1.05213e-06 5.22665e-09 -1.05177e-06
82 -0.00101574 -1.0267e-06 5.03776e-09 -1.02635e-06
83 -0.00100351 -1.00218e-06 4.85786e-09 -1.00184e-06
84 -0.000991567 -9.78528e-07 4.68643e-09 -9.78203e-07
85 -0.000979907 -9.55704e-07 4.52297e-09 -9.55391e-07
86 -0.000968519 -9.3367e-07 4.36702e-09 -9.33367e-07
87 -0.000957392 -9.12389e-07 4.21816e-09 -9.12096e-07
88 -0.000946517 -8.91827e-07 4.07599e-09 -8.91545e-07
89 -0.000935887 -8.71952e-07 3.94014e-09 -8.7168e-07
90 -0.000925494 -8.52735e-07 3.81026e-09 -8.52471e-07
91 -0.000915328 -8.34146e-07 3.68603e-09 -8.33891e-07
92 -0.000905383 -8.16158e-07 3.56713e-09 -8.15911e-07
93 -0.000895653 -7.98746e-07 3.4533e-09 -7.98507e-07
94 -0.000886129 -7.81886e-07 3.34426e-09 -7.81654e-07
95 -0.000876805 -7.65553e-07 3.23976e-09 -7.65329e-07
96 -0.000867676 -7.49727e-07 3.13957e-09 -7.4951e-07
97 -0.000858735 -7.34387e-07 3.04347e-09 -7.34176e-07
98 -0.000849976 -7.19512e-07 2.95125e-09 -7.19308e-07
99 -0.000841394 -7.05086e-07 2.86272e-09 -7.04887e-07
100 -0.000832983 -6.91088e-07 2.7777e-09 -6.90896e-07