Problem 181: Investigating in How Many Ways Objects of Two Different Colours Can Be Grouped

Project Euler Problem 181 Solution

EulerSolve provides an optimized solution for Project Euler Problem 181, Investigating in How Many Ways Objects of Two Different Colours Can Be Grouped, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We have $B$ black objects and $W$ white objects. A single group is determined only by its colour composition $(b,w)$ with $b,w \ge 0$ and $b+w>0$. A complete grouping is therefore a multiset of such compositions whose total number of black objects is $B$ and whose total number of white objects is $W$. Reordering the groups does not create a new solution. For Project Euler, the target instance is $60$ black objects and $40$ white objects, but the derivation used by the implementations works for general $B$ and $W$. The central question is not which object goes into which group, but how many times each possible colour composition occurs. Mathematical Approach The cleanest formulation is to treat every possible group composition as a type and then count how many copies of each type are used. Encoding a grouping by multiplicities Let $$\mathcal{T}=\{(t_b,t_w)\mid 0 \le t_b \le B,\ 0 \le t_w \le W,\ (t_b,t_w)\ne(0,0)\}.$$ For each type $t=(t_b,t_w)\in\mathcal{T}$, let $k_t$ be the number of groups having exactly $t_b$ black objects and $t_w$ white objects....

Detailed mathematical approach

Problem Summary

We have $60$ black objects and $40$ white objects and want to split all of them into nonempty groups. As the example of the statement shows, objects of the same colour cannot be told apart, and neither the order of the groups nor the order inside a group matters: three black objects $\mathtt B$ and one white object $\mathtt W$ can be grouped in $7$ ways, $\mathtt{(BBBW)}$, $\mathtt{(B,BBW)}$, $\mathtt{(B,B,BW)}$, $\mathtt{(B,B,B,W)}$, $\mathtt{(B,BB,W)}$, $\mathtt{(BBB,W)}$ and $\mathtt{(BB,BW)}$. We must find in how many ways sixty black and forty white objects can be grouped in the same way. ^[1]

The whole solution in one sentence: a grouping is just an unordered list of group types $(b,w)$, meaning $b$ black and $w$ white objects in a group, whose black and white parts add up to $60$ and $40$, and we count such lists with a table of $61\cdot41$ numbers into which the $2500$ possible types are brought one at a time.

Mathematical Approach

We describe groupings by pairs of numbers, check small cases by hand and then count with a table.

1. A group is a pair of numbers

Because objects of one colour cannot be told apart, a group is described completely by its type $(b,w)$: it holds $b$ black and $w$ white objects, where $b,w\ge0$ and $(b,w)\ne(0,0)$. A grouping is a collection of types in which the order does not matter but a type may occur several times, a multiset; the number of occurrences of a type is its multiplicity. ^[2] For example, $\mathtt{(B,B,BW)}$ is the multiset $\{(1,0),(1,0),(1,1)\}$, in which $(1,0)$ has multiplicity $2$. The seven groupings of the statement, in the order listed there, are $\{(3,1)\}$, $\{(1,0),(2,1)\}$, $\{(1,0),(1,0),(1,1)\}$, $\{(1,0),(1,0),(1,0),(0,1)\}$, $\{(1,0),(2,0),(0,1)\}$, $\{(3,0),(0,1)\}$ and $\{(2,0),(1,1)\}$.

Conversely, every multiset of types whose black parts add up to $B$ and whose white parts add up to $W$ describes exactly one grouping. Let $p(B,W)$ be the number of these multisets, with $p(0,0)=1$ for the empty one; we need $p(60,40)$. They are known as the partitions of the bipartite number $(B,W)$. ^[5]

2. Small cases by hand

With black objects only, a grouping of $n$ objects writes $n$ as a sum of positive group sizes without regard to order, a partition of $n$, ^[3] so $p(n,0)$ is the partition number $p(n)$: $1,1,2,3,5,7,11$ for $n=0,1,\dots,6$. Exchanging the two colours shows that $p(B,W)=p(W,B)$.

With one white object, its group holds some number $k$ of black objects, $0\le k\le n$, and the other $n-k$ black objects form any grouping of their own. Hence

$$p(n,1)=\sum_{k=0}^{n}p(n-k),$$

