Problem 194: Coloured Configurations

Project Euler Problem 194 Solution

EulerSolve provides an optimized solution for Project Euler Problem 194, Coloured Configurations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The configuration is a chain of $a+b$ colored units. Exactly $a$ of them are type A and $b$ are type B. Adjacent units share an ordered boundary pair of vertices, and every edge in the unit graph imposes the usual coloring rule that its endpoints must have different colors. For Project Euler 194 the parameters are $a=25$, $b=75$, and $c=1984$, and the final answer is required modulo $10^8$. A direct count of whole-chain colorings would mix together 100 units and an enormous number of color assignments. The implementations avoid that by isolating one local counting problem: how many valid extensions does a single unit have when the two colors on its incoming boundary are already fixed and distinct? Mathematical Approach Fix an ordered incoming boundary pair $(x,y)$ with $x \ne y$. A unit contributes five unfixed vertices: two outgoing boundary vertices and three internal vertices. Because the constraints only say that adjacent vertices must have different colors, the actual names of $x$ and $y$ never matter; only the fact that they are distinct matters. The Boundary-Pair Invariant Both unit types start from a distinct ordered pair $(x,y)$ and end at another distinct ordered pair. The outgoing top boundary color must differ from the incoming top one, and the outgoing bottom boundary color must differ from the outgoing top one....

Detailed mathematical approach

Problem Summary

Each of the two units $A$ and $B$ is a small graph with a vertex at each corner of a rectangle and three more vertices on a vertical line through its middle: the upper inner vertex is joined to both top corners, the lower one to both bottom corners, and the middle one to the other two. Unit $A$ also has the four sides of the rectangle as edges; unit $B$ lacks the bottom side. The units stand in a row, glued along their vertical sides: the right side of each unit is the left side of the next. A configuration of type $(a,b,c)$ is such a row of $a$ units $A$ and $b$ units $B$, in any order, whose vertices are coloured from a set of $c$ colours so that the two ends of every edge differ; not every colour has to be used. $N(a,b,c)$ is the number of these configurations; the figure of the statement shows the row $A,B,B,A$ coloured with four colours. We are given $N(1,0,3)=24$, $N(0,2,4)=92928$ and $N(2,2,3)=20736$ and must find the last $8$ digits of $N(25,75,1984)$. ^[1]

The whole solution in one sentence: neighbouring units share one edge, whose ends always have different colours, and renaming the colours shows that the number of ways to finish a unit does not depend on those two colours, so $N(a,b,c)$ is the number $\binom{a+b}{a}$ of orders of the units times $c(c-1)\,A(c)^a\,B(c)^b$, where $A(c)$ and $B(c)$, the numbers of ways to finish one unit of each kind, are polynomials of degree $5$ that we find by counting.

Mathematical Approach

1. Names for the vertices of a unit

Call the corners of a unit $t_0$ (top left), $b_0$ (bottom left), $t_1$ (top right) and $b_1$ (bottom right), and its inner vertices $u$ (upper), $m$ (middle) and $w$ (lower). Unit $B$ has nine edges, the sides $t_0b_0$, $t_0t_1$, $t_1b_1$ and the edges $ut_0$, $ut_1$, $um$, $mw$, $wb_0$, $wb_1$; unit $A$ also has the bottom side $b_0b_1$. A row of $k$ units has $2+5k$ vertices; the figure of the statement, with $k=4$, has $22$ vertices and $35$ edges. We write $x,y$ for the colours of $t_0,b_0$ and $p,q$ for those of $t_1,b_1$. Let $A(c)$ be the number of ways to colour $t_1,b_1,u,m,w$ in a unit $A$ whose left side already has two colours $x\ne y$, and $B(c)$ the same number for a unit $B$; by Section 3 they do not depend on $x$ and $y$.

2. One unit with three colours

Take the colours $1,2,3$ and fix $x=1$, $y=2$. With three colours, a vertex joined to two vertices of different colours has exactly one colour left. The corner $t_1$ is joined to $t_0$, so $p$ is $2$ or $3$; the corner $b_1$ is joined to $t_1$ and, in unit $A$, to $b_0$, so $q\ne p$ and $q\ne2$. Unit $A$ thus allows three pairs $(p,q)$:

