# Binomial coefficient

I wrote a C++ function to calculate the binomial coefficient, trying to avoid overflows as much as possible.

/*!
* @brief Calculates the binomial coefficient indexed by n and k
*
* This implementation is based on the algorithm described here: http://blog.plover.com/math/choose.html
* and the relevant follow-up: http://blog.plover.com/math/choose-2.html
*
* @param n                     The number of elements
* @param k                     The number of elements in each subset
* @return                      The number of k-element subsets in a n-element set
*
* @throw std::overflow_error   If the computation caused an overflow
*/
uint64_t binomial_coefficient(uint64_t n, uint64_t k)
{
if (k > n)
return 0;

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

uint64_t r = 1;

for (uint64_t d = 1; d <= k; d++)
{
uint64_t mult = n;

bool divided = true;

if (mult % d == 0)
mult /= d;
else if (r % d == 0)
r /= d;
else
divided = false;

const uint64_t r_mult = r * mult;
if (r_mult / mult != r)
throw std::overflow_error("Overflow");

r = r_mult;

if (!divided)
r /= d;

n--;
}

return r;
}


I have implemented the algorithm described here with a slight modification (try to divide n or r by d before multiplying n by r) as described in the article follow-up.

I also added an overflow check, which I believe should be reliable, as unsigned arithmetic overflow should not yield undefined behaviour.

Do you see any blatant or subtle bug in my implementation, and/or any room for improvement (besides lookup tables for common values and alike, which I am not willing to pursue at the moment)?

• Just so you know, binomial coefficients have a very naturally occurring lookup table: Pascal's Triangle. Using the triangle is actually a better way to calculate binomial coefficients for smaller n and k, and you get memoization for free if you implement it right. Commented Jul 28, 2014 at 5:22
• did you look at Boost? Commented Jul 28, 2014 at 13:15
• @TemplateRex: Would that also suggest that the OP's function should return a floating-point type? Commented Jul 28, 2014 at 17:28
• @Jamal good point, I wrote my own for that reason. My comment was to let the OP make his requirements / Rationale explicit Commented Jul 28, 2014 at 17:34
• @TemplateRex: I was trying to come up with something simpler (not incorporating Boost since I'm not too familiar with it), but I couldn't quite think of anything good. Commented Jul 28, 2014 at 17:36

• An exception message should tell some details of the overflow (at least values of n and k)

• The function begs to be templated.

The rest is not a review, but rather a long comment.

I don't think the algorithm achieves the declared goal (to avoid overflows as much as possible): an overflow occurs in the intermediate result (r_mult), even though an actual result (r) could very well fit uint64_t after division.

A root of the problem lies in an opportunistic approach. Factors are reduced only if one of them is divisible by d. The fact is, they can be reduced even if they are not. Surely they are divisible by a greatest common factor:

    uint64_t result = 1;

for (uint64_t d = 1; d <= k; d++)
{
uint64_t num = n - k + d;
uint64_t denom = d;

uint64_t common_factor = gcd(result, denom);

result /= common_factor;
denom /= common_factor;

assert(num % denom == 0);
num /= denom;

if (result > UINT64_MAX / num) {
throw std::overflow_error("Overflow");
}

result *= num;
}

return result;


PS: it could be mathematically proven that the assertion never triggers.

• Regarding the first point, would the if block here also need to display n and k (probably via std::cerr)? Or can all of that be passed into std::overflow_error()? Commented Jul 28, 2014 at 19:43
• vnp: thank you! @Jamal: I do think that one shall never print anything in a function supposed to do... computations (don't touch my streams!). Exceptions exist for the very purpose of handling errors at the right time, don't they.
– gd1
Commented Jul 28, 2014 at 20:16
• @gd1: That is true. I just haven't been using exception classes much yet. Commented Jul 28, 2014 at 20:17

Mostly minor things, but could still be helpful:

• Functions should generally be named in a verb form, so rename the function to something like calculate_binomial_coefficient().

• This is a little hard to read:

if (r_mult / mult != r)


and could use more parenthesis for clarity:

if ((r_mult / mult) != r)

• Consider using curly braces for these single-line conditional statements. It could also help with avoiding some additional bugs, if you end up adding anything to them.

• Since this is C++ and not C, use std::uint64_t. Try to use all the equivalents that are found in the std namespace.

• I think using the std namespace on typedefs is unnecessarily verbose. Otherwise I agree. Commented Jul 28, 2014 at 2:05
• @Rapptz: In which way does the fact that a type is actually a typedef change the way it should be qualified? Commented Jul 28, 2014 at 9:48