The binomial coefficient calculation is very reliable.

Question

What is the best method for calculating a binomial coefficient in C ++? I saw code fragments, but it seemed to me that it is always possible only in a certain region. I need a very, very, very reliable calculation. I tried this with the gamma function:

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

but it already differs at n = 8, k = 2 from 1 (and at n = 30, k = 2 it falls). I only need to calculate at least n = 3000 for k = 2.

Answer 1

this is

 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)/(nk); // path to k=n-1 is faster } 

It requires O (min {k, nk}) operations, is reliable, and can be executed at compile time.

Explanation Binomial coefficients are defined as B(n,k)=k!(nk)!/n! whence we see that B(n,k)=B(n,nk) . We can also obtain the recurrence relations n*B(n,k)=(nk)*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 result is also trivial (n*(n-1))/2 .

Of course, you can do this with other integer types as well. If you need to know a binomial coefficient that exceeds the largest representable integer type, you should use approximate methods: use double instead. In this case, the use of a gamma function is preferred.

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

Answer 2

You can use an asymptotically more efficient recurrence formula:

 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 binom(n - 1, k - 1) * n / k; // recursive } 

To calculate C (n, k), only the operations O (k) are used .