$(p,q)=(2,1)$: $u$ avoids $1,2$ and $w$ avoids $2,1$, so $u=w=3$, and $m$ is $1$ or $2$: two colourings;
$(p,q)=(2,3)$: $u=3$, $w=1$ and $m=2$: one colouring;
$(p,q)=(3,1)$: $u=2$, $w=3$ and $m=1$: one colouring.

So $A(3)=4$. Unit $B$, without the edge $b_0b_1$, also allows $(p,q)=(3,2)$: then $u=2$, $w$ is $1$ or $3$ and $m$ takes the colour left over, two more colourings, so $B(3)=6$. Each of the $3\cdot2=6$ colourings of the left side leaves $4$ ways by Section 3, so one unit $A$ has $24$ colourings: the given $N(1,0,3)=24$.

3. Renaming colours

Claim: if $x\ne y$ and $x'\ne y'$, a unit has as many colourings in which $t_0,b_0$ have the colours $x,y$ as colourings in which they have the colours $x',y'$. Proof: take a permutation $\sigma$ of the $c$ colours with $\sigma(x)=x'$ and $\sigma(y)=y'$; it exists because each pair consists of two different colours. Applying $\sigma$ to the colour of every vertex keeps the ends of every edge different, since $\sigma$ is one-to-one, and turns the colourings that start with $x,y$ into those that start with $x',y'$; $\sigma^{-1}$ undoes this, so both sets have the same size.

Hence $A(c)$ and $B(c)$ are well defined, and one unit has $c(c-1)A(c)$ or $c(c-1)B(c)$ colourings. These are values of the chromatic polynomials of the units, which count proper colourings, with different colours at the ends of every edge, from $c$ distinguishable colours. ^[2]

4. Gluing the units

Fix an order of the $a+b$ units and colour the row from left to right. The leftmost side can be coloured in $c(c-1)$ ways. When a unit is reached, its left side is the right side of the previous unit and, being an edge, already has two different colours. Among the vertices coloured so far, the five new vertices of the unit are joined only to that side, so by Section 3 there are $A(c)$ or $B(c)$ ways to continue, whatever was chosen before. By the rule of product ^[3] each order gives $c(c-1)A(c)^aB(c)^b$ colourings, and an order is fixed by choosing which $a$ of the $a+b$ places hold a unit $A$, in $\binom{a+b}{a}$ ways: ^[4]

$$\begin{aligned}N(a,b,c)&=\binom{a+b}{a}\,c(c-1)\\&\quad\cdot A(c)^a\,B(c)^b.\end{aligned}$$

5. The inner path of a unit

To find $A(c)$ and $B(c)$, fix all four corners and colour the path $u,m,w$. Let $U$ be the set of colours allowed for $u$, all except $x$ and $p$, and $W$ the set allowed for $w$, all except $y$ and $q$. For each colour of $u$ in $U$ and each colour of $w$ in $W$, the vertex $m$ must avoid both: it has $c-2$ choices when they differ and $c-1$ when they are equal, which happens once for every colour in $U\cap W$. So the path has

$$(c-2)\,|U|\,|W|+|U\cap W|$$

colourings, where $|U|=c-2$ because $x\ne p$, $|W|$ is $c-2$ if $q\ne y$ and $c-1$ if $q=y$, and $|U\cap W|=c-|\{x,y,p,q\}|$.

6. Adding over the right corners

The corner colours satisfy $p\ne x$ and $q\ne p$, and unit $A$ also needs $q\ne y$. Call a colour new if it is neither $x$ nor $y$. Every allowed pair $(p,q)$ falls into exactly one of five cases, with the path counts of Section 5:

