Detailed mathematical approach

Problem Summary

Let \(G_{R,C}\) be the \(R\times C\) grid graph, and let \(F(R,C,q)\) denote its chromatic polynomial, meaning the number of proper colorings of the grid with \(q\) available colors. The target quantity is

$S(R,C,N)=\sum_{x=1}^{N}F(R,C,x)\pmod{10^9+7}.$

In the actual problem one has \(R=9\), \(C=10\), and a very large upper limit \(N=1112131415\). Brute-force coloring is hopeless, and even evaluating the chromatic polynomial separately for every \(x\le N\) is impossible. The successful strategy is therefore split into two layers: first compute the full polynomial \(F(R,C,q)\), then sum that polynomial over \(q=1,2,\dots,N\) by turning the result into a short list of power sums.

Mathematical Approach

Write \(D=RC\). Since the grid has \(D\) vertices, \(F(R,C,q)\) is a polynomial in \(q\) of degree at most \(D\). The implementation constructs this polynomial by sweeping a narrow frontier through the grid and storing only the equality pattern among colors that are still relevant to future cells.

Step 1: Sweep the Grid with a Narrow Frontier

The grid is processed one cell at a time. At any moment, only a small set of boundary vertices can still interact with unprocessed cells; these vertices form the frontier. If \(w=\min(R,C)\), then the frontier never contains more than \(w\) vertices, so the exponential part of the search depends on the smaller dimension instead of the full \(RC\).

A frontier state does not remember actual color numbers. It remembers only which frontier positions are forced to carry the same color. In other words, a state is a set partition of the active frontier positions. Two labelings that differ only by renaming the blocks describe the same coloring constraints, so they are canonically relabeled and merged into a single state.

Step 2: Add Each Edge by Deletion-Contraction

Whenever the sweep reaches a new adjacency, the chromatic polynomial satisfies the standard identity

$P_{G+e}(q)=P_G(q)-P_{G/e}(q).$

The first term counts all colorings before the new restriction is enforced. The second term removes the colorings in which the two endpoints receive the same color. Inside the frontier representation, “the endpoints receive the same color” means that their two frontier blocks are identified. Therefore every edge update contributes one positive copy of the current state and one negative copy of the state obtained by merging the two relevant blocks.

If the endpoints are already in the same block, the two terms cancel completely. That is exactly right: a partial coloring in which adjacent vertices are already equal cannot survive after the edge is inserted.

Step 3: Remove Frontier Vertices at the Correct Time

After a frontier vertex has received all edges that can still affect future cells, it may leave the frontier. At that point there are two cases.

If its color class appears nowhere else on the frontier, then that abstract class will never interact with the future again. It can therefore be assigned any of the \(q\) concrete colors, contributing a multiplicative factor \(q\).

If the same class still appears elsewhere on the frontier, then the choice of that concrete color is not free yet; the remaining representative still carries the constraint forward. In that case the state is simply shrunk and canonically relabeled, with no extra factor.

Step 4: Recover the Whole Chromatic Polynomial

The process starts from the empty frontier with coefficient \(1\). As the sweep advances, it repeatedly performs three local actions: introduce a new frontier vertex, add the horizontal and vertical edges that now become known, and retire frontier vertices that can no longer interact with the unseen part of the grid.

Each state carries a polynomial in \(q\), and coefficients of identical canonical states are added modulo \(10^9+7\). After the entire grid has been swept and the last frontier vertices have been removed, only the empty state remains. Its coefficient array is exactly

$F(R,C,q)=\sum_{d=0}^{D} a_d q^d.$

Step 5: Convert the Final Sum into Power Sums

Once the coefficients \(a_d\) are known, the required quantity becomes

$S(R,C,N)=\sum_{x=1}^{N}F(R,C,x)=\sum_{d=0}^{D} a_d\sum_{x=1}^{N}x^d \pmod{10^9+7}.$

So the remaining task is to evaluate the power sums \(\sum_{x=1}^{N}x^d\) for \(0\le d\le D\). For fixed \(d\), Faulhaber's theorem tells us that this is a polynomial in \(N\) of degree \(d+1\). Therefore it is enough to know its values at \(N=0,1,\dots,d+1\), compute those small prefix sums directly, and then recover the value at the huge target \(N\) by Lagrange interpolation modulo the prime \(10^9+7\).

Step 6: Worked Example on the \(2\times 2\) Grid

The \(2\times 2\) grid is a 4-cycle, so its chromatic polynomial is

$F(2,2,q)=q(q-1)^2(q-2)=q^4-4q^3+6q^2-3q.$

This matches the small checkpoints used by the implementation:

$F(2,2,3)=18,\qquad F(2,2,20)=130340.$

The summation stage then becomes

$S(2,2,N)=\sum_{x=1}^{N}\left(x^4-4x^3+6x^2-3x\right),$

which is exactly a linear combination of the power sums for degrees \(1,2,3,4\). The full \(9\times 10\) problem uses the same idea, only with degree up to \(90\) instead of \(4\).

How the Code Works

The C++, Python, and Java implementations follow the same computational pipeline. They first orient the sweep so that the frontier width is as small as possible. They then maintain a dictionary from canonical frontier partitions to coefficient arrays of length \(RC+1\). Introducing a new cell appends a new singleton block to the frontier, inserting a grid edge applies the positive-minus-merged transition dictated by deletion-contraction, and retiring a frontier vertex either multiplies the polynomial by \(q\) or simply removes that position, depending on whether its color class is still present elsewhere on the frontier.

After the empty-state coefficients \((a_0,a_1,\dots,a_D)\) have been obtained, the implementation computes the final answer degree by degree. For each \(d\), it builds the short table of prefix values of \(\sum_{x=1}^{n}x^d\), precomputes factorials and inverse factorials modulo the prime modulus, and evaluates the corresponding degree-\((d+1)\) polynomial at \(N\) using Lagrange interpolation. The last step is the modular accumulation of \(a_d\) times that power sum.

Complexity Analysis

Let \(w=\min(R,C)\), let \(D=RC\), and let \(T_w\) denote the number of reachable canonical frontier states during the sweep. Every transition updates a coefficient array of length \(D+1\), so the dynamic program costs roughly \(O(RC\cdot T_w\cdot D)\) modular arithmetic, with the real constant determined by how many states are created and merged at each sweep position. Memory usage is \(O(T_w\cdot D)\).

The post-processing phase is much smaller. Evaluating all power sums by interpolation takes \(O(D^2)\) time and \(O(D)\) extra memory. Here \(D=90\), so this final stage is cheap; the true bottleneck is the frontier state space, not the interpolation.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=544
2. Chromatic polynomial: Wikipedia — Chromatic polynomial
3. Deletion-contraction formula: Wikipedia — Deletion-contraction formula
4. Lagrange polynomial: Wikipedia — Lagrange polynomial
5. Faulhaber's formula: Wikipedia — Faulhaber's formula