Multinomial coefficients with overflow control

David Serrano Martínez

Rate me:

4.92/5 (4 votes)

15 Dec 2013CPOL3 min read

23.2K

146

There is no safety in numbers, or in anything else. (James Thurber)

Download multinomial.zip - 6.6 KB

Introduction

Calculation of multinomial coefficients is often necessary in scientific and statistic computations.

Look at this ball set:

We could wonder how many different ways we can arrange these 10 balls in a row, regarding solely ball colors and not ball numbers. Thus, the following two arrangements would be considered equal and we would just count once:

and

We have 5 blue balls, 3 red balls and 2 green balls. The solution of this problem can be noted as:

multi(5, 3, 2) = 2520

In a more general way, let:

multi(r1, r2,..., rm) a multinomial coefficient with ri non negative. Its calculation is defined as:

multi(r<sub>1</sub>, 
r<sub>2</sub>,..., 
r<sub>m</sub>) = (r<sub>1</sub> + r<sub>2</sub> +...+ r<sub>m</sub>)! / (r<sub>1</sub>! r<sub>2</sub>!... r<sub>m</sub>!)

The sign ! notes the factorial operator. The sum r1 + r2 +...+ rm is said to be the order of the multinomial coefficient.

Although conceptually very simple, its resolution can be tricky because both numerator and denominator may easily overflow even when the result itself does not.

The aim of this article is providing a fast and efficient way of computing multinomial coefficients. Possible overflows will be monitored during calculations. Thus, a exception will be thrown in case a given numeric type cannot hold the result.

Because the problem of money change is inextricably related with the calculation of multinomial coefficients, I will also address it. The description of this adjacent problem could be like this:

How many different ways can a given amount of money (rem) be gathered with coins of face value in the set {1, 2, 3,..., top}? The solution of this problem will be noted as: parti(rem, top).

Background

One of the main Internet resources on this topic is Dave Barber's Efficient Calculation of Multinomial Coefficients. His description on this subject is so comprehensive and accurate that I am just going to focus on the subtle differences between both implementations. Thus, consider this article as a complement of his (and thanks Dave for your support!).

Using the code (C++11)

C++

#include <cstdlib>
#include <iostream>
#include <stdexcept>
#include "multinomial.h"

 
int main(int argn, char *argc[])
{
    // How many different ways can I gather 21 euros with coins of face value
    // in the set: {1, 2, 3, 4 and 5 euros} (if all those coins existed)?
    try
    {
        size_t result_1=multinomial::parti<size_t>(21, 5);
    }
    catch (std::overflow_error& ex)
    {
        std::cout << ex.what() << std::endl;
    }
 
 
    // Solving the multinomial coefficient: multi(9, 8, 4).
    try
    {
        unsigned long long result_2=multinomial::multi<unsigned long long>({9, 8, 4});
    }
    catch (std::overflow_error& ex)
    {
        std::cout << ex.what() << std::endl;
    }
 
 
    system("PAUSE");
 
    return 0;
}

As you can see, a numeric type for the result must be specified in both cases as a template parameter. If such numeric types cannot hold the result, an overflow exception will be thrown.

Whereas Dave Barber's implementation relies on the client to use an ad hoc numeric type with overflow control, this implementation can use whatever built-in numeric type you wish because the overflow control is hardcoded with little overhead.

Calculations are carried out by recursion, storing the intermediate results in static vector caches. These are the recursive goberning rules:

multi(r1, r2,..., rm) = multi(r1-1, r2,..., rm) + multi(r1, r2-1,..., rm) +...+ multi(r1, r2,..., rm-1)
parti(rem, top) = parti(rem, top - 1) + parti(rem - top, top)

Please, visit Dave Barber's page for more detailed information.

Overflow control

The overflow control mechanism has been written in the file overflow.h. The code is as follows:

C++