the partial sums $1,2,4,7,12,19,30$ of the partition numbers. ^[4] For $n=3$ this is $3+2+1+1=7$, the count of the statement: $k=0$ gives $\mathtt{(B,B,B,W)}$, $\mathtt{(B,BB,W)}$ and $\mathtt{(BBB,W)}$; $k=1$ gives $\mathtt{(B,B,BW)}$ and $\mathtt{(BB,BW)}$; $k=2$ gives $\mathtt{(B,BBW)}$; and $k=3$ gives $\mathtt{(BBBW)}$.

The nine groupings $\mathtt{(BBWW)}$, $\mathtt{(B,BWW)}$, $\mathtt{(W,BBW)}$, $\mathtt{(BB,WW)}$, $\mathtt{(BW,BW)}$, $\mathtt{(B,B,WW)}$, $\mathtt{(W,W,BB)}$, $\mathtt{(B,W,BW)}$ and $\mathtt{(B,B,W,W)}$ show that $p(2,2)=9$. Larger cases are tedious to list by hand; the method below gives, for example, $p(3,2)=16$, $p(3,3)=31$, $p(5,2)=47$ and $p(4,3)=57$, and the diagonal values $1,2,9,31,109,339,1043$ of $p(n,n)$ for $n=0,1,\dots,6$. ^[6]

3. The generating function

Let $k_{ij}\ge0$ be the number of groups of type $(i,j)$. A grouping is the same as a choice of all these multiplicities with

$$\sum_{(i,j)\ne(0,0)}k_{ij}\,(i,j)=(B,W).$$

A generating function stores a whole table of numbers as the coefficients of one expression; $[x^by^w]$ denotes the coefficient of $x^by^w$. ^[7] For example, the type $(1,1)$ gets the factor $1+xy+x^2y^2+\dots$, whose terms stand for zero, one, two, … groups $\mathtt{BW}$; multiplying terms adds their exponents, just as putting groups together adds up their objects. In general, each type $(i,j)$ gets the factor

$$1+x^iy^j+x^{2i}y^{2j}+\dots=\frac1{1-x^iy^j},$$

whose term $x^{ki}y^{kj}$ stands for $k$ groups of this type. Multiplying out the product of all factors picks one term from every factor, all but finitely many of them equal to $1$, that is, one multiplicity for every type, and the picked terms multiply to $x^by^w$, where $(b,w)$ is the total of the chosen groups. Hence

$$\begin{aligned}F(x,y)&=\sum_{b,w\ge0}p(b,w)\,x^by^w\\&=\prod_{(i,j)\ne(0,0)}\frac1{1-x^iy^j},\end{aligned}$$

and we need $[x^{60}y^{40}]F$. Setting $y=0$ keeps only the factors with $j=0$ and gives Euler's product for the partition numbers, ^[3] and the classical money-changing example works in the same way, with one factor for each kind of coin. ^[7] A factor with $i\gt B$ or $j\gt W$ contributes only its constant term $1$ to $[x^By^W]F$, since its other terms have too large an exponent; so only the $T=(B+1)(W+1)-1$ types inside the box matter, $2500$ of them for $(60,40)$, and terms with an exponent above $B$ in $x$ or above $W$ in $y$ can be dropped.

4. Bringing in one type at a time

Put the types in a fixed order $t_1,t_2,\dots,t_T$, and let $N_s(b,w)$ be the number of multisets that use only the types $t_1,\dots,t_s$ and add up to $(b,w)$. Without any type only the empty multiset exists, so $N_0(0,0)=1$ and $N_0(b,w)=0$ otherwise. If $t_s=(i,j)$, then

$$N_s(b,w)=\sum_{k\ge0}N_{s-1}(b-ki,w-kj),$$

where entries with a negative index count as $0$. To see this, sort the multisets counted on the left by the number $k$ of copies of $t_s$ they contain. Removing those copies leaves a multiset of the earlier types with the total $(b-ki,w-kj)$; adding $k$ copies to such a multiset reverses this, and different values of $k$ give different multisets. So both sides count the same multisets, whatever the order of the types, and after the last type $N_T(B,W)=p(B,W)$. In the language of Section 3, step $s$ multiplies the product of the first $s-1$ factors by the factor of $t_s$. Building each table from the previous one, so that every partial count is computed only once, is dynamic programming. ^[8]

5. Worked example: three black objects and one white object

For $B=3$ and $W=1$ the program takes the seven types in the order $(0,1)$, $(1,0)$, $(1,1)$, $(2,0)$, $(2,1)$, $(3,0)$, $(3,1)$. We write the table as two rows, the row $w=0$ and the row $w=1$, each listing the entries for $b=0,1,2,3$. It starts as $1,0,0,0$ and $0,0,0,0$.

$(0,1)$: a lone white object can be added to the empty multiset, so the row $w=1$ becomes $1,0,0,0$.
$(1,0)$: any number of lone black objects can be added, so every entry becomes the sum of itself and all entries to its left in its row: $1,1,1,1$ and $1,1,1,1$.
$(1,1)$: $N(b,1)$ gains $N(b-1,0)$, and the row $w=1$ becomes $1,2,2,2$.
$(2,0)$: $N(b,w)$ gains $N(b-2,w)$, giving $1,1,2,2$ and $1,2,3,4$.
$(2,1)$: $N(b,1)$ gains $N(b-2,0)$, giving $1,2,4,5$.
$(3,0)$: $N(3,w)$ gains $N(0,w)$, giving $1,1,2,3$ and $1,2,4,6$.
$(3,1)$: $N(3,1)$ gains $N(0,0)=1$, giving $1,2,4,7$.

The final rows are $p(b,0)=1,1,2,3$ and $p(b,1)=1,2,4,7$ for $b=0,1,2,3$, as in Section 2. The entry for $(3,1)$ grows as $0,1,2,4,5,6,7$, and each increase is a grouping that uses the new type: $\mathtt{(B,B,B,W)}$ after $(1,0)$, $\mathtt{(B,B,BW)}$ after $(1,1)$, $\mathtt{(B,BB,W)}$ and $\mathtt{(BB,BW)}$ after $(2,0)$, $\mathtt{(B,BBW)}$ after $(2,1)$, $\mathtt{(BBB,W)}$ after $(3,0)$ and $\mathtt{(BBBW)}$ after $(3,1)$.

6. One table, updated from the top down

The program overwrites a single table. For the type $(i,j)$ it walks through the cells $(b,w)$ with $b$ from $B$ down to $0$ and, for each $b$, $w$ from $W$ down to $0$, and adds the value of the current cell, the base, to the cells $(b+ki,w+kj)$, $k=1,2,\dots$, that lie in the table. Every target is visited before its base: it has a larger $b$ if $i\ge1$, and the same $b$ with a larger $w$ if $i=0$. So a cell receives additions only from bases visited after it, and when it is read as a base it still holds its old value from $N_{s-1}$. Each cell thus receives $N_{s-1}(b-ki,w-kj)$ exactly once for every $k\ge1$, on top of its own old value, and ends with $N_s(b,w)$ as in Section 4.

The direction matters. Take the single type $(1,0)$ and the cells $b=0,1,2$ of the row $w=0$, holding $1,0,0$; the correct result is $1,1,1$. Walking upwards, the base $b=0$ adds $1$ to $b=1$ and to $b=2$, and then the base $b=1$, which already holds its new value $1$, adds it to $b=2$ once more. The result $1,1,2$ counts the grouping $\mathtt{(B,B)}$ twice: once as two copies taken together and once as one copy added to another.

7. A shorter update

Removing one copy instead of all copies gives a recurrence with a single term: a multiset counted by $N_s(b,w)$ either contains no copy of $t_s=(i,j)$, or removing one copy leaves a multiset counted by $N_s(b-i,w-j)$, and adding one copy reverses this. Hence

$$\begin{aligned}N_s(b,w)={}&N_{s-1}(b,w)\\&+N_s(b-i,\,w-j).\end{aligned}$$

Here the walk must go upwards, because the cell $(b-i,w-j)$ has to hold its new value already. A type $(i,j)$ needs one addition for each of the $(B-i+1)(W-j+1)$ cells with $b\ge i$ and $w\ge j$. Over $i=0,1,\dots,B$ the first factor sums to $1+2+\dots+(B+1)=\binom{B+2}{2}$, the second likewise to $\binom{W+2}{2}$, and their product also counts the excluded type $(0,0)$, whose $(B+1)(W+1)$ cells we subtract. So the whole computation needs

$$\begin{aligned}&\binom{B+2}{2}\binom{W+2}{2}\\&\quad-(B+1)(W+1)\end{aligned}$$

additions, $1\,625\,650$ for $(60,40)$, against $2\,508\,743$ in the loops of the C++ program. Both give the same table in a moment.

8. How large the numbers get