case 1, $p=y$ and $q=x$: one pair, with $(c-2)^3+(c-2)$ colourings;
case 2, $p=y$ and $q$ new: $c-2$ pairs, with $(c-2)^3+(c-3)$ each;
case 3, $p$ new and $q=x$: $c-2$ pairs, with $(c-2)^3+(c-3)$ each;
case 4, $p$ and $q$ new, $q\ne p$: $(c-2)(c-3)$ pairs, with $(c-2)^3+(c-4)$ each;
case 5, $p$ new and $q=y$, only in unit $B$: $c-2$ pairs, with $(c-1)(c-2)^2+(c-3)$ each.

Cases 1 to 4 make up unit $A$, with $1+2(c-2)+(c-2)(c-3)$${}=c^2-3c+3$ pairs, so

$$\begin{aligned}A(c)&=(c-2)^3(c^2-3c+3)\\&\quad+(c-2)(c^2-5c+7).\end{aligned}$$

Case 5 adds $(c-2)^3(c-1)+(c-2)(c-3)$ colourings. With the terms of $A(c)$, the parts with the factor $(c-2)^3$ add up to $(c-2)^3(c^2-2c+2)$ and the others to $(c-2)(c^2-4c+4)=(c-2)^3$, so

$$B(c)=(c-2)^3(c^2-2c+3).$$

Multiplied out,

$$\begin{aligned}A(c)&=c^5-9c^4+34c^3\\&\quad-69c^2+77c-38,\\B(c)&=c^5-8c^4+27c^3\\&\quad-50c^2+52c-24.\end{aligned}$$

Both formulas hold for every $c\ge2$; for $c\le3$ some cases are simply empty. Case 5 counts the colourings of unit $B$ in which $b_0$ and $b_1$ share a colour, those that the bottom side of unit $A$ forbids; this is an instance of the deletion–contraction rule for chromatic polynomials. ^[2]

7. Worked example: four colours, and the given values

For $c=4$ the cases 1 to 4 contribute $1\cdot10$, $2\cdot9$, $2\cdot9$ and $2\cdot8$ colourings, so $A(4)=62$ from $7$ pairs $(p,q)$. Case 5 adds $2\cdot13=26$, so $B(4)=88$ from $9$ pairs. The closed forms agree: $A(4)=8\cdot7+2\cdot3$ and $B(4)=8\cdot11$. For $c=3$ the cases give $2+1+1+0=4$, and case 5 adds $2$, the colourings found by hand in Section 2. The formula of Section 4 now reproduces the three given values, where the first factor $6$ of the last line is $\binom42$:

$$\begin{aligned}N(1,0,3)&=6\cdot4=24,\\N(0,2,4)&=12\cdot88^2=92928,\\N(2,2,3)&=6\cdot6\cdot4^2\cdot6^2=20736.\end{aligned}$$

These values also settle two points of reading: all $6$ orders of two units $A$ and two units $B$ count separately, even mirror images such as $A,A,B,B$ and $B,B,A,A$; and units are never turned upside down, because allowing a turned unit $B$, which lacks its top side instead, would multiply $N(0,2,4)$ by $4$.

8. Six values fix each polynomial

The published programs do not use these closed forms. They start from the brute-force values for $c=2,3,\dots,7$,

$$\begin{aligned}A&:\ 0,\ 4,\ 62,\ 372,\ 1396,\ 3980,\\B&:\ 0,\ 6,\ 88,\ 486,\ 1728,\ 4750,\end{aligned}$$

and rebuild each polynomial from them. This works because both are polynomials of degree $5$, as Section 6 shows. Two different polynomials of degree at most $5$ cannot agree at six points, since their difference would have six roots, so six values determine each polynomial. ^[5]

9. Forward differences and Newton's formula

Write the values of $A$ in a row, below it their forward differences $\Delta f(c)=f(c+1)-f(c)$, below these the differences of the differences, and so on: ^[6]

$0,\ 4,\ 62,\ 372,\ 1396,\ 3980$,
$4,\ 58,\ 310,\ 1024,\ 2584$,
$54,\ 252,\ 714,\ 1560$,
$198,\ 462,\ 846$,
$264,\ 384$,
$120$.