#ifndef OVERFLOW_H
#define OVERFLOW_H
 
 
#include <limits>
#include <stdexcept>
#include <type_traits>

 
// Which numeric type is the greatest?
template <typename numeric_type_A, typename numeric_type_B>
struct greater_numeric_type
{
    typedef typename std::conditional<
        (std::numeric_limits<numeric_type_A>::digits <
            std::numeric_limits<numeric_type_B>::digits),
        numeric_type_B,
        numeric_type_A
        >::type type;
};
 
 
// Which numeric type is the smallest?
template <typename numeric_type_A, typename numeric_type_B>
struct lesser_numeric_type
{
    typedef typename std::conditional<
        std::is_same<
            typename greater_numeric_type<numeric_type_A, numeric_type_B>::type,
            numeric_type_A
            >::value,
        numeric_type_B,
        numeric_type_A
        >::type type;
};
 
 
// Does a + b overflow in the sum_numeric_type type?
template <
    typename numeric_type_A,
    typename numeric_type_B,
    typename sum_numeric_type=
        typename greater_numeric_type<numeric_type_A, numeric_type_B>::type
        >
bool is_sum_overflowed(
    const numeric_type_A& a,
    const numeric_type_B& b,
    bool is_a_overflowed=false,
    bool is_b_overflowed=false
    )
{
    // If a or b are overflowed in their own numeric types, true is returned.
    if(is_a_overflowed || is_b_overflowed)
        return true;
 
    if(a<numeric_type_A(0) || b<numeric_type_B(0))
        throw std::invalid_argument(
            "Function: is_sum_overflowed<>(...). Arguments must be non negative."
            );
 
    const sum_numeric_type& max_value=
        std::numeric_limits<sum_numeric_type>::max();
 
    // If a or b are overflowed in sum_numeric_type type, true is returned.
    if(a>max_value || b>max_value)
        return true;
 
    return max_value-a<b;
}
 
 
// Does a * b overflow in the product_numeric_type type?
// (This function is not used in the code. It is showed here for completeness).
template <
    typename numeric_type_A,
    typename numeric_type_B,
    typename product_numeric_type=
        typename greater_numeric_type<numeric_type_A, numeric_type_B>::type
        >
bool is_product_overflowed(
    const numeric_type_A& a,
    const numeric_type_B& b,
    bool is_a_overflowed=false,
    bool is_b_overflowed=false
    )
{
    // If a or b are overflowed in their own numeric types, true is returned.
    if(is_a_overflowed || is_b_overflowed)
        return true;
 
    if(a<numeric_type_A(0) || b<numeric_type_B(0))
        throw std::invalid_argument(
            "Function: is_product_overflowed<>(...). Arguments must be non negative."
            );
 
    const product_numeric_type& max_value=
        std::numeric_limits<product_numeric_type>::max();
 
    // If a or b are overflowed in product_numeric_type type, true is returned.
    if(a>max_value || b>max_value)
        return true;
 
    return max_value/a<b;
}
 
 
#endif

Overflow control can only be instantiated for numeric types supporting the std::numeric_limits<> trait. Fortunately all numeric built-in types support this trait.

This overflow control is intended to be reused in other projects. As an example, if we need to know if the sum of two numbers a and b (both of them of type size_t) overflows in a int type, we simply write:

C++

bool overflowed_sum=is_sum_overflowed<size_t, size_t, int>(a, b);

The calculations

These are carried out in multinomial.h. Although comments are in Spanish, key comments and tables are in English so understanding is not prevented at all. I encourage you to visit Dave Barber's page for further in-depth knowledge.

Every overflowable number is represented by a standard pair: the first element is the number itself and the second one is a boolean representing its state (true for overflowed and false for underflowed). This state is updated after any operation (additions).

Intermediate results are stored in static vector caches, so many subsequent calculations are immediate.

Source code

Please find attached the source code above. This has been written with Code::Blocks and compiled with MinGW GCC g++ v.4.8.1-4 in C++11.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Written By

David Serrano Martínez

Systems Engineer

Spain

I work as a senior industrial engineer for Public Administration in Spain. I have experience in developing software for ballistic computations. I like maths and programming and, above all, riding my mountain bike. Contact me at davidalvi (at gmail dot com).

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.