A fast and accurate way to approximate 
n!  for large n.

We have the approximation in the next cell from equation 6.1.38 in Abramowitz M. & Stegun I. A. 1972 Handbook of mathematical functions with formulas, graphs and mathematical tables. New York, NY: Dover

n! ≈ (Sqrt[2  π]) (n^(n + 
              1/2)) (E^(-n))

An "more useful" approximation was found by finding a rational function that is a very good approximation for the expression in the next cell.

((((n^(n + 1/2)) (E^(-n))))/(n!))

The rational approximation for the previous cell is used inplace of (Sqrt[2  π]) in the function NFactorial below. 
In this definition of NFactorial I avoid undesired overflow by using  Exp[Log[n](n+1/2)-n]  instead of ((n^(n + 1/2)) (e^(-n))).

Off[General::spell1];

ClearAll[NFactorial];

Attributes[NFactorial] = {Listable, 
          NumericFunction};

NFactorial[
        n_Integer?
          Positive] := 
With[{nn = 
            SetPrecision[n, 25]}, 
SetPrecision[
          Exp[Log[nn]*(nn + 0.5`50) - 
                nn]*(0.252954663446343580963342351 + 
                  n*(-25.34877082357403134347634 + 
                        40.68750025249009458918000*
                          n))/
(1 + 
                  n*(-11.465360115610115571696880 + 
                        16.231964134562266197858360*n)), 
          21.0]
]

I found that this approximation for Factorial is correct to 21 decimal places for all n greater than 49106.  In the cells below NFactorial tries to approximate Factorial at 44000000 and 45000000 respectively. I am using Mathematica 4 in which case Factorial[45000000] is larger than $MaxNumber. However, if you are using Mathematica 6 Factorial will not overflow until about Factorial[86000000].

NFactorial[44000000]

NFactorial[45000000]

In many cases my NFactorial algorithm is faster than the built-in algorithm for approximating Factorial[n]. Using Mathematica 4 my algorithm is faster than the built-in algorithm for all  (n > 60000). I don't have the results of extensive timing tests using Mathematica 6, but one user told me it took 1.28 Seconds to evaluate the next cell using Mathematica 6.

N[Factorial[800000]] // Timing

The timing tests below suggest that computing  NFactorial[n] with Mathematica 4  takes less than 0.000441 Seconds for essentially any integer over 60000.

First[Timing[NFactorial[Range[60000, 60000 + 10000]]]]/10000

First[Timing[NFactorial[Range[44000000, 44000000 + 10000]]]]/10000

The next cell shows how close the rational function used above gets to (Sqrt[2  π])  in the limit as  n->∞.  When using Mathematica 4, (n!) results in an overflow for (n) greater than 44787863.  I haven't tried to see how well  NFactorial[n]  approximates Factorial for (n) larger than  44787863, but I suspect it is correct to 21 decimal places for all (n) greater than 49105.

ration = (0.252954663446343580963342351 + 
              n*(-25.34877082357403134347634 + 
                    40.68750025249009458918000*
                      n))/
(1 + 
              n*(-11.465360115610115571696880 + 
                    16.231964134562266197858360*
                      n));

Limit[ration - (Sqrt[2  π]), {n -> ∞}]