The first entries of the rows are the $k$-th differences $\Delta^kA(2)=0,4,54,198,264,120$ for $k=0,\dots,5$; for $B$ they are $0,6,76,240,288,120$. The last one is $120=5!$: every polynomial of degree $5$ with leading coefficient $1$ has the constant fifth difference $5!$, and all its higher differences vanish. ^[6] Newton's forward-difference formula rebuilds the polynomial from the first entries,

$$A(c)=\sum_{k=0}^{5}\Delta^kA(2)\binom{c-2}{k},$$

and $B(c)$ in the same way. The reason: with $n=c-2$, Pascal's rule ^[4] gives $\Delta\binom nk=\binom n{k-1}$, so at $n=0$ the $j$-th difference of the right side is its $j$-th weight $\Delta^jA(2)$. Both sides thus share the first column of the difference table, which fixes the whole table by additions and hence the values at $c=2,\dots,7$; by Section 8 they are the same polynomial. ^[6] At $c=10$, for example, the six terms are $0$, $32$, $1512$, $11088$, $18480$ and $6720$, whose sum $37832$ is $A(10)$, as a brute-force count confirms.

10. Only the last eight digits

$N(25,75,1984)$ has $1679$ digits, but we need only its remainder modulo $10^8$, and remainders may be taken after every multiplication, since the remainder of a product depends only on the remainders of the factors. Divisions are different: $12/2=6$, but reducing $12$ modulo $10$ first would give $2/2=1$. So the binomial coefficient $\binom{100}{25}\approx2.43\cdot10^{23}$ and the weights $\binom{c-2}{k}$, which are built by exact divisions, are computed exactly and reduced afterwards. The powers come from square-and-multiply: ^[7] squaring once per binary digit of the exponent gives $B(c)$, $B(c)^2$, $B(c)^4$ and so on, and the powers that belong to a digit $1$ are multiplied into the result. For $b=75=1001011_2$ this takes $7$ squarings and $4$ multiplications instead of $74$.

How the Code Works

