10
\$\begingroup\$

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

\$\endgroup\$
  • 2
    \$\begingroup\$ 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. \$\endgroup\$ – mleyfman Jul 28 '14 at 5:22
  • \$\begingroup\$ did you look at Boost? \$\endgroup\$ – TemplateRex Jul 28 '14 at 13:15
  • \$\begingroup\$ @TemplateRex: Would that also suggest that the OP's function should return a floating-point type? \$\endgroup\$ – Jamal Jul 28 '14 at 17:28
  • \$\begingroup\$ @Jamal good point, I wrote my own for that reason. My comment was to let the OP make his requirements / Rationale explicit \$\endgroup\$ – TemplateRex Jul 28 '14 at 17:34
  • \$\begingroup\$ @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. \$\endgroup\$ – Jamal Jul 28 '14 at 17:36
6
\$\begingroup\$
  • 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.

\$\endgroup\$
  • \$\begingroup\$ 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()? \$\endgroup\$ – Jamal Jul 28 '14 at 19:43
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – gd1 Jul 28 '14 at 20:16
  • \$\begingroup\$ @gd1: That is true. I just haven't been using exception classes much yet. \$\endgroup\$ – Jamal Jul 28 '14 at 20:17
8
\$\begingroup\$

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.

\$\endgroup\$
  • 2
    \$\begingroup\$ I think using the std namespace on typedefs is unnecessarily verbose. Otherwise I agree. \$\endgroup\$ – Rapptz Jul 28 '14 at 2:05
  • 4
    \$\begingroup\$ @Rapptz: In which way does the fact that a type is actually a typedef change the way it should be qualified? \$\endgroup\$ – Nobody Jul 28 '14 at 9:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.