Calculating the Amount of Combinations

From Andreas' answer:

Here's an ancient algorithm which is exact and doesn't overflow unless the result is to big for a long long

unsigned long long
choose(unsigned long long n, unsigned long long k) {
    if (k > n) {
        return 0;
    }
    unsigned long long r = 1;
    for (unsigned long long d = 1; d <= k; ++d) {
        r *= n--;
        r /= d;
    }
    return r;
}

This algorithm is also in Knuth's "The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms" I think.

UPDATE: There's a small possibility that the algorithm will overflow on the line:

r *= n--;

for very large n. A naive upper bound is sqrt(std::numeric_limits<long long>::max()) which means an n less than rougly 4,000,000,000.

Consider n == 67 and k == 33. The above algorithm overflows with a 64 bit unsigned long long. And yet the correct answer is representable in 64 bits: 14,226,520,737,620,288,370. And the above algorithm is silent about its overflow, choose(67, 33) returns:

8,829,174,638,479,413

A believable but incorrect answer.

However the above algorithm can be slightly modified to never overflow as long as the final answer is representable.

The trick is in recognizing that at each iteration, the division r/d is exact. Temporarily rewriting:

r = r * n / d;
--n;

For this to be exact, it means if you expanded r, n and d into their prime factorizations, then one could easily cancel out d, and be left with a modified value for n, call it t, and then the computation of r is simply:

// compute t from r, n and d
r = r * t;
--n;

A fast and easy way to do this is to find the greatest common divisor of r and d, call it g:

unsigned long long g = gcd(r, d);
// now one can divide both r and d by g without truncation
r /= g;
unsigned long long d_temp = d / g;
--n;

Now we can do the same thing with d_temp and n (find the greatest common divisor). However since we know a-priori that r * n / d is exact, then we also know that gcd(d_temp, n) == d_temp, and therefore we don't need to compute it. So we can divide n by d_temp:

unsigned long long g = gcd(r, d);
// now one can divide both r and d by g without truncation
r /= g;
unsigned long long d_temp = d / g;
// now one can divide n by d/g without truncation
unsigned long long t = n / d_temp;
r = r * t;
--n;

Cleaning up:

unsigned long long
gcd(unsigned long long x, unsigned long long y)
{
    while (y != 0)
    {
        unsigned long long t = x % y;
        x = y;
        y = t;
    }
    return x;
}

unsigned long long
choose(unsigned long long n, unsigned long long k)
{
    if (k > n)
        throw std::invalid_argument("invalid argument in choose");
    unsigned long long r = 1;
    for (unsigned long long d = 1; d <= k; ++d, --n)
    {
        unsigned long long g = gcd(r, d);
        r /= g;
        unsigned long long t = n / (d / g);
        if (r > std::numeric_limits<unsigned long long>::max() / t)
           throw std::overflow_error("overflow in choose");
        r *= t;
    }
    return r;
}

Now you can compute choose(67, 33) without overflow. And if you try choose(68, 33), you'll get an exception instead of a wrong answer.

The following routine will compute the n-choose-k, using the recursive definition and memoization. The routine is extremely fast and accurate:

inline unsigned long long n_choose_k(const unsigned long long& n,
                                     const unsigned long long& k)
{
   if (n  < k) return 0;
   if (0 == n) return 0;
   if (0 == k) return 1;
   if (n == k) return 1;
   if (1 == k) return n;       
   typedef unsigned long long value_type;
   value_type* table = new value_type[static_cast<std::size_t>(n * n)];
   std::fill_n(table,n * n,0);
   class n_choose_k_impl
   {
   public:

      n_choose_k_impl(value_type* table,const value_type& dimension)
      : table_(table),
        dimension_(dimension)
      {}

      inline value_type& lookup(const value_type& n, const value_type& k)
      {
         return table_[dimension_ * n + k];
      }

      inline value_type compute(const value_type& n, const value_type& k)
      {
         if ((0 == k) || (k == n))
            return 1;
         value_type v1 = lookup(n - 1,k - 1);
         if (0 == v1)
            v1 = lookup(n - 1,k - 1) = compute(n - 1,k - 1);
         value_type v2 = lookup(n - 1,k);
         if (0 == v2)
            v2 = lookup(n - 1,k) = compute(n - 1,k);
         return v1 + v2;
      }

      value_type* table_;
      value_type dimension_;
   };
   value_type result = n_choose_k_impl(table,n).compute(n,k);
   delete [] table;
   return result;
}

Using a dirty trick with a long double, it is possible to get the same accuracy as Howard Hinnant (and probably more):

unsigned long long n_choose_k(int n, int k)
{
    long double f = n;
    for (int i = 1; i<k+1; i++)
        f /= i;
    for (int i=1; i<k; i++)
        f *= n - i;

    unsigned long long f_2 = std::round(f);

    return f_2;
}

The idea is to divide first by k! and then to multiply by n(n-1)...(n-k+1). The approximation through the double can be avoided by inverting the order of the for loop.

Improves Howard Hinnant's answer (in this question) a little bit: Calling gcd() per loop seems a bit slow. We could aggregate the gcd() call into the last one, while making the most use of the standard algorithm from Knuth's book "The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms":

const uint64_t u64max = std::numeric_limits<uint64_t>::max();
uint64_t choose(uint64_t n, uint64_t k)
{
    if (k > n)
        throw std::invalid_argument(std::string("invalid argument in ") + __func__);

    if (k > n - k)
        k = n - k;

    uint64_t r = 1;
    uint64_t d;
    for (d = 1; d <= k; ++d) {
        if (r > u64max / n)
            break;
        r *= n--;
        r /= d;
    }

    if (d > k)
        return r;

    // Let N be the original n,
    // n is the current n (when we reach here)
    // We want to calculate C(N,k),
    // Currently we already calculated the r value so far:
    // r = C(N, n) = C(N, N-n) = C(N, d-1)
    // Note that N-n = d-1
    // In addition we know the following identity formula:
    //  C(N,k) = C(N,d-1) * C(N-d+1, k-d+1) / C(k, k-d+1)
    //         = C(N,d-1) * C(n, k-d+1) / C(k, k-d+1)
    // Using this formula, we effectively reduce the calculation,
    // while recursively use the same function.
    uint64_t b = choose(n, k-d+1);
    if (b == u64max) {
        return u64max;  // overflow
    }

    uint64_t c = choose(k, k-d+1);
    if (c == u64max) {
        return u64max;  // overflow
    }

    // Now, the combinatorial should be r * b / c
    // We can use gcd() to calculate this:
    // We Pick b for gcd: b < r almost (if not always) in all cases
    uint64_t g = gcd(b, c);
    b /= g;
    c /= g;
    r /= c;

    if (r > u64max / b)
        return u64max;   // overflow

    return r * b;
}

Note that the recursive depth is normally 2 (I don't really see a case goes to 3, the combinatorial reducing is quite decent.), i.e. calling choose() for 3 times, for non-overflow cases.

Replace uint64_t with unsigned long long if you prefer it.

One of SHORTEST way :

int nChoosek(int n, int k){
    if (k > n) return 0;
    if (k == 0) return 1;
    return nChoosek(n - 1, k) + nChoosek(n - 1, k - 1);
}