Binomial coefficient

Question

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)?

Answer 1

• 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.

Answer 2

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.