Author: David Max MA DPhil (2017-11-08)
A Fortran-77 program to test Stirling's formula for computing factorials
Factorials n! can be calculated directly for small n. For large n, however, the computation becomes cumbersome.
Stirling's formula is a commonly-used approximation for computing (the log of) factorial n for large n.
Program STIRLING is designed to test the accuracy of Stirling's approximation for a few values of n.
There are two forms of the approximation, the first of which is simply nln(n)-n. Program STIRLING computes the second, more sophisticated form.
The program is written in Fortran-77 and compiles okay using the gcc compiler (i.e., gfortran).
The computations are done using 8-byte precision (REAL*8).
STIRLING Fortran program listing
Fortran program for download
The Fortran-77 source file, together with a batch file to compile it and a file containing output results, are available to download in either tar or zip archive form.
On clicking on one of the links, your browser should ask you what you want to do with the selected archive file.
First, save the file into a working directory (or folder, in Windows-speak) of your choice.
Under Linux/Unix, the tar archive file can be unpacked using the shell command...
tar -xf STIRLING.tar.
(tar is also available for Windows).
On Windows the zip archive file can be unpacked by right clicking on the file and choosing "Extract All...".
If gcc has been installed, the program file STIRLING.exe can be built on Windows by running the batch file MAKE.STIRLING2.bat from the command line.
The batch file provided may be a useful starting place for compiling the source with an alternative compiler (Intel, Salford, etc).
Results for n up to 50
The table below gives some results from program STIRLING.
- Column 1 contains relevant value of n.
- Column 2 contains the approximate values for ln(n!) using the more sophisticated form of the approximation.
- Column 3 is the estimate for ln(n!) obtained by recursive calculation starting from n = 1.
- Column 4 is the approximate fractional error in the Stirling formula estimate for ln(n!). This gets smaller as n increases.
- Column 5 shows, for completeness, the value of ln(n!) obtained from Abramovich and Stegun [1] (p.274) (figure there is multiplied by ln(10)).
| n | STIRLING | ln(n!) recursive | frac. error | A. & S. |
|---|---|---|---|---|
| 1 | -0.081061 | 0.000000 | 0.000000 | 0.000000 |
| 2 | 0.651806 | 0.693147 | -0.059642 | 0.693147 |
| 3 | 1.764082 | 1.791759 | -0.015447 | 1.791759 |
| 4 | 3.157263 | 3.178054 | -0.006542 | 3.178054 |
| 5 | 4.770847 | 4.787492 | -0.003477 | 4.787492 |
| 6 | 6.565375 | 6.579251 | -0.002109 | 6.579251 |
| 7 | 8.513265 | 8.525161 | -0.001395 | 8.525161 |
| 8 | 10.594192 | 10.604603 | -0.000982 | 10.604603 |
| 9 | 12.792572 | 12.801827 | -0.000723 | 12.799525 |
| 10 | 15.096082 | 15.104413 | -0.000552 | 15.104412 |
| 11 | 17.494734 | 17.502308 | -0.000433 | 17.502308 |
| 12 | 19.980272 | 19.987214 | -0.000347 | 19.987215 |
| 13 | 22.545755 | 22.552164 | -0.000284 | 22.552164 |
| 14 | 25.185270 | 25.191221 | -0.000236 | 25.191220 |
| 15 | 27.893717 | 27.899271 | -0.000199 | 27.899272 |
| 16 | 30.666652 | 30.671860 | -0.000170 | 30.671861 |
| 17 | 33.500172 | 33.505073 | -0.000146 | 33.505075 |
| 18 | 36.390816 | 36.395445 | -0.000127 | 36.395445 |
| 19 | 39.335499 | 39.339884 | -0.000111 | 39.339885 |
| 20 | 42.331450 | 42.335616 | -0.000098 | 42.335617 |
| 21 | 45.376171 | 45.380139 | -0.000087 | 45.380139 |
| 22 | 48.467394 | 48.471181 | -0.000078 | 48.471182 |
| 23 | 51.603053 | 51.606676 | -0.000070 | 51.606675 |
| 24 | 54.781257 | 54.784729 | -0.000063 | 54.784730 |
| 25 | 58.000272 | 58.003605 | -0.000057 | 58.003606 |
| 26 | 61.258497 | 61.261702 | -0.000052 | 61.261702 |
| 27 | 64.554452 | 64.557539 | -0.000048 | 64.557539 |
| 28 | 67.886767 | 67.889743 | -0.000044 | 67.889744 |
| 29 | 71.254166 | 71.257039 | -0.000040 | 71.257039 |
| 30 | 74.655459 | 74.658236 | -0.000037 | 74.658236 |
| 31 | 78.089535 | 78.092224 | -0.000034 | 78.092224 |
| 32 | 81.555355 | 81.557959 | -0.000032 | 81.557960 |
| 33 | 85.051942 | 85.054467 | -0.000030 | 85.054468 |
| 34 | 88.578377 | 88.580828 | -0.000028 | 88.580828 |
| 35 | 92.133795 | 92.136176 | -0.000026 | 92.136176 |
| 36 | 95.717380 | 95.719695 | -0.000024 | 95.719694 |
| 37 | 99.328360 | 99.330612 | -0.000023 | 99.330613 |
| 38 | 102.966006 | 102.968199 | -0.000021 | 102.968198 |
| 39 | 106.629624 | 106.631760 | -0.000020 | 106.631760 |
| 40 | 110.318556 | 110.320640 | -0.000019 | 110.320640 |
| 41 | 114.032179 | 114.034212 | -0.000018 | 114.034212 |
| 42 | 117.769897 | 117.771881 | -0.000017 | 117.771881 |
| 43 | 121.531144 | 121.533082 | -0.000016 | 121.533082 |
| 44 | 125.315377 | 125.317271 | -0.000015 | 125.317270 |
| 45 | 129.122082 | 129.123934 | -0.000014 | 129.123934 |
| 46 | 132.950763 | 132.952575 | -0.000014 | 132.952576 |
| 47 | 136.800950 | 136.802723 | -0.000013 | 136.802724 |
| 48 | 140.672188 | 140.673924 | -0.000012 | 140.673924 |
| 49 | 144.564043 | 144.565744 | -0.000012 | 144.565744 |
| 50 | 148.476100 | 148.477767 | -0.000011 | 148.477767 |
The results for the 'recursive' column agree with the Abramovich and Stegun table.
The small differences might be the result of using lower (8-byte) precision here than in the computations for the A. & S. tables.
The form of the Stirling formula computed here agrees with the accurate value to 5.5 parts in 104 with n = 10.
At n = 50, agreement has reached 1.1 parts in 105.
For comparison, the very simple form of the approximation, nln(n)-n, returns 13.036 (for n = 10) and 145.601 (n = 50).
References
[1] Abramowitz, M. & Stegun, I.A. 1965. Handbook of Mathematical Functions. Dover, New York.