The C++ program ^[8] starts in main, which reads the optional arguments --a=, --b=, --c= and --mod= (by default the problem's parameters and the modulus $10^8$), runs run_checkpoints unless --skip-checkpoints is given, and prints the result of count_configurations_mod, padded with leading zeros to eight digits for the modulus $10^8$.

count_unit_extensions_bruteforce fixes the left side as the colours $0$ and $1$, as Section 3 allows, and five nested loops try every colour for right_top, right_bottom, upper, middle and lower, the vertices $t_1,b_1,u,m,w$, skipping every colour of an already coloured neighbour; the test against left_bottom in the loop for right_bottom, the bottom side, is made only for unit $A$. build_forward_differences calls it for $c=2,\dots,7$, replaces the row of six values five times by its differences and keeps the first entry of each row, the numbers of Section 9. evaluate_from_forward_differences evaluates Newton's formula with $n=c-2$, updating $\binom nk=\binom n{k-1}(n-k+1)/k$, a division that is always exact. ^[4] unit_extensions builds the two tables once and keeps them in static variables.

run_checkpoints compares the interpolated values with fresh brute-force counts for every $c$ from $2$ to $10$, where $c=8,9,10$ lie outside the six sample points and so test the interpolation beyond them, and checks the three given values exactly with count_configurations_exact_small. count_configurations_mod reduces $\binom{a+b}{a}$, computed exactly by choose_exact with the same multiplicative rule in 128-bit integers, ^[9] and multiplies it by $c(c-1)$ and by the powers of $A(c)$ and $B(c)$, reduced by normalize_mod, from pow_mod, the method of Section 10; mul_mod forms every product in unsigned 128-bit arithmetic, so it cannot overflow.

For the problem's parameters every 128-bit intermediate value stays below $2^{83}$. The parsers check only that each argument fits its integer type, not that later arithmetic fits: $a=b=64$ overflows a signed 128-bit product in choose_exact, and $c$ above about $4.4\cdot10^7$ overflows the interpolation, which is undefined behaviour in C++ and in our tests gave a wrong result.

The Python and Java versions compute the same formula for the fixed parameters, without arguments or checks, from the six values of each unit written into the source. eval_poly and evalPoly form the first column of differences in place, updating from right to left, and evaluate Newton's formula exactly, with Python's built-in integers or Java's BigInteger. ^[9] Python uses math.comb and the three-argument pow; Java builds the binomial coefficient with BigInteger in a loop of $a$ steps and forms every modular product with it as well. Neither pads the printed remainder with zeros.

Complexity and Verification

The brute force tries at most $c^5$ colourings for each $c\le10$, fewer than half a million in all, and the interpolation is a fixed amount of work. Counting each operation on fixed-width integers as one step, the C++ binomial coefficient takes $\min(a,b)$ iterations and each power one per binary digit of its exponent, so for positive $a$ and $b$ the time is $O(\min(a,b)+\log a+\log b)$ and the memory $O(1)$. Java's binomial loop takes $a$ iterations and Python calls math.comb; at these sizes their exact integers cost next to nothing. All three programs finish in well under a second.

The C++ checkpoints pass, and the three programs print the same eight digits as an exact evaluation of the formula of Section 4 with Python's unbounded integers. Further checks for this article, not part of the published solutions:

a direct count of all colourings of the whole glued graph, without the product formula, reproduces $N(a,b,c)$ in $42$ cases, $c=2,3$ with $1\le a+b\le4$, $c=4$ with $1\le a+b\le3$ and $c=5$ with $1\le a+b\le2$, among them the three given values;
a left-to-right count with the actual colours of every shared side, which avoids the renaming argument, agrees for every order of at most $6$ units and every $c$ from $2$ to $7$;
the closed forms of Section 6 match brute-force values for $c$ from $2$ to $12$, and Newton's formula matches them for every $c$ up to $2000$.

Footnotes and References

The references below supply the background facts used above; the case analysis of Section 6 and the reading of the statement checked in Section 7 are derived in this article.

Project Euler 194 — Coloured Configurations. The official statement glues the units A and B, shown as pictures, along their vertical edges; a configuration of type (a, b, c) is a graph built of a units A and b units B whose vertices are coloured using up to c colours so that no two adjacent vertices have the same colour. It gives N(1, 0, 3) = 24, N(0, 2, 4) = 92928 and N(2, 2, 3) = 20736 and asks for the last 8 digits of N(25, 75, 1984).
Wikipedia — Chromatic polynomial. It counts the proper vertex colourings of a graph with x colours, is a monic polynomial of degree n for a graph with n vertices, and satisfies the deletion–contraction recurrence P(G, x) = P(G − uv, x) − P(G / uv, x).
Wikipedia — Rule of product. With a ways of doing something and b ways of doing another thing, there are a · b ways of doing both.
Wikipedia — Binomial coefficient. C(n, k) counts the k-element subsets of an n-element set. It satisfies Pascal's rule C(n, k) = C(n − 1, k − 1) + C(n − 1, k), the recursive form C(n, k) = C(n, k − 1)(n − k + 1)/k of the multiplicative formula, and C(n, k) = C(n, n − k).
Wikipedia — Newton polynomial. The interpolating polynomial is unique, because the only polynomial of degree at most k with more than k roots is zero.
Wikipedia — Finite difference. The n-th forward difference of a polynomial of degree n with leading coefficient a is the constant a·hⁿ·n!, and further differences are 0; Newton's series f(x) = Σ C(x − x₀, k) Δᵏ[f](x₀) holds for any polynomial f.
Wikipedia — Modular exponentiation. The right-to-left binary method squares the base once per bit of the exponent and multiplies it into the result for every bit 1, reducing modulo m after each product, in O(log exponent) steps.
C++ — Euler194.cpp. The linked immutable C++ revision, with count_unit_extensions_bruteforce, evaluate_from_forward_differences and run_checkpoints, is the implementation record for this article.
GCC manual — 128-bit Integers: __int128 and unsigned __int128 exist on targets with an integer mode wide enough to hold 128 bits. Python — numeric types: integers have unlimited precision. Java — BigInteger: immutable arbitrary-precision integers.

Problem 194 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

struct Options {
    int a = 25;
    int b = 75;
    int c = 1984;
    u64 mod = 100000000ULL;
    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;
    }

    i64 parsed = 0;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<i64>(ch - '0');
        if (parsed > static_cast<i64>(std::numeric_limits<int>::max())) {
            return false;
        }
    }

    value = static_cast<int>(parsed);
    return true;
}

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

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

    u64 parsed = 0ULL;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(ch - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_arguments(const 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, "--a=", options.a)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--b=", options.b)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--c=", options.c)) {
            continue;
        }
        if (parse_u64_after_prefix(arg, "--mod=", options.mod)) {
            continue;
        }

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

    return true;
}

