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.
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));
}