Arbitrary precision arithmetic with very big factorials

Question

This is a mathematical problem, not programming to be something useful!

I want to count factorials of very big numbers (10^n where n>6). I reached to arbitrary precision, which is very helpful in tasks like 1000!. But it obviously dies(StackOverflowException :) ) at much higher values. I'm not looking for a direct answer, but some clues on how to proceed further.

    static BigInteger factorial(BigInteger i)
    {
        if (i < 1)
            return 1;
        else
            return i * factorial(i - 1);
    }
    static void Main(string[] args)
    {
        long z = (long)Math.Pow(10, 12);
        Console.WriteLine(factorial(z));
        Console.Read();
    }

Would I have to resign from System.Numerics.BigInteger? I was thinking of some way of storing necessary data in files, since RAM will obviously run out. Optimization is at this point very important. So what would You recommend?

Also, I need values to be as precise as possible. Forgot to mention that I don't need all of these numbers, just about 20 last ones.

Answer 1

As other answers have shown, the recursion is easily removed. Now the question is: can you store the result in a BigInteger, or are you going to have to go to some sort of external storage?

The number of bits you need to store n! is roughly proportional to n log n. (This is a weak form of Stirling's Approximation.) So let's look at some sizes: (Note that I made some arithmetic errors in an earlier version of this post, which I am correcting here.)

(10^6)! takes order of 2 x 10^6 bytes = a few megabytes
(10^12)! takes order of 3 x 10^12 bytes = a few terabytes
(10^21)! takes order of 10^22 bytes = ten billion terabytes

A few megs will fit into memory. A few terabytes is easily within your grasp but you'll need to write a memory manager probably. Ten billion terabytes will take the combined resources of all the technology companies in the world, but it is doable.

Now consider the computation time. Suppose we can perform a million multiplications per second per machine and that we can parallelize the work out to multiple machines somehow.

(10^6)! takes order of one second on one machine
(10^12)! takes order of 10^6 seconds on one machine = 
                     10 days on one machine = 
                     a few minutes on a thousand machines.
(10^21)! takes order of 10^15 seconds on one machine = 
                        30 million years on one machine =
                        3 years on 10 million machines
                        1 day on 10 billion machines (each with a TB drive.)

So (10^6)! is within your grasp. (10^12)! you are going to have to write your own memory manager and math library, and it will take you some time to get an answer. (10^21)! you will need to organize all the resources of the world to solve this problem, but it is doable.

Or you could find another approach.

Answer 2

The solution is easy: Calculate the factorials without using recursion, and you won't blow out your stack.

I.e. you're not getting this error because the numbers are too large, but because you have too many levels of function calls. And fortunately, for factorials there's no reason to calculate them recursively.

Once you've solved your stack problem, you can worry about whether your number format can handle your "very big" factorials. Since you don't need the exact values, use one of the many efficient numeric approximations (which you can count on to get all of the most significant digits right). The most common one is Stirling's approximation:

n! ~ n^n e^{-n} sqrt(2 \pi n)

The image is from this page, where you'll find discussion and a second, more accurate formula (although "in most cases the difference is quite small", they say). Of course this number is still too large for you to store, but now you can work with logarithms and drop the unimportant digits before you extract the number. Or use the Wikipedia version of the approximation, which is already expressed as a logarithm.

Answer 3

Unroll recursion:

static BigInteger factorial(BigInteger n)
{
    BigInteger res = 1;
    for (BigInteger i = 2; i <= n; ++i)
        res *= i;
    return res;
}