std::string to_string_i128(i128 value) {
    if (value == 0) {
        return "0";
    }

    bool negative = false;
    u128 magnitude = 0;
    if (value < 0) {
        negative = true;
        magnitude = static_cast<u128>(-value);
    } else {
        magnitude = static_cast<u128>(value);
    }

    std::string digits;
    while (magnitude > 0) {
        const int digit = static_cast<int>(magnitude % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        magnitude /= 10U;
    }

    if (negative) {
        digits.push_back('-');
    }
    std::reverse(digits.begin(), digits.end());
    return digits;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) %
                            static_cast<u128>(mod));
}

u64 pow_mod(u64 base, int exp, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }

    u64 result = 1ULL % mod;
    base %= mod;

    while (exp > 0) {
        if ((exp & 1) != 0) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1;
    }

    return result;
}

i128 choose_exact(int n, int k) {
    if (k < 0 || k > n) {
        return 0;
    }
    k = std::min(k, n - k);
    i128 result = 1;
    for (int i = 1; i <= k; ++i) {
        result = (result * static_cast<i128>(n - k + i)) / static_cast<i128>(i);
    }
    return result;
}

i128 pow_i128(i128 base, int exp) {
    i128 result = 1;
    while (exp > 0) {
        if ((exp & 1) != 0) {
            result *= base;
        }
        base *= base;
        exp >>= 1;
    }
    return result;
}

u64 count_unit_extensions_bruteforce(const int c, const bool is_type_a) {
    if (c < 2) {
        return 0ULL;
    }

    constexpr int left_top = 0;
    constexpr int left_bottom = 1;

    u64 count = 0ULL;

    for (int right_top = 0; right_top < c; ++right_top) {
        if (right_top == left_top) {
            continue;
        }

        for (int right_bottom = 0; right_bottom < c; ++right_bottom) {
            if (right_bottom == right_top) {
                continue;
            }
            if (is_type_a && right_bottom == left_bottom) {
                continue;
            }

            for (int upper = 0; upper < c; ++upper) {
                if (upper == left_top || upper == right_top) {
                    continue;
                }

                for (int middle = 0; middle < c; ++middle) {
                    if (middle == upper) {
                        continue;
                    }

                    for (int lower = 0; lower < c; ++lower) {
                        if (lower == middle || lower == left_bottom ||
                            lower == right_bottom) {
                            continue;
                        }
                        ++count;
                    }
                }
            }
        }
    }

    return count;
}

std::array<i128, 6> build_forward_differences(const bool is_type_a) {
    std::array<i128, 6> row{};
    for (int i = 0; i < 6; ++i) {
        row[static_cast<std::size_t>(i)] = static_cast<i128>(
            count_unit_extensions_bruteforce(i + 2, is_type_a));
    }

    std::array<i128, 6> differences{};
    differences[0] = row[0];
    for (int level = 1; level < 6; ++level) {
        for (int i = 0; i < 6 - level; ++i) {
            row[static_cast<std::size_t>(i)] =
                row[static_cast<std::size_t>(i + 1)] -
                row[static_cast<std::size_t>(i)];
        }
        differences[static_cast<std::size_t>(level)] = row[0];
    }

    return differences;
}

i128 evaluate_from_forward_differences(const std::array<i128, 6>& differences,
                                       const int c) {
    if (c < 2) {
        return 0;
    }

    const i64 n = static_cast<i64>(c) - 2LL;
    i128 choose_n_k = 1;
    i128 value = 0;

    for (int k = 0; k < 6; ++k) {
        if (k > 0) {
            choose_n_k = (choose_n_k * static_cast<i128>(n - (k - 1))) /
                         static_cast<i128>(k);
        }
        value += differences[static_cast<std::size_t>(k)] * choose_n_k;
    }

    return value;
}

