Expanding a bit on Gene's comments ... N! = LCM (or product) of the individual prime powers P^E for P<=N. The exponent E = ( N - digitsum(N base P) ) / (P-1). The portion of the LCM attributable to P > N/2 (where each E=1) is roughly e^(N/2), somewhat less than the approx 2^N attributable to P=2. Any method for fixing up an exact N! will either need to use most of the primes <N, or have a good fraction of the high- order bits known. The usual Stirling's formula needs (very very roughly) N terms (and N Bernoulli numbers) for O(N) or O(N logN) bits, so it's not much savings for computing N! exactly. But having a good way to get 10 or 20 digits of X! is handy. Since the formulas always seem to include sqrts, maybe the right approach is to try approximating N!^2? Rich ________________________________ From: math-fun-bounces+rschroe=sandia.gov@mailman.xmission.com on behalf of R. William Gosper Sent: Sat 10/22/2005 1:00 PM To: math-fun@mailman.xmission.com Subject: [math-fun] another Stirling's tweak 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 _______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun