Calculate binomial coffeficient very reliably

Question

What is the best method to calculate a binomial coefficient in C++? I've seen some code fragments but it seemed to me that it is always only doable in some specific region. I need a very, very, very reliable calculation. I tried it with a gamma function:

unsigned n=N;
unsigned k=2;
number = tgammal(n + 1) / (tgammal(k + 1) * tgammal(n - k + 1));

but it differs already at n=8, k=2 of 1 (and by n=30, k=2 it crashes).. I "only" need a calculation of about at least n=3000 with k=2.

Answer 1

This

constexpr inline size_t binom(size_t n, size_t k) noexcept
{
    return
      (        k> n  )? 0 :          // out of range
      (k==0 || k==n  )? 1 :          // edge
      (k==1 || k==n-1)? n :          // first
      (     k+k < n  )?              // recursive:
      (binom(n-1,k-1) * n)/k :       //  path to k=1   is faster
      (binom(n-1,k) * n)/(n-k);      //  path to k=n-1 is faster
}

requires O(min{k,n-k}) operations, is reliable and can be done at compile time.

Explanation The binomial coefficients are defined as B(n,k)=k!(n-k)!/n!, from which we see that B(n,k)=B(n,n-k). We can also obtain the recurrence relations n*B(n,k)=(n-k)*B(n-1,k)=k*B(n-1,k-1). Moreover, the result is trivial for k=0,1,n,n-1.

For k=2, the result is also trivially (n*(n-1))/2.

Of course, you can do that also with other integer types. If you need to know a binomial coefficient which exceeds the largest representable integer type, you must use approximate methods: using double instead. In this case, using the gamma function is preferrable

#include <cmath>
inline double logbinom(double n, double k) noexcept
{
    return std::lgamma(n+1)-std::lgamma(n-k+1)-std::lgamma(k+1);
}
inline double binom(double n, double k) noexcept
{
    return std::exp(logbinom(n,k));
}