i128 unit_extensions(const int c, const bool is_type_a) {
    static const std::array<i128, 6> a_differences =
        build_forward_differences(true);
    static const std::array<i128, 6> b_differences =
        build_forward_differences(false);

    return evaluate_from_forward_differences(is_type_a ? a_differences
                                                       : b_differences,
                                             c);
}

u64 normalize_mod(const i128 value, const u64 mod) {
    i128 reduced = value % static_cast<i128>(mod);
    if (reduced < 0) {
        reduced += static_cast<i128>(mod);
    }
    return static_cast<u64>(reduced);
}

i128 count_configurations_exact_small(const int a, const int b, const int c) {
    if (c < 2) {
        return 0;
    }

    const i128 a_extensions = unit_extensions(c, true);
    const i128 b_extensions = unit_extensions(c, false);

    i128 result = choose_exact(a + b, a);
    result *= static_cast<i128>(c) * static_cast<i128>(c - 1);
    result *= pow_i128(a_extensions, a);
    result *= pow_i128(b_extensions, b);
    return result;
}

u64 count_configurations_mod(const int a,
                             const int b,
                             const int c,
                             const u64 mod) {
    if (mod == 0ULL) {
        throw std::runtime_error("Modulus must be positive");
    }
    if (c < 2) {
        return 0ULL;
    }

    const i128 a_extensions = unit_extensions(c, true);
    const i128 b_extensions = unit_extensions(c, false);

    u64 result = static_cast<u64>(choose_exact(a + b, a) % static_cast<i128>(mod));

    const u64 c_mod = static_cast<u64>(static_cast<i128>(c) % static_cast<i128>(mod));
    const u64 c_minus_one_mod =
        static_cast<u64>(static_cast<i128>(c - 1) % static_cast<i128>(mod));
    result = mul_mod(result, mul_mod(c_mod, c_minus_one_mod, mod), mod);
    result = mul_mod(result, pow_mod(normalize_mod(a_extensions, mod), a, mod), mod);
    result = mul_mod(result, pow_mod(normalize_mod(b_extensions, mod), b, mod), mod);

    return result;
}

bool check_equal(const i128 actual,
                 const i128 expected,
                 const std::string& label) {
    if (actual == expected) {
        return true;
    }
    std::cerr << "Checkpoint failed for " << label << ": expected "
              << to_string_i128(expected) << ", got " << to_string_i128(actual)
              << '\n';
    return false;
}

bool run_checkpoints() {
    for (int c = 2; c <= 10; ++c) {
        const i128 a_interp = unit_extensions(c, true);
        const i128 b_interp = unit_extensions(c, false);
        const i128 a_brute =
            static_cast<i128>(count_unit_extensions_bruteforce(c, true));
        const i128 b_brute =
            static_cast<i128>(count_unit_extensions_bruteforce(c, false));

        if (!check_equal(a_interp, a_brute,
                         "unit A interpolation at c=" + std::to_string(c))) {
            return false;
        }
        if (!check_equal(b_interp, b_brute,
                         "unit B interpolation at c=" + std::to_string(c))) {
            return false;
        }
    }

    if (!check_equal(count_configurations_exact_small(1, 0, 3), 24,
                     "N(1,0,3)")) {
        return false;
    }
    if (!check_equal(count_configurations_exact_small(0, 2, 4), 92928,
                     "N(0,2,4)")) {
        return false;
    }
    if (!check_equal(count_configurations_exact_small(2, 2, 3), 20736,
                     "N(2,2,3)")) {
        return false;
    }

    return true;
}

} // namespace

