Efficient calculation of multinomial coefficients.
Version of Friday 13 December 2013.
Dave Barber's other pages.

Thanks to Ralf Anderson for error detection. Also see David Serrano's version.

 Quick reference. The C++ code is here; you can copy and paste it into your favorite editor. To find a multinomial coefficient, invoke: ` multinomial::multi(std::vector const & vec)` where the components of vec are the arguments to the multinomial function. The return value is of type result_type. To find the number of partitions of an integer rem, where no element is larger than top, invoke: ` multinomial::parti(size_t const rem, size_t const top)` The return value is of type size_t.

Detailed explanation.

When mathematicians solve problems in combinatorics, they often need to calculate multinomial coefficients. Unfortunately, the most direct definition of these numbers requires factorials, which are often inconveniently large.

The multinomial coefficient function, here written multi (…), takes several arguments which are nonnegative integers. It is defined in plain English as

 the factorial of the sum of the arguments divided by the product of the factorials of the individual arguments

For instance, with three or five arguments:

• multi (a, b, c) ≡ (a + b + c)! ÷ (a! × b! × c!)
• multi (a, b, c, d, e) ≡ (a + b + c + d + e)! ÷ (a! × b! × c! × d! × e!)
The zero-argument case is handled by convention: multi ( ) ≡ 1.

The sum of the arguments is the order of the multi (…) invocation.

This function is commutative: the sequence in which the arguments are written makes no difference.

The purpose of this web page is to provide a efficient scheme for calculating multi (…) in the computer programming language C++. A similar implementation can be developed in other languages.

A numerical overflow problem can be demonstrated with as few as three arguments. For instance:

• multi (7, 4, 2) = (7 + 4 + 2)! ÷ (7! × 4! × 2!) = 6,227,020,800 ÷ 241,920 = 25,740

Even though the result (25,740) is of reasonable size, the numerator (6,227,020,800) will not fit into the 32-bit unsigned integer offered on many computers, where the maximum value is 4,294,967,295.

Of little help is the 64-bit unsigned integer available on most newer machines, because it has a maximum of 18,446,744,073,709,551,615, which is less than 21!, as seen here:

• multi (9, 8, 4) = (9 + 8 + 4)! ÷ (9! × 8! × 4!) = 51,090,942,171,709,440,000 ÷ 351,151,718,400 = 145,495,350

The result fits into a 32-bit integer, the denominator into a 64-bit, but the numerator into neither.

Extended-precision integer packages are readily available, and they solve the overflow problem, but often at a considerable speed penalty. Writing computer code to cancel common factors from numerator and denominator is possible, but complicated, and offers little hope of fast performance.

The key to efficient calculation is to exploit a recursive property of multinomial coefficients. Although it is awkward to write in the general case, the four-argument formula as an example ought to make it clear:

• multi (a, b, c, d) = multi (a − 1, b, c, d) + multi (a, b − 1, c, d) + multi (a, b, c − 1, d) + multi (a, b, c, d − 1)

More concretely:

• multi (5, 9, 8, 4) = multi (4, 9, 8, 4) + multi (5, 8, 8, 4) + multi (5, 9, 7, 4) + multi (5, 9, 8, 3)

If an invocation of multi (…) would have a negative argument, the invocation is replaced by zero, and recursion of that branch ends.

The first virtue of this scheme is that no intermediate result requires larger storage than the ultimate answer. A second benefit is that calculations involve neither multiplication nor division, but only addition — computers can typically add faster than they can multiply, and multiply faster than they can divide. If a large result is anticipated and extended-precision integers are chosen, the time savings with recursion can be substantial.

One disadvantage of the recursive formula, if it is implemented naïvely, is that it leads to extensive repetition in the calculation. For instance, both multi (4, 9, 8, 4) and multi (5, 8, 8, 4) need multi (4, 8, 8, 4) as an intermediate result. This redundancy can be eliminated, however, by establishing a cache.

In this implementation the cache is initialized with the multinomial coefficient of order zero, namely multi ( ) = 1. Later, when a multinomial coefficient of order n is requested,

• If the required multinomial coefficient is already in the cache, no calculation is necessary; it is simply retrieved.
• If the required multinomial coefficient is not in the cache, then all the multinomial coefficients of order n are calculated and encached — but only after ensuring that all multinomial coefficients of lower order are in the cache.

For instance, if the first request is for multi (2, 1, 1), which is of order 4, the following will be calculated and encached, in this sequence:

• order 1: multi (1)
• order 2: multi (1, 1), multi (2)
• order 3: multi (1, 1, 1), multi (2, 1), multi (3)
• order 4: multi (1, 1, 1, 1), multi (2, 1, 1), multi (2, 2), multi (3, 1), multi (4)

If a second request is for multi (1, 1, 1), no calculation is necessary.

If a third request is for multi (3, 2), the seven multinomial coefficients of order 5 will be calculated and encached, but the multinomial coefficients of order 4 and below will not need to be recalculated.

The first 45 entries in the multi cache, if that many are needed, look like this:

