At 20:50:54 Mon, Dec 8 1997, I gave an improved Sirling's approximation - z - y z 2 1 z! ~ S(z) := sqrt(2 pi) e (z + y) sqrt(z + y + -) 6 where y is somewhat arbitrary. The case y=0 closely resembles (and thoroughly trounces) the traditional formula. However, there is a clearly best value of y, ~.96841318776609465, that minimizes max |log(S(z)/z!)|. 0<=z<oo The worst cases are at z=0 and z = w :~ .7662945329650573, [w,y] being a solution of | 2 1 | log (z + y + -) | d 6 | [-- (log(z!) - ---------------- - z log(z + y) + z + y) | , dz 2 | |z = w 2 1 2 1 log((y + -) (y + w + -)) 6 6 w! - -------------------------- - w log(y + w) + 2 y + log(-----) + w]. 2 2 %pi S(0)/0! ~ 1.00021756216683 ~ w!/S(w), i.e., .02% max error. z! = round(S(z)) for z = 0..7. S(8) ~ 40319.05, which could be fixed by a parity argument. It might be possible to combine this S(z) approximation with fancier divisibility arguments (and a slightly different, possibly variable y) to get a competitive (exact) algorithm for large integer z. For a numerical Gamma function, one can choose the y (e.g. .545525525713106) that minimaxes error (in this case, 3e-7) over, in this case, 9 <= z < oo. Then for small z use the standard trick, z! = (z+9)!/(z+9)/.../(z+1). Strangely, there is another, minutely inferior y ~ .9523531774380170, with max error near w=15.79 vs 15.45 . --rwg