int main(int argc, char** argv) {
    try {
        Options options;
        if (!parse_arguments(argc, argv, options)) {
            return 1;
        }

        if (options.a < 0 || options.b < 0 || options.c < 0) {
            std::cerr << "Parameters --a, --b, and --c must be non-negative.\n";
            return 1;
        }

        if (options.run_checkpoints && !run_checkpoints()) {
            return 1;
        }

        const u64 answer =
            count_configurations_mod(options.a, options.b, options.c, options.mod);

        if (options.mod == 100000000ULL) {
            std::cout << std::setw(8) << std::setfill('0') << answer << '\n';
        } else {
            std::cout << answer << '\n';
        }
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

# Problem 194: Coloured Configurations
# Count colourings of a graph with a=25 type-A nodes, b=75 type-B nodes, c=1984, mod 10^8.
# Formula: C(a+b, a) * c * (c-1) * ext_a(c)^a * ext_b(c)^b mod M

from math import comb

def solve():
    a, b, c = 25, 75, 1984
    M = 10**8

    # Brute force values for c=2..7:
    # type_a: 0, 4, 62, 372, 1396, 3980
    # type_b: 0, 6, 88, 486, 1728, 4750
    # These are degree-5 polynomials. Use Newton forward differences.

    def eval_poly(vals, c_val):
        """Evaluate polynomial at c_val using Newton forward differences from values at 2,3,...,7"""
        diffs = vals[:]
        for level in range(1, len(diffs)):
            for i in range(len(diffs)-1, level-1, -1):
                diffs[i] -= diffs[i-1]
        result = 0
        binom = 1
        for i, d in enumerate(diffs):
            result += binom * d
            if i < len(diffs)-1:
                binom = binom * (c_val - 2 - i) // (i + 1)
        return result

    ext_a = eval_poly([0, 4, 62, 372, 1396, 3980], c) % M
    ext_b = eval_poly([0, 6, 88, 486, 1728, 4750], c) % M

    result = comb(a+b, a) % M
    result = result * (c % M) % M * ((c-1) % M) % M
    result = result * pow(ext_a, a, M) % M
    result = result * pow(ext_b, b, M) % M
    print(result)

solve()

Java

import java.math.BigInteger;

public class Euler194 {
    static long M = 100000000L;
    static BigInteger BM = BigInteger.valueOf(M);

    // Brute-force extension values at c=2,...,7:
    // Type A: 0, 4, 62, 372, 1396, 3980
    // Type B: 0, 6, 88, 486, 1728, 4750
    static long evalPoly(long[] vals, int c) {
        long[] d = vals.clone();
        for (int lev = 1; lev < d.length; lev++)
            for (int i = d.length - 1; i >= lev; i--)
                d[i] -= d[i - 1];
        BigInteger result = BigInteger.ZERO;
        BigInteger binom = BigInteger.ONE;
        for (int i = 0; i < d.length; i++) {
            result = result.add(binom.multiply(BigInteger.valueOf(d[i])));
            if (i < d.length - 1)
                binom = binom.multiply(BigInteger.valueOf(c - 2 - i)).divide(BigInteger.valueOf(i + 1));
        }
        return result.mod(BM).longValue();
    }

    static long mulmod(long a, long b) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BM).longValue();
    }

    static long powmod(long base, int exp) {
        long r = 1;
        base = ((base % M) + M) % M;
        for (; exp > 0; exp >>= 1) {
            if ((exp & 1) == 1)
                r = mulmod(r, base);
            base = mulmod(base, base);
        }
        return r;
    }

    public static void main(String[] args) {
        int a = 25, b = 75, c = 1984;
        long eA = evalPoly(new long[] { 0, 4, 62, 372, 1396, 3980 }, c);
        long eB = evalPoly(new long[] { 0, 6, 88, 486, 1728, 4750 }, c);
        BigInteger comb = BigInteger.ONE;
        for (int i = 1; i <= a; i++)
            comb = comb.multiply(BigInteger.valueOf(a + b - i + 1)).divide(BigInteger.valueOf(i));
        long result = comb.mod(BM).longValue();
        result = mulmod(result, c % M);
        result = mulmod(result, (c - 1) % M);
        result = mulmod(result, powmod(eA, a));
        result = mulmod(result, powmod(eB, b));
        System.out.println(result);
    }
}