Every number in the table is at most $p(60,40)$. Once a type has been processed, a cell holds $N_s(b,w)\le p(b,w)$, since it counts only some of the groupings of $b$ black and $w$ white objects; while the type is being processed, the cell holds $N_{s-1}(b,w)$ plus some of the nonnegative terms of Section 4, so it never exceeds $N_s(b,w)$. Moreover $p(b,w)\le p(b+1,w)$, because adding a group of one black object to a grouping is a one-to-one map: it is undone by removing one such group. In the same way $p(b,w)\le p(b,w+1)$. An exact computation with unbounded integers shows $p(60,40)\lt2^{63}$, so the signed $64$-bit integers of the C++ and Java programs cannot overflow. The numbers grow fast, however: every $p(b,w)$ with $b+w\le114$ is below $2^{63}$, but $p(60,55)$ and $p(57,58)$ are not, so the C++ options described below give exact results for every $B+W\le114$.

How the Code Works

The C++ program ^[9] starts in main, which reads the options with parse_arguments: the defaults in the struct Options are $60$ and $40$, --black= and --white= with decimal digits change them, --skip-checkpoints switches the checks off, and any other argument gives exit code $1$. The digits are accumulated in an int without a range check, so a value beyond its range is undefined behaviour. ^[10] Unless switched off, run_checkpoints then requires solve(3, 1) to return $7$ and solve to agree with brute_solve for $(4,3)$ and $(5,2)$, where both give $57$ and $47$, with exit code $2$ otherwise; finally main prints the result of solve.

build_group_types lists the $T$ types, the number of black objects in the outer loop and of white objects in the inner loop, without $(0,0)$: the order of Section 5. solve keeps the table of Section 6 in dp, $B+1$ rows of $W+1$ numbers of type i64, an alias of long long, with dp[0][0] set to $1$. For each type (tb, tw) it walks downwards over b and w and skips bases equal to $0$, which would add nothing. It computes max_repeat, the largest $k$ whose target stays in the table, separately for tb == 0, for tw == 0 and for both positive, so that it never divides by zero, and adds the base for $k=1,2,\dots$ up to max_repeat.

brute_dfs works from the other end: for the next type it tries every possible number of copies and recurses on the remaining objects, returning $1$ when all types have been used and nothing remains. The results for the states (idx, rem_b, rem_w) are kept in an unordered_map, so no state is expanded twice: the top-down form of dynamic programming called memoization. ^[8] It shares the type list with solve but not the in-place walk, so the checkpoints test Section 6 against Section 4.

The Python and Java programs contain only the table computation. Both fix $40$ black and $60$ white objects, as the header comment of the Python file also says, which exchanges the colours of the statement. Since $p(40,60)=p(60,40)$ by Section 2, they print the same number, and their final tables are the transposes of the C++ table. Python builds the type list with a list comprehension and has unbounded integers; Java runs through the types in two nested loops and uses a two-dimensional array of signed $64$-bit long values, which suffice by Section 8.

Complexity and Verification

There are $T=(B+1)(W+1)-1$ types, and for each of them the walk visits all $(B+1)(W+1)$ cells: $2500\cdot2501=6\,252\,500$ visits for $(60,40)$. A base can take $k$ copies of $(i,j)$ only if $b+ki\le B$ and $w+kj\le W$, so most types fit only a few times above most cells, and the loops over $k$ add only $2\,508\,743$ times in C++. In general, for a fixed $k$ at most $(B/k)(W/k)$ two-colour types fit $k$ times, each above at most $(B+1)(W+1)$ bases, and the sum of all $1/k^2$ is below $2$, as $1/k^2\le1/(k-1)-1/k$ for $k\ge2$ and these differences add up to $1$; so these types cause fewer than $2BW(B+1)(W+1)$ additions. The at most $(B+W)/k$ one-colour types that fit $k$ times add $O(n^3\log n)$ when $B$ and $W$ are proportional to $n$. So the method needs $\Theta(n^4)$ cell visits and additions and $\Theta(n^2)$ memory: one table of $2501$ numbers plus, in C++ and Python, the list of $2500$ types, about $40$ kilobytes in C++. On a desktop computer a whole run, start-up and the C++ checkpoints included, takes about $0.01$ seconds in C++, $0.04$ seconds in Java and one second in Python.

All three programs were run for this article and print the same number. A brute force that follows the statement literally gives exactly $p(b,w)$ for every $b+w\le12$: for every $n\le12$ it runs through all partitions of $n$ numbered objects into groups, $5\,034\,585$ partitions in all, colours the first $b$ objects black for every $b$ and counts the distinct multisets of group types. With exact integers, the shorter update of Section 7, the transposed table for $40$ black and $60$ white objects and the recurrence

