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.

How **accurate** does it need to be? `return 1` is very reliable, it's just entirely inaccurate . — MSalters, Jun 23 '17 at 10:23
Do you need the value for k=2 only? cause that would sure be very easy — Harshit Singhal, Jun 23 '17 at 10:37
@MSalters Good question. I just meant the output has to be correct, so I guess therefore I mean accurate? — Ben, Jun 23 '17 at 10:40
@Ben, Could you please answer? I understand that `1 + 2` isn't a binomial coefficient, but that's irrelevant to the question whether or not you consider this calculation reliable. You're asking for reliable solution, and the term "reliable" isn't clear here, so I cannot be sure that my answer is a good one. Help me help you. — alexeykuzmin0, Jun 23 '17 at 10:44
ans for k=2 would simply be n*(n-1)/2...am i missing something here? — Harshit Singhal, Jun 23 '17 at 10:45

Walter · Accepted Answer · 2018-12-19T12:30:35.003

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

I know it's a while but I came across this solution again and I think I didn't understand how it works or forgot it. How do I get to the formula "(binom(n-1,k-1) * n)/k"? And how can I call a function I'm operating in actually? — Ben, Dec 19 '18 at 12:20
@Ben. I added an explanation for the algebra. The call is recursive: the function calls itself with different argument types. Recursive algorithms are very frequent and allow for some elegance. — Walter, Dec 19 '18 at 12:32

score 4 · Answer 2 · answered Jun 23 '17 at 10:30

4

You may use asymptotically more effective recurrent 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
}

This will use only O(k) operations to calculate C(n, k).

answered Jun 23 '17 at 10:30

alexeykuzmin0

6,344
2
28
51

Is it a strange co-incidence that this code looks almost identical to that in my answer? Even the comments? Only difference is the last line which uses integer division and faster for large k. – Walter Jun 23 '17 at 10:32
@Walter, The only difference is that I use another recurrent formula: `C(n, k) = n/k * C(n - 1, k - 1)` instead of `C(n, k) = C(n - 1, k - 1) + C(n - 1, k)`. My solution runs in **O(k)** operations, while yours in **O(nk)**. So, I propose ~3000 times faster solution with no cons compared to yours. – alexeykuzmin0 Jun 23 '17 at 10:41
Indeed, my recursion is not good for large `n`. I have adopted your method but improved in case `k>n/2`. – Walter Jun 23 '17 at 10:55

Calculate binomial coffeficient very reliably

2 Answers2