> >>> int(math.exp(math.lgamma(41)))
> 815915283247882431423526575245034982027017846784L
> [...]
> Prelude> fac 40
> 815915283247897734345611269596115894272000000000
Another problem is that lgamma runs in constant time but it’s accurate only in a finite range. To extend the range the polynomial must have more terms, so the “accurate” version doesn’t run in constant time. And for enough precision, you must use the “long” version of the floating point numbers. The complete analysis is complex and I don’t have enough time to do it.
The recursive functions runs in linear time only for small numbers. It uses only integer arithmetic that is faster. But if the number are big enough you must use “long” integers, and the time is probably quadratic. So it’s not clear which version is better, neither for small numbers nor for big numbers and asymptotic behavior.