Calculator of combinations without repetition

Question

I am interested to find an absolute value (not an approximation) of "combination without repetition" for given \$n\$ and \$k\$, or \$\binom{n}{k}\$.

The brute force solution would look like this

private static ulong Factorial(int x)
{
    ulong res = 1;
    while (x > 1)
    {
        res *= (ulong)x--;
    }
    return res;
}

public static int Combination0(int k, int n)
{
    k = Math.Min(k, n - k);
    if (n < 2 || k < 1) return 1;
    if (k == 1) return n;
    return (int)(Factorial(n) / (Factorial(k) * Factorial(n - k)));
}

We can slightly optimize this solution, by finding \$\prod_{n\geq i>k}{i}\$ instead of \$\frac{n!}{(n-k)!}\$.

private static ulong Factorial(int x, int until = 0)
{
    ulong res = 1;
    while (x > until)
    {
        res *= (ulong)x--;
    }
    return res;
}

public static int Combination1(int k, int n)
{
    k = Math.Min(k, n - k);
    if (n < 2 || k < 1) return 1;
    if (k == 1) return n;
    return (int)(Factorial(n, n - k) / Factorial(k));
}

But these two solutions have one significant problem - we are limited by ulong.MaxValue, which is more than \$20!\$, but less than \$21!\$.

Another way to find the number of combinations, which doesn't have the previously described problem, is the Pascal's triangle.

public static int Combination2(int k, int n)
{
    k = Math.Min(k, n - k);
    if (n < 2 || k < 1) return 1;
    if (k == 1) return n;
    int[] triangle = new int[k + 1];
    triangle[0] = 1;

    // expanding
    int i = 0;
    for (; i < k; i++)
    {
        for (int j = i + 1; j > 0; j--)
        {
            triangle[j] += triangle[j - 1];
        }
    }

    // progressing
    for (; i < n - k; i++)
    {
        for (int j = k; j > 0; j--)
        {
            triangle[j] += triangle[j - 1];
        }
    }

    // collapsing
    for (; i < n; i++)
    {
        int until = k - (n - i);
        for (int j = k; j > until; j--)
        {
            triangle[j] += triangle[j - 1];
        }
    }
    return triangle[k];
}

But the problem is that Combination2 is significantly slow.

I would appreciate any comments and suggestions for an improvement.

---

Update

@quasar and @henrik-hansen suggested the way to prevent overflow by calculating \$\prod_{0 \leq i < k}{\frac{n-i}{i+1}}\$.

Answer 1

You should forget about factorial when it comes to Combinations (n, k). Instead you can use the formula: n(n-1)(n-2)...(n-k+1)/(1*2*3*...*k). You start with n and then iterate over x = 1 .. k - 1 and successively multiply with (n-x) and at the same time reduce by dividing with x. All in all it ends up like this:

public ulong Combinations(ulong n, ulong k)
{
  ulong count = n;

  for (ulong x = 1; x <= k - 1; x++)
  {
    count = count * (n - x) / x;
  }

  return count / k;
}

In this way you prevent overflow from intermediate factorial calculations.

Answer 2

One trick is to keep the partial products in small, \${n}\choose{k}\$ so they don't overflow.

I iteratively multiply \$n/(n-k)\$ by \$(n-1)/(n-k-1)\$, cache the result in an accumulator, multiply that by \$(n-2)/(n-k-2)\$ and so forth.

#include <iostream>

template <class T>
T choose(T n, T k)
{
    T accum = 1;
    T m = n;
    for (T i = 1; i <= m - k; i++)
    {
        accum = accum * n / (n - k);
        n--;
    }
    return accum;
}

int main()
{
    std::cout << std::fixed;
    long double n = 50, k = 25;
    std::cout << "\nLDBL_MAX" << LDBL_MAX;
    long double result = choose(n, k);
    std::cout << "\nC(" << n << "," << k << ") = " << result;
    std::cin.clear();
    std::cin.get();
    return 0;
}

Answer 3

I think the question of which method is faster depends on your intended use. FWIW I've found that if I use any C(N,k) in a program then I tend to use a lot of them. In this case the fastest way can be to make a table of all the C(N,k) (0<=k<=N, N<=Nmax) at program start, and then the speed of the routine for calculating C(N,k) is less important. Anyway, I think that the fastest way to do this is to use the recurrence relation; my tests show it's around 100 times faster than Henrik's method for Nmax = 57.

The recurrence relation method is also able to compute a few more rows than the others. I reckon overflow strikes Henrik's method around N=58, whereas the recurrence relation is good up to N=68

Finally in production code I would want to see at least a claim, ideally a reference to a proof, that where there are integer divisions the left hand side is indeed a multiple of the right hand side. This is the case in Henrik's method, but as Henrik points out, not in Quasars's. The problem with not having such a claim is that, say in debugging, someone might stare at the code and wonder if they are losing accuracy.