How can I calculate a factorial in C# using a library call?

Question

I need to calculate the factorial of numbers up to around 100! in order to determine if a series of coin flip-style data is random, as per this Wikipedia entry on Bayesian probability. As you can see there, the necessary formula involves 3 factorial calculations (but, interestingly, two of those factorial calculations are calculated along the way to the third).

I saw this question here, but I'd think that integer is going to get blown out pretty quickly. I could also make a function that is more intelligent about the factorial calculation (ie, if I have 11!/(7!3!), as per the wiki example, I could go to (11*10*9*8)/3!), but that smacks of premature optimization to me, in the sense that I want it to work, but I don't care about speed (yet).

So what's a good C# library I can call to calculate the factorial in order to get that probability? I'm not interested in all the awesomeness that can go into factorial calculation, I just want the result in a way that I can manipulate it. There does not appear to be a factorial function in the Math namespace, hence the question.

Answer 1

There has been a previous question on a similar topic. Someone there linked the Fast Factorial Functions web site, which includes some explanations of efficient algorithms and even C# source code.

Answer 2

You could try Math.NET - I haven't used that library, but they do list Factorial and Logarithmic Factorial.

Answer 3

Do you want to calculate factorials, or binomial coefficients?

It sounds like you want to calculate binomial coefficients - especially as you mention 11!/(7!3!).

There may be a library that can do this for you, but as a (presumably) programmer visiting stack overflow there's no reason not to write one yourself. It's not too complicated.

To avoid memory overflow, don't evaluate the result until all common factors are removed.

This algorithm still needs to be improved, but you have the basis for a good algorithm here. The denominator values need to be split into their prime factors for the best result. As it stands, this will run for n = 50 quite quickly.

float CalculateBinomial(int n, int k)
{
    var numerator = new List<int>();
    var denominator = new List<int>();
    var denominatorOld = new List<int>();

    // again ignore the k! common terms
    for (int i = k + 1; i <= n; i++)
        numerator.Add(i);

    for (int i = 1; i <= (n - k); i++)
    {
        denominator.AddRange(SplitIntoPrimeFactors(i));
    }

    // remove all common factors
    int remainder;                
    for (int i = 0; i < numerator.Count(); i++)
    {
        for (int j = 0; j < denominator.Count() 
            && numerator[i] >= denominator[j]; j++)
        {
            if (denominator[j] > 1)
            {
                int result = Math.DivRem(numerator[i], denominator[j], out remainder);
                if (remainder == 0)
                {
                    numerator[i] = result;
                    denominator[j] = 1;
                }
            }
        }
    }

    float denominatorResult = 1;
    float numeratorResult = 1;

    denominator.RemoveAll(x => x == 1);
    numerator.RemoveAll(x => x == 1);

    denominator.ForEach(d => denominatorResult = denominatorResult * d);
    numerator.ForEach(num => numeratorResult = numeratorResult * num);

    return numeratorResult / denominatorResult;
}

static List<int> Primes = new List<int>() { 2, 3, 5, 7, 11, 13, 17, 19, 
    23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };

List<int> SplitIntoPrimeFactors(int x)
{
    var results = new List<int>();
    int remainder = 0;

    int i = 0;
    while (!Primes.Contains(x) && x != 1)
    {
        int result = Math.DivRem(x, Primes[i], out remainder);
        if (remainder == 0)
        {
            results.Add(Primes[i]);
            x = result;
            i = 0;
        }
        else
        {
            i++;
        }
    }
    results.Add(x);
    return results;
}

I can estimate n = 110, k = 50 (returns 6x10^31) but cannot run n = 120, k = 50.