$$\begin{aligned}b\,p(b,w)&=\sum_{u=1}^{b}\sum_{v=0}^{w}c(u,v)\\&\qquad\cdot p(b-u,\,w-v)\end{aligned}$$

all reproduce the whole $61\times41$ table; the recurrence is used row by row for $b\ge1$, starting from the row $p(0,w)$ that it gives with the colours exchanged. Here $c(u,v)$ is the sum of $u/m$ over the positive common divisors $m$ of $u$ and $v$, for $v=0$ all divisors of $u$; the recurrence follows by comparing the coefficients of $x\,\partial F/\partial x$ and $F\cdot x\,\partial(\ln F)/\partial x$, since $\ln F$ is the sum of $x^{mi}y^{mj}/m$ over all types $(i,j)$ and all $m\ge1$. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The reduction to multisets of types, the in-place update and its proof, the operation counts and the size bounds are derived in this article.

Project Euler 181 — Grouping Two Different Coloured Objects. The official statement lists the 7 ways to group three black objects B and one white object W and asks in how many ways sixty black and forty white objects can be grouped.
Wikipedia — Multiset. A multiset is a set in which an element may occur several times, its multiplicity; the order of listing does not matter.
Wikipedia — Partition (number theory). A partition of n writes n as a sum of positive integers without regard to order; p(0), …, p(10) are 1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, and the generating function of p(n) is the product of the factors 1/(1 − qʲ) over all j ≥ 1.
OEIS A000070. The partial sums p(0) + p(1) + … + p(n) of the partition numbers: 1, 2, 4, 7, 12, 19, 30, …
OEIS A054225. The bipartite partitions, that is, the partitions of the pairs (n − k, k) into sums of pairs, with their generating function, a product of the factors 1/(1 − x^(i−j) y^j). The entry cites MacMahon (1890) and Auluck (1953).
OEIS A002774. The diagonal p(n, n), the bipartite partitions of n white and n black objects: 1, 2, 9, 31, 109, 339, 1043, …
Wikipedia — Generating function. A generating function encodes a sequence, or in several variables an array, as the coefficients of a formal power series, [xⁿ] denoting the coefficient of xⁿ; in its money-changing example one factor for each kind of coin counts the ways to pay n cents, ignoring the order of the coins.
Wikipedia — Dynamic programming. Dynamic programming breaks a problem into simpler sub-problems and solves each of them only once, bottom-up in a table or top-down by storing the results of a recursion, which is called memoization.
C++ — Euler181.cpp. The linked immutable C++ revision contains build_group_types, solve, brute_dfs and run_checkpoints. It is the implementation record for this article.
cppreference — Arithmetic operators. The overflow of a signed integer operation is undefined behaviour, while unsigned integer arithmetic is performed modulo 2ⁿ.

Problem 181 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using i64 = long long;

struct Options {
    int black = 60;
    int white = 40;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }

    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--black=", options.black) ||
            parse_int_after_prefix(arg, "--white=", options.white)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.black >= 0 && options.white >= 0;
}

std::vector<std::pair<int, int>> build_group_types(const int black, const int white) {
    std::vector<std::pair<int, int>> types;
    types.reserve(static_cast<std::size_t>((black + 1) * (white + 1) - 1));
    for (int b = 0; b <= black; ++b) {
        for (int w = 0; w <= white; ++w) {
            if (b == 0 && w == 0) {
                continue;
            }
            types.emplace_back(b, w);
        }
    }
    return types;
}

i64 solve(const int black, const int white) {
    const std::vector<std::pair<int, int>> types = build_group_types(black, white);

    std::vector<std::vector<i64>> dp(
        static_cast<std::size_t>(black + 1),
        std::vector<i64>(static_cast<std::size_t>(white + 1), 0));
    dp[0][0] = 1;

    for (const auto& [tb, tw] : types) {
        for (int b = black; b >= 0; --b) {
            for (int w = white; w >= 0; --w) {
                const i64 base = dp[static_cast<std::size_t>(b)][static_cast<std::size_t>(w)];
                if (base == 0) {
                    continue;
                }

                int max_repeat = 0;
                if (tb == 0) {
                    max_repeat = (tw == 0) ? 0 : (white - w) / tw;
                } else if (tw == 0) {
                    max_repeat = (black - b) / tb;
                } else {
                    max_repeat = std::min((black - b) / tb, (white - w) / tw);
                }

                for (int k = 1; k <= max_repeat; ++k) {
                    dp[static_cast<std::size_t>(b + k * tb)][static_cast<std::size_t>(w + k * tw)] += base;
                }
            }
        }
    }

    return dp[static_cast<std::size_t>(black)][static_cast<std::size_t>(white)];
}

