A technique for precise computation with factorials in a digital computer
Fernando Jose Corbató · 1961
Computation involving factorials in a digital computer has various hazards. Foremost of the hazards is the rapid build-up of magnitude with argument n, for 10: @@@@ 106 and 20: @@@@ 1018. Clearly a fixed-point representation will not suffice to avoid overflow except for trivial values of n; a floating-point representation will offer a reprieve for the range overflow but will contain only a limited number of significant figures. For most purposes computing indirectly with a “compression-function” such as the logarithm will eliminate range overflow problems; however the accuracy at best will be limited since the error produced by the n-2 additions required to form ln(n!) will be the sum of the errors in the individual logarithms of the integers and will be transmitted through the exponential function to the result. These accuracy limitations are usually not a problem when the terms involving factorials all have like signs and no cancellations occur. There are instances, however, when terms involving rational fractions of factorials are of similar magnitude and alternate in sign. Such a case arises in the decomposition of products of spherical harmonics.