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

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. The rest of the unit only sees the incoming pair through adjacency constraints.

This is the invariant that makes the chain multiplicative. When one unit is glued to the next, its right boundary becomes the next unit's left boundary, and that shared boundary is again just a distinct ordered pair. So every unit can be counted independently once the local extension count is known.

Local Extension Counts for Type A and Type B

Let $A(c)$ be the number of valid colorings of one type-A unit when the incoming boundary colors are fixed and distinct, and let $B(c)$ be the analogous number for one type-B unit.

The two unit graphs differ in exactly one place. In type A, the outgoing lower boundary vertex is also adjacent to the incoming lower boundary vertex, so it is forbidden to reuse that color. Type B omits that edge, which is why $B(c)$ is larger.

Since a unit has five unfixed vertices, these local counts are degree-5 coloring polynomials in $c$. The implementations determine them exactly by brute force for $c=2,3,4,5,6,7$:

$$\bigl(A(2),A(3),A(4),A(5),A(6),A(7)\bigr)=(0,4,62,372,1396,3980),$$

$$\bigl(B(2),B(3),B(4),B(5),B(6),B(7)\bigr)=(0,6,88,486,1728,4750).$$

Six exact values are enough to recover a degree-5 polynomial uniquely, so Newton forward interpolation gives the closed formulas used implicitly by the implementations:

$$A(c)=c^5-9c^4+34c^3-69c^2+77c-38,$$

$$B(c)=c^5-8c^4+27c^3-50c^2+52c-24=(c-2)^3(c^2-2c+3).$$

A Worked Example: The Case $c=3$

Take three colors $\{x,y,z\}$ and fix the incoming boundary as $(x,y)$. A complete enumeration of one unit gives

$$A(3)=4,\qquad B(3)=6.$$

The difference is already visible at this smallest nontrivial palette size. Type B has two extra extensions because its outgoing lower boundary is allowed to reuse the color $y$; type A forbids exactly those colorings through its extra lower edge. This tiny example captures the only structural distinction between the two building blocks.

Assembling the Whole Configuration

Once the local counts are known, the full answer factorizes. First choose which of the $a+b$ positions are occupied by type-A units. That contributes

$$\binom{a+b}{a}.$$

Then choose the ordered color pair on the very first boundary. Because the two colors must be distinct, this contributes

$$c(c-1).$$

Every type-A unit contributes a factor $A(c)$, every type-B unit contributes a factor $B(c)$, and the location of the unit does not matter because every interface is just a distinct ordered pair and the color names are symmetric. Therefore

$$\boxed{N(a,b,c)=\binom{a+b}{a}\,c(c-1)\,A(c)^a B(c)^b.}$$

For the Euler instance $(a,b,c)=(25,75,1984)$, that exact integer is finally reduced modulo $10^8$.

How the Code Works

Exact Enumeration of One Unit

The C++, Python, and Java implementations begin by fixing two distinct incoming boundary colors and brute-forcing the five remaining vertices for $c=2,\dots,7$. That produces the six exact values of $A(c)$ and $B(c)$ listed above. The C++ implementation also checks the interpolated values against fresh brute-force counts for a few additional small palette sizes.

Forward-Difference Evaluation at the Target $c$

Instead of hard-coding symbolic algebra, the implementations evaluate the degree-5 polynomials through Newton forward differences. If $n=c-2$, then

$$A(c)=\sum_{k=0}^{5}\Delta^kA(2)\binom{n}{k},\qquad B(c)=\sum_{k=0}^{5}\Delta^kB(2)\binom{n}{k}.$$

For these two unit counts, the first entries of the forward-difference tables are

$$A:\ (0,4,54,198,264,120),\qquad B:\ (0,6,76,240,288,120).$$

Because the underlying functions really are degree-5 polynomials, this interpolation is exact rather than approximate.

Final Modular Product

After obtaining $A(c)$ and $B(c)$, the implementations compute $\binom{a+b}{a}$, multiply by $c(c-1)$, and raise the two local counts to the powers $a$ and $b$ with fast modular exponentiation. The final stage runs modulo $10^8$, while the small checkpoint calculations may use exact integer arithmetic first and reduce only at the end.

Complexity Analysis

The local brute-force phase is constant-sized: it evaluates only six palette sizes, and each unit has only five unfixed vertices. Building the forward-difference table and evaluating the resulting polynomial are therefore $O(1)$.

For parameterized inputs, forming $\binom{a+b}{a}$ by multiplicative accumulation costs $O(\min(a,b))$ arithmetic steps, and the two modular powers cost $O(\log a+\log b)$. Extra memory usage is $O(1)$ beyond a few fixed tables of sample values and differences.

Footnotes and References

Problem page: https://projecteuler.net/problem=194
Graph coloring: Wikipedia - Graph coloring
Chromatic polynomial: Wikipedia - Chromatic polynomial
Newton polynomial: Wikipedia - Newton polynomial
Finite difference: Wikipedia - Finite difference
Binomial coefficient: Wikipedia - Binomial coefficient
Modular exponentiation: Wikipedia - Modular exponentiation

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);
    }
}