i64 brute_dfs(const int idx, const int rem_b, const int rem_w,
              const std::vector<std::pair<int, int>>& types,
              std::unordered_map<std::uint64_t, i64>& memo) {
    if (rem_b < 0 || rem_w < 0) {
        return 0;
    }
    if (idx == static_cast<int>(types.size())) {
        return (rem_b == 0 && rem_w == 0) ? 1 : 0;
    }

    const std::uint64_t key = (static_cast<std::uint64_t>(idx) << 32) |
                              (static_cast<std::uint64_t>(rem_b) << 16) |
                              static_cast<std::uint64_t>(rem_w);
    const auto it = memo.find(key);
    if (it != memo.end()) {
        return it->second;
    }

    const int tb = types[static_cast<std::size_t>(idx)].first;
    const int tw = types[static_cast<std::size_t>(idx)].second;

    int max_repeat = 0;
    if (tb == 0) {
        max_repeat = rem_w / tw;
    } else if (tw == 0) {
        max_repeat = rem_b / tb;
    } else {
        max_repeat = std::min(rem_b / tb, rem_w / tw);
    }

    i64 ways = 0;
    for (int k = 0; k <= max_repeat; ++k) {
        ways += brute_dfs(
            idx + 1,
            rem_b - k * tb,
            rem_w - k * tw,
            types,
            memo);
    }

    memo.emplace(key, ways);
    return ways;
}

i64 brute_solve(const int black, const int white) {
    const std::vector<std::pair<int, int>> types = build_group_types(black, white);
    std::unordered_map<std::uint64_t, i64> memo;
    return brute_dfs(0, black, white, types, memo);
}

bool run_checkpoints() {
    if (solve(3, 1) != 7LL) {
        std::cerr << "Checkpoint failed for 3 black and 1 white" << '\n';
        return false;
    }
    if (solve(4, 3) != brute_solve(4, 3)) {
        std::cerr << "Checkpoint failed for brute cross-check at (4,3)" << '\n';
        return false;
    }
    if (solve(5, 2) != brute_solve(5, 2)) {
        std::cerr << "Checkpoint failed for brute cross-check at (5,2)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.black, options.white) << '\n';
    return 0;
}

Python

# Problem 181: Investigating in how many ways objects of two different colours can be grouped.
# Partition (40 black, 60 white) into groups of (b,w) objects.

def solve():
    black, white = 40, 60
    # Group types: all (b,w) with b+w > 0
    types = [(b, w) for b in range(black+1) for w in range(white+1) if b+w > 0]
    dp = [[0]*(white+1) for _ in range(black+1)]
    dp[0][0] = 1
    for tb, tw in types:
        for b in range(black, -1, -1):
            for w in range(white, -1, -1):
                base = dp[b][w]
                if base == 0: continue
                if tb == 0: mr = (white-w)//tw if tw > 0 else 0
                elif tw == 0: mr = (black-b)//tb
                else: mr = min((black-b)//tb, (white-w)//tw)
                for k in range(1, mr+1):
                    dp[b+k*tb][w+k*tw] += base
    print(dp[black][white])

solve()

Java

public class Euler181 {
    public static void main(String[] args) {
        int B = 40, W = 60;
        long[][] dp = new long[B + 1][W + 1];
        dp[0][0] = 1;
        for (int tb = 0; tb <= B; tb++)
            for (int tw = 0; tw <= W; tw++) {
                if (tb == 0 && tw == 0)
                    continue;
                for (int b = B; b >= 0; b--)
                    for (int w = W; w >= 0; w--) {
                        long base = dp[b][w];
                        if (base == 0)
                            continue;
                        int mr;
                        if (tb == 0)
                            mr = tw == 0 ? 0 : (W - w) / tw;
                        else if (tw == 0)
                            mr = (B - b) / tb;
                        else
                            mr = Math.min((B - b) / tb, (W - w) / tw);
                        for (int k = 1; k <= mr; k++)
                            dp[b + k * tb][w + k * tw] += base;
                    }
            }
        System.out.println(dp[B][W]);
    }
}