Cache
Value
stored
Multinomial
coefficient of …
01( )
11(1)
22(1, 1)
31(2)
46(1, 1, 1)
53(2, 1)
61(3)
724(1, 1, 1, 1)
812(2, 1, 1)
96(2, 2)
104(3, 1)
111(4)
12120(1, 1, 1, 1, 1)
1360(2, 1, 1, 1)
1430(2, 2, 1)
Cache
Value
stored
Multinomial
coefficient of …
1520(3, 1, 1)
1610(3, 2)
175(4, 1)
181(5)
19720(1, 1, 1, 1, 1, 1)
20360(2, 1, 1, 1, 1)
21180(2, 2, 1, 1)
2290(2, 2, 2)
23120(3, 1, 1, 1)
2460(3, 2, 1)
2520(3, 3)
2630(4, 1, 1)
2715(4, 2)
286(5, 1)
291(6)
Cache
Value
stored
Multinomial
coefficient of …
305040(1, 1, 1, 1, 1, 1, 1)
312520(2, 1, 1, 1, 1, 1)
321260(2, 2, 1, 1, 1)
33630(2, 2, 2, 1)
34840(3, 1, 1, 1, 1)
35420(3, 2, 1, 1)
36210(3, 2, 2)
37140(3, 3, 1)
38210(4, 1, 1, 1)
39105(4, 2, 1)
4035(4, 3)
4142(5, 1, 1)
4221(5, 2)
437(6, 1)
441(7)

The C++ function to invoke multi (…) is called multinomial::multi is and is declared

```    namespace multinomial {
template <typename result_type>
result_type multi(std::vector<size_t> const &)
}```

Presumably, result_type will be an unsigned integer, perhaps of extended precision. For many purposes size_t will be satisfactory, even though its range can vary from machine to machine. A disadvantage of size_t, however, is that the programmer cannot expect to be alerted in case of overflow. If exactitude is not required, result_type could instead be floating point.

Any zeroes that happen to be in the argument std::vector<size_t> do not affect the answer. This turns out to be convenient because when a std::vector<size_t> is constructed, its components are by default set to size_t(0) which equals zero.

Among other things, namespace multinomial contains these items for convenience:

```    namespace multinomial {
typedef std::vector<size_t> SVI;
void view (std::ostream &, SVI const &);
}```

With those declarations, this sample code:

```    using namespace multinomial;
SVI v(4);
v.at(0) = 7;
v.at(1) = 9;
v.at(2) = 4;
v.at(3) = 6;
std::cout << "\n multi ";
view (std::cout, v);
std::cout << " = " << multi<unsigned long long int>(v);```
gives this output:
`    multi (7, 9, 4, 6) = 12760912164000`

`    multi<double>(v)`
an answer such as 1.27609e+13 would have ensued.

The cache is stored in a std::vector<result_type> in class multinomial::combo. Within the vector, multinomial coefficients of order n are stored at lower addresses than multinomial coefficients of order n + 1, as shown in the table above. The vector grows as needed, and it is not necessary for the programmer to predict its ultimate size.

Efficient subscripting of this vector is a nontrivial matter, because the multinomial coefficients of order n are of the same quantity as the integer partitions of n, and partition quantity functions are complicated to figure. Class multinomial::index uses recursion to generate the values, and stores them in its own cache, which is separate from that of class multinomial::combo.

Programmers who are interested only in finding multinomial coefficients may ignore class multinomial::index. However, programmers who are enumerating partitions can easily use, through a wrapper function parti, the code of class index that calculates partition quantities:

```    namespace multinomial {
size_t parti (size_t const rem, size_t const top);
}```

This tells how many partitions of the integer rem exist, with no component greater than top. For example multinomial::parti (7, 4) equals 11, because there are 11 ways to break 7 down into addends that lie between 1 and 4:

• 1 + 1 + 1 + 1 + 1 + 1 + 1
• 2 + 1 + 1 + 1 + 1 + 1
• 2 + 2 + 1 + 1 + 1
• 2 + 2 + 2 + 1
• 3 + 1 + 1 + 1 + 1
• 3 + 2 + 1 + 1
• 3 + 2 + 2
• 3 + 3 + 1
• 4 + 1 + 1 + 1
• 4 + 2 + 1
• 4 + 3

The recursion formula is

`    parti (rem, top) = parti (rem - top, top) + parti (rem, top - 1)`
with two conditions:
• If top > rem, then parti (rem, rem) is substituted for parti (rem, top).
• If top = 0, then 0 is substituted for parti (rem, top), and that branch of the recursion ends.

The first 45 entries in the parti cache, if that many are needed, look like this:

Cache
Value
stored
(rem,
top)
01(0, 0)
10(1, 0)
21(1, 1)
30(2, 0)
41(2, 1)
52(2, 2)
60(3, 0)
71(3, 1)
82(3, 2)
93(3, 3)
100(4, 0)
111(4, 1)
123(4, 2)
134(4, 3)
145(4, 4)
Cache
Value
stored
(rem,
top)
150(5, 0)
161(5, 1)
173(5, 2)
185(5, 3)
196(5, 4)
207(5, 5)
210(6, 0)
221(6, 1)
234(6, 2)
247(6, 3)
259(6, 4)
2610(6, 5)
2711(6, 6)
280(7, 0)
291(7, 1)
Cache
Value
stored
(rem,
top)
304(7, 2)
318(7, 3)
3211(7, 4)
3313(7, 5)
3414(7, 6)
3515(7, 7)
360(8, 0)
371(8, 1)
385(8, 2)
3910(8, 3)
4015(8, 4)
4118(8, 5)
4220(8, 6)
4321(8, 7)
4422(8, 8)