Solution

Answer: 640432376

Let

$$F(r,c,n)$$

be the number of proper colourings of the $r\times c$ grid graph using colours from ${1,\dots,n}$, where adjacent means sharing an edge.

We are asked for

$$S(r,c,n)=\sum_{k=1}^{n}F(r,c,k)$$

and specifically

$$S(9,10,1112131415)\pmod{10^9+7}.$$

1. Mathematical analysis

The key observation is:

$F(r,c,n)$ is exactly the chromatic polynomial of the $r\times c$ grid graph evaluated at $n$.

For a graph with $N=rc$ vertices, the chromatic polynomial has degree $N$.

Thus:

$F(9,10,n)$ is a polynomial in $n$ of degree $90$.
Therefore

$$S(9,10,n)=\sum_{k=1}^{n}F(9,10,k)$$

is a polynomial in $n$ of degree $91$.

So the strategy is:

Compute the chromatic polynomial $P(q)=F(9,10,q)$.
Evaluate it at enough small integers.
Use Lagrange interpolation to evaluate $S(9,10,1112131415)$.

The difficult part is step 1.

Frontier / transfer-DP

We process the grid cell-by-cell.

At any moment, only the current “frontier” matters.

For a grid with $r=9$, the frontier has size 9.

Instead of storing actual colours, we store only equality structure.

Example:

[5,7,5,9]

is canonically represented as

[0,1,0,2]

because only equality relations matter.

This drastically compresses the state space.

Polynomial DP

Suppose the current frontier contains $m$ distinct colours.

When colouring the next cell:

reusing an existing frontier colour contributes coefficient $1$,
introducing a brand-new colour contributes a factor

$$(q-m),$$

because there are $q-m$ colours not currently on the frontier.

Thus every DP value becomes a polynomial in $q$.

After all cells are processed, we obtain the full chromatic polynomial.

Why interpolation works

Since $S(n)$ has degree $91$, knowing its values at

$$n=0,1,2,\dots,91$$

determines it uniquely.

We compute those values modulo

$$M=10^9+7$$

and then evaluate the polynomial at

$$1112131415$$

using standard Lagrange interpolation over a finite field.

2. Python implementation

MOD = 1_000_000_007
EMPTY = -1

def poly_eval(coeffs, x, mod=MOD):
    """Evaluate polynomial sum coeffs[i] * x^i."""
    x %= mod
    res = 0
    p = 1
    for a in coeffs:
        res = (res + a * p) % mod
        p = (p * x) % mod
    return res

def lagrange_eval(values, x, mod=MOD):
    """
    Given values[i] = f(i) for i=0..n
    where deg(f) <= n,
    evaluate f(x) modulo mod.
    """
    n = len(values) - 1

    if 0 <= x <= n:
        return values[x] % mod

    x %= mod

    fac = [1] * (n + 1)
    for i in range(1, n + 1):
        fac[i] = fac[i - 1] * i % mod

    invfac = [1] * (n + 1)
    invfac[n] = pow(fac[n], mod - 2, mod)

    for i in range(n, 0, -1):
        invfac[i - 1] = invfac[i] * i % mod

    pre = [1] * (n + 2)
    suf = [1] * (n + 2)

    for i in range(n + 1):
        pre[i + 1] = pre[i] * (x - i) % mod

    for i in range(n, -1, -1):
        suf[i] = suf[i + 1] * (x - i) % mod

    ans = 0

    for i in range(n + 1):
        num = pre[i] * suf[i + 1] % mod

        den = invfac[i] * invfac[n - i] % mod
        if (n - i) & 1:
            den = (-den) % mod

        ans = (ans + values[i] * num % mod * den) % mod

    return ans % mod

def chromatic_poly_grid(r, c, mod=MOD):
    """
    Compute chromatic polynomial of an r x c grid graph.
    Returns coefficients in power basis.
    """

    canon_cache = {}

    def canon(state):
        """Canonical relabelling of colours."""
        if state in canon_cache:
            return canon_cache[state]

        mp = {}
        nxt = 0
        out = []

        for x in state:
            if x == EMPTY:
                out.append(EMPTY)
            else:
                if x not in mp:
                    mp[x] = nxt
                    nxt += 1
                out.append(mp[x])

        res = tuple(out)
        canon_cache[state] = res
        return res

    trans_cache = {}

    def trans(st, i):
        """
        Set frontier position i to every possible label.
        """
        key = (st, i)

        if key in trans_cache:
            return trans_cache[key]

        mx = max(st)
        m = 0 if mx == EMPTY else mx + 1

        base = list(st)

        nexts = [None] * (m + 1)

        for lab in range(m + 1):
            base[i] = lab
            nexts[lab] = canon(tuple(base))

        base[i] = st[i]

        res = (m, nexts)
        trans_cache[key] = res
        return res

    forget_cache = {}

    def forget(st, i):
        key = (st, i)

        if key in forget_cache:
            return forget_cache[key]

        base = list(st)
        base[i] = EMPTY

        res = canon(tuple(base))
        forget_cache[key] = res

        return res

    def add_poly(tgt, poly):
        if len(tgt) < len(poly):
            tgt.extend([0] * (len(poly) - len(tgt)))

        for j, v in enumerate(poly):
            tgt[j] = (tgt[j] + v) % mod

    def add_q_minus_m_mul(tgt, poly, m):
        """
        tgt += (q - m) * poly
        """
        need = len(poly) + 1

        if len(tgt) < need:
            tgt.extend([0] * (need - len(tgt)))

        mm = m % mod

        for j, v in enumerate(poly):
            tgt[j] = (tgt[j] - mm * v) % mod
            tgt[j + 1] = (tgt[j + 1] + v) % mod

    start = tuple([EMPTY] * r)

    dp = {start: [1]}

    for col in range(c):
        for row in range(r):

            i = row

            has_up = row > 0
            has_left = col > 0

            ndp = {}

            for st, poly in dp.items():

                m, nexts = trans(st, i)

                up = st[i - 1] if has_up else -2
                left = st[i] if has_left else -2

                # reuse existing frontier colour
                for lab in range(m):

                    if lab == up or lab == left:
                        continue

                    ns = nexts[lab]

                    tgt = ndp.setdefault(ns, [])
                    add_poly(tgt, poly)

                # introduce new colour
                ns = nexts[m]

                tgt = ndp.setdefault(ns, [])
                add_q_minus_m_mul(tgt, poly, m)

            dp = ndp

    # forget remaining frontier
    for i in range(r):

        ndp = {}

        for st, poly in dp.items():

            ns = forget(st, i)

            tgt = ndp.setdefault(ns, [])
            add_poly(tgt, poly)

        dp = ndp

    return dp[start]

def solve():

    # sanity checks from problem statement

    p22 = chromatic_poly_grid(2, 2)

    assert poly_eval(p22, 3) == 18
    assert poly_eval(p22, 20) == 130340

    p34 = chromatic_poly_grid(3, 4)

    assert poly_eval(p34, 6) == 102923670

    p44 = chromatic_poly_grid(4, 4)

    s44 = sum(poly_eval(p44, k) for k in range(1, 16)) % MOD

    assert s44 == 325951319

    # target
    r, c, n = 9, 10, 1112131415

    p = chromatic_poly_grid(r, c)

    # S(n) has degree 91
    prefix = [0] * 92

    acc = 0

    for k in range(1, 92):
        acc = (acc + poly_eval(p, k)) % MOD
        prefix[k] = acc

    ans = lagrange_eval(prefix, n)

    print(ans)

solve()

3. Code walkthrough

`canon(state)`

Replaces arbitrary colour labels by canonical labels:

[7,3,7,9] -> [0,1,0,2]

This removes irrelevant symmetry.

`dp`

dp[state] stores a polynomial in $q$.

The coefficient list represents:

$$a_0 + a_1 q + a_2 q^2 + \cdots$$

Reusing colours

If we reuse an existing frontier colour, the number of choices is fixed, so we simply copy the polynomial.

New colour

If there are currently $m$ active colours on the frontier, then introducing a new colour contributes:

$$(q-m)$$

Hence:

add_q_minus_m_mul(...)

implements multiplication by $(q-m)$.

Interpolation

Since $S(n)$ has degree $91$, values at $0..91$ uniquely determine it.

lagrange_eval() evaluates that polynomial at

$$n = 1112131415$$

modulo $10^9+7$.

4. Verification

The program reproduces all given checks:

$F(2,2,3)=18$
$F(2,2,20)=130340$
$F(3,4,6)=102923670$
$S(4,4,15)\equiv325951319\pmod{10^9+7}$

Therefore the implementation is consistent with the statement.

Running the full computation gives:

$$S(9,10,1112131415)\equiv 640432376 \pmod{10^9+7}$$

Answer: 640432376

Project Euler Problem 544