Problem 544: Chromatic Conundrum

Project Euler Problem 544 Solution

EulerSolve provides an optimized solution for Project Euler Problem 544, Chromatic Conundrum, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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

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

Problem 544 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <thread>
#include <unordered_map>
#include <algorithm>
#include <string>
#include <vector>

using namespace std;

constexpr int64_t MOD = 1000000007LL;
constexpr int MAX_DEG = 90;
constexpr int MAX_FRONTIER = 10;

using Poly = array<int64_t, MAX_DEG + 1>;

struct State {
    uint8_t size = 0;
    uint8_t num_labels = 0;
    array<uint8_t, MAX_FRONTIER + 1> labels{};
};

static uint64_t encode_state(const State& s) {
    uint64_t key = s.size;
    for (uint8_t i = 0; i < s.size; ++i) {
        key = (key << 4) | s.labels[i];
    }
    return key;
}

static State canonicalize(const State& s) {
    State out = s;
    array<uint8_t, 16> remap;
    remap.fill(0xFF);
    uint8_t next = 0;
    for (uint8_t i = 0; i < out.size; ++i) {
        uint8_t lbl = out.labels[i];
        if (remap[lbl] == 0xFF) remap[lbl] = next++;
        out.labels[i] = remap[lbl];
    }
    out.num_labels = next;
    return out;
}

static inline void poly_add(Poly& dest, const Poly& src) {
    for (int i = 0; i <= MAX_DEG; ++i) {
        int64_t v = dest[i] + src[i];
        if (v >= MOD) v -= MOD;
        dest[i] = v;
    }
}

static inline void poly_sub(Poly& dest, const Poly& src) {
    for (int i = 0; i <= MAX_DEG; ++i) {
        int64_t v = dest[i] - src[i];
        if (v < 0) v += MOD;
        dest[i] = v;
    }
}

static inline void poly_mul_q(Poly& p) {
    for (int i = MAX_DEG - 1; i >= 0; --i) {
        p[i + 1] = p[i];
    }
    p[0] = 0;
}

struct StateMap {
    vector<State> states;
    vector<Poly> polys;
    unordered_map<uint64_t, size_t> index;

    void reserve(size_t n) {
        states.reserve(n);
        polys.reserve(n);
        index.reserve(n * 2);
    }

    void add_state(const State& s, const Poly& p, bool negate) {
        uint64_t key = encode_state(s);
        auto it = index.find(key);
        if (it == index.end()) {
            states.push_back(s);
            Poly poly = p;
            if (negate) {
                for (int i = 0; i <= MAX_DEG; ++i) {
                    if (poly[i] != 0) poly[i] = MOD - poly[i];
                }
            }
            polys.push_back(poly);
            index.emplace(key, states.size() - 1);
        } else {
            Poly& dest = polys[it->second];
            if (negate) {
                poly_sub(dest, p);
            } else {
                poly_add(dest, p);
            }
        }
    }

    void merge_from(const StateMap& other) {
        for (size_t i = 0; i < other.states.size(); ++i) {
            add_state(other.states[i], other.polys[i], false);
        }
    }
};

constexpr size_t kParallelThreshold = 2000;

template <typename Func>
static StateMap transform_parallel(const StateMap& cur, int threads, size_t reserve_hint, Func func) {
    if (threads <= 1 || cur.states.size() < kParallelThreshold) {
        StateMap out;
        out.reserve(reserve_hint);
        for (size_t i = 0; i < cur.states.size(); ++i) {
            func(cur.states[i], cur.polys[i], out);
        }
        return out;
    }

    int t = min<int>(threads, static_cast<int>(cur.states.size()));
    vector<StateMap> locals(static_cast<size_t>(t));
    vector<thread> workers;
    workers.reserve(static_cast<size_t>(t));

    size_t total = cur.states.size();
    size_t chunk = (total + static_cast<size_t>(t) - 1) / static_cast<size_t>(t);

    for (int ti = 0; ti < t; ++ti) {
        size_t start = static_cast<size_t>(ti) * chunk;
        size_t end = min(total, start + chunk);
        if (start >= end) continue;
        locals[static_cast<size_t>(ti)].reserve((end - start) * 2);
        workers.emplace_back([&, start, end, ti]() {
            StateMap& out = locals[static_cast<size_t>(ti)];
            for (size_t i = start; i < end; ++i) {
                func(cur.states[i], cur.polys[i], out);
            }
        });
    }

    for (auto& w : workers) w.join();

    StateMap merged;
    merged.reserve(reserve_hint);
    for (const auto& local : locals) {
        merged.merge_from(local);
    }
    return merged;
}

static StateMap add_vertex_end(const StateMap& cur) {
    StateMap out;
    out.reserve(cur.states.size());
    for (size_t i = 0; i < cur.states.size(); ++i) {
        State ns = cur.states[i];
        ns.labels[ns.size] = ns.num_labels;
        ns.size++;
        ns.num_labels++;
        out.add_state(ns, cur.polys[i], false);
    }
    return out;
}

static StateMap add_edge(const StateMap& cur, int pos_a, int pos_b, int threads) {
    size_t reserve_hint = cur.states.size() * 2;
    auto func = [pos_a, pos_b](const State& st, const Poly& poly, StateMap& out) {
        out.add_state(st, poly, false);
        State merged = st;
        uint8_t la = merged.labels[static_cast<size_t>(pos_a)];
        uint8_t lb = merged.labels[static_cast<size_t>(pos_b)];
        if (la != lb) {
            for (uint8_t i = 0; i < merged.size; ++i) {
                if (merged.labels[i] == lb) merged.labels[i] = la;
            }
            merged = canonicalize(merged);
        }
        out.add_state(merged, poly, true);
    };
    return transform_parallel(cur, threads, reserve_hint, func);
}

static StateMap remove_pos_swap(const StateMap& cur, int pos, int threads) {
    size_t reserve_hint = cur.states.size();
    auto func = [pos](const State& st, const Poly& poly, StateMap& out) {
        uint8_t label = st.labels[static_cast<size_t>(pos)];
        bool unique = true;
        for (uint8_t i = 0; i < st.size; ++i) {
            if (i != pos && st.labels[i] == label) {
                unique = false;
                break;
            }
        }
        State ns = st;
        uint8_t last = static_cast<uint8_t>(ns.size - 1);
        ns.labels[static_cast<size_t>(pos)] = ns.labels[last];
        ns.size--;
        ns = canonicalize(ns);
        Poly p = poly;
        if (unique) poly_mul_q(p);
        out.add_state(ns, p, false);
    };
    return transform_parallel(cur, threads, reserve_hint, func);
}

static StateMap remove_last(const StateMap& cur, int threads) {
    size_t reserve_hint = cur.states.size();
    auto func = [](const State& st, const Poly& poly, StateMap& out) {
        uint8_t pos = static_cast<uint8_t>(st.size - 1);
        uint8_t label = st.labels[pos];
        bool unique = true;
        for (uint8_t i = 0; i + 1 < st.size; ++i) {
            if (st.labels[i] == label) {
                unique = false;
                break;
            }
        }
        State ns = st;
        ns.size--;
        ns = canonicalize(ns);
        Poly p = poly;
        if (unique) poly_mul_q(p);
        out.add_state(ns, p, false);
    };
    return transform_parallel(cur, threads, reserve_hint, func);
}

static Poly chromatic_poly_grid(int rows, int cols, int threads) {
    if (rows > cols) swap(rows, cols);
    StateMap cur;
    cur.reserve(1);

    State empty;
    Poly base{};
    base[0] = 1;
    cur.add_state(empty, base, false);

    int frontier = 0;
    for (int c = 0; c < cols; ++c) {
        for (int r = 0; r < rows; ++r) {
            cur = add_vertex_end(cur);
            frontier += 1;
            int new_pos = frontier - 1;
            if (c > 0) {
                cur = add_edge(cur, r, new_pos, threads);
            }
            if (r > 0) {
                cur = add_edge(cur, r - 1, new_pos, threads);
            }
            if (c > 0) {
                cur = remove_pos_swap(cur, r, threads);
                frontier -= 1;
            }
        }
    }

    while (frontier > 0) {
        cur = remove_last(cur, threads);
        frontier -= 1;
    }

    Poly result{};
    for (size_t i = 0; i < cur.states.size(); ++i) {
        if (cur.states[i].size == 0) {
            result = cur.polys[i];
            break;
        }
    }
    return result;
}

static int64_t pow_mod(int64_t a, int64_t e) {
    int64_t res = 1 % MOD;
    a %= MOD;
    while (e > 0) {
        if (e & 1) res = (res * a) % MOD;
        a = (a * a) % MOD;
        e >>= 1;
    }
    return res;
}

static int64_t eval_poly(const Poly& coeffs, int max_deg, int64_t n) {
    int64_t res = 0;
    int64_t x = n % MOD;
    for (int d = max_deg; d >= 0; --d) {
        res = (res * x + coeffs[d]) % MOD;
    }
    return res;
}

static int64_t lagrange_eval(const vector<int64_t>& y,
                             int64_t x,
                             const vector<int64_t>& fact,
                             const vector<int64_t>& inv_fact) {
    int m = static_cast<int>(y.size()) - 1;
    if (x <= m) return y[static_cast<size_t>(x)];

    vector<int64_t> prefix(static_cast<size_t>(m) + 1);
    vector<int64_t> suffix(static_cast<size_t>(m) + 1);
    prefix[0] = 1;
    for (int i = 0; i < m; ++i) {
        int64_t term = x - i;
        if (term < 0) term += MOD;
        prefix[static_cast<size_t>(i + 1)] = (prefix[static_cast<size_t>(i)] * term) % MOD;
    }
    suffix[static_cast<size_t>(m)] = 1;
    for (int i = m - 1; i >= 0; --i) {
        int64_t term = x - (i + 1);
        if (term < 0) term += MOD;
        suffix[static_cast<size_t>(i)] = (suffix[static_cast<size_t>(i + 1)] * term) % MOD;
    }

    int64_t res = 0;
    for (int i = 0; i <= m; ++i) {
        int64_t num = prefix[static_cast<size_t>(i)] * suffix[static_cast<size_t>(i)] % MOD;
        int64_t denom = inv_fact[static_cast<size_t>(i)] * inv_fact[static_cast<size_t>(m - i)] % MOD;
        int64_t term = y[static_cast<size_t>(i)] * num % MOD * denom % MOD;
        if ((m - i) & 1) term = MOD - term;
        res += term;
        if (res >= MOD) res -= MOD;
    }
    return res;
}

static int64_t sum_pows(int d, int64_t n, const vector<int64_t>& fact, const vector<int64_t>& inv_fact) {
    if (n <= 0) return 0;
    if (d == 0) return n % MOD;
    int m = d + 1;
    vector<int64_t> y(static_cast<size_t>(m) + 1);
    y[0] = 0;
    for (int i = 1; i <= m; ++i) {
        y[static_cast<size_t>(i)] = (y[static_cast<size_t>(i - 1)] + pow_mod(i, d)) % MOD;
    }
    int64_t x = n % MOD;
    return lagrange_eval(y, x, fact, inv_fact);
}

static int64_t compute_S(const Poly& coeffs, int max_deg, int64_t n,
                         const vector<int64_t>& fact, const vector<int64_t>& inv_fact) {
    int64_t res = 0;
    for (int d = 0; d <= max_deg; ++d) {
        if (coeffs[d] == 0) continue;
        int64_t sum = sum_pows(d, n, fact, inv_fact);
        res = (res + coeffs[d] * sum) % MOD;
    }
    return res;
}

static bool run_validation(int threads, const vector<int64_t>& fact, const vector<int64_t>& inv_fact) {
    bool ok = true;

    {
        Poly p = chromatic_poly_grid(2, 2, threads);
        int64_t f3 = eval_poly(p, 4, 3);
        int64_t f20 = eval_poly(p, 4, 20);
        if (f3 != 18) {
            cerr << "Validation failed: F(2,2,3) got " << f3 << ", expected 18\n";
            ok = false;
        }
        if (f20 != 130340) {
            cerr << "Validation failed: F(2,2,20) got " << f20 << ", expected 130340\n";
            ok = false;
        }
    }

    {
        Poly p = chromatic_poly_grid(3, 4, threads);
        int64_t f6 = eval_poly(p, 12, 6);
        if (f6 != 102923670) {
            cerr << "Validation failed: F(3,4,6) got " << f6 << ", expected 102923670\n";
            ok = false;
        }
    }

    {
        Poly p = chromatic_poly_grid(4, 4, threads);
        int64_t s15 = compute_S(p, 16, 15, fact, inv_fact);
        if (s15 != 325951319) {
            cerr << "Validation failed: S(4,4,15) got " << s15 << ", expected 325951319\n";
            ok = false;
        }
    }

    if (ok) {
        cerr << "Validation checkpoints passed.\n";
    }
    return ok;
}

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    bool validate = true;
    unsigned hw = thread::hardware_concurrency();
    int threads = hw ? static_cast<int>(hw) : 1;
    threads = max(1, min(threads, 8));

    for (int i = 1; i < argc; ++i) {
        string arg(argv[i]);
        if (arg == "--no-validate") {
            validate = false;
        } else if (arg == "--threads" && i + 1 < argc) {
            threads = max(1, stoi(argv[++i]));
        }
    }

    vector<int64_t> fact(static_cast<size_t>(MAX_DEG) + 3, 1);
    vector<int64_t> inv_fact(static_cast<size_t>(MAX_DEG) + 3, 1);
    for (size_t i = 1; i < fact.size(); ++i) fact[i] = fact[i - 1] * static_cast<int64_t>(i) % MOD;
    inv_fact.back() = pow_mod(fact.back(), MOD - 2);
    for (size_t i = fact.size() - 1; i > 0; --i) {
        inv_fact[i - 1] = inv_fact[i] * static_cast<int64_t>(i) % MOD;
    }

    if (validate && !run_validation(min(threads, 2), fact, inv_fact)) {
        return 1;
    }

    constexpr int R = 9;
    constexpr int C = 10;
    constexpr int64_t N = 1112131415LL;
    Poly poly = chromatic_poly_grid(R, C, threads);
    int64_t answer = compute_S(poly, R * C, N, fact, inv_fact);
    cout << answer << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


if __name__ == "__main__":
    print(solve())

Java

import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;

public class Euler544 {
    static final long MOD = 1000000007L;
    static final int MAX_DEG = 90;

    static long encodeState(int size, byte[] labels) {
        long key = size;
        for (int i = 0; i < size; i++) {
            key = (key << 4) | labels[i];
        }
        return key;
    }

    static byte[] decodeLabels(long key, int size) {
        byte[] labels = new byte[size];
        for (int i = size - 1; i >= 0; i--) {
            labels[i] = (byte) (key & 0xF);
            key >>= 4;
        }
        return labels;
    }

    static int decodeSize(long key, int originalSize) {
        return originalSize;
    }

    static class CanonicalResult {
        int size;
        int numLabels;
        long key;

        CanonicalResult(int size, int numLabels, long key) {
            this.size = size;
            this.numLabels = numLabels;
            this.key = key;
        }
    }

    static CanonicalResult canonicalize(int size, byte[] labels) {
        byte[] remap = new byte[16];
        Arrays.fill(remap, (byte) -1);
        byte next = 0;
        byte[] out = new byte[size];
        for (int i = 0; i < size; i++) {
            byte lbl = labels[i];
            if (remap[lbl] == -1)
                remap[lbl] = next++;
            out[i] = remap[lbl];
        }
        return new CanonicalResult(size, next, encodeState(size, out));
    }

    static void polyAdd(long[] dest, long[] src) {
        for (int i = 0; i <= MAX_DEG; i++) {
            dest[i] = (dest[i] + src[i]) % MOD;
        }
    }

    static void polySub(long[] dest, long[] src) {
        for (int i = 0; i <= MAX_DEG; i++) {
            dest[i] = (dest[i] - src[i] + MOD) % MOD;
        }
    }

    static long[] polyMulQ(long[] p) {
        long[] out = new long[MAX_DEG + 1];
        System.arraycopy(p, 0, out, 1, MAX_DEG);
        return out;
    }

    static class StateInfo {
        int size;
        int numLabels;
        long[] poly;

        StateInfo(int size, int numLabels, long[] poly) {
            this.size = size;
            this.numLabels = numLabels;
            this.poly = poly;
        }
    }

    static class StateMap {
        Map<Long, StateInfo> map = new HashMap<>();

        void addState(int size, int numLabels, long key, long[] p, boolean negate) {
            long[] pToAdd;
            if (negate) {
                pToAdd = new long[MAX_DEG + 1];
                for (int i = 0; i <= MAX_DEG; i++) {
                    pToAdd[i] = p[i] == 0 ? 0 : MOD - p[i];
                }
            } else {
                pToAdd = p;
            }

            StateInfo info = map.get(key);
            if (info == null) {
                map.put(key, new StateInfo(size, numLabels, pToAdd.clone()));
            } else {
                polyAdd(info.poly, pToAdd);
            }
        }
    }

    static StateMap addVertexEnd(StateMap cur) {
        StateMap out = new StateMap();
        for (Map.Entry<Long, StateInfo> entry : cur.map.entrySet()) {
            long key = entry.getKey();
            StateInfo info = entry.getValue();
            byte[] labels = decodeLabels(key, info.size);
            byte[] newLabels = Arrays.copyOf(labels, info.size + 1);
            newLabels[info.size] = (byte) info.numLabels;
            long newKey = encodeState(info.size + 1, newLabels);
            out.addState(info.size + 1, info.numLabels + 1, newKey, info.poly, false);
        }
        return out;
    }

    static StateMap addEdge(StateMap cur, int posA, int posB) {
        StateMap out = new StateMap();
        for (Map.Entry<Long, StateInfo> entry : cur.map.entrySet()) {
            long key = entry.getKey();
            StateInfo info = entry.getValue();
            out.addState(info.size, info.numLabels, key, info.poly, false);

            byte[] labels = decodeLabels(key, info.size);
            byte la = labels[posA];
            byte lb = labels[posB];
            if (la != lb) {
                for (int i = 0; i < info.size; i++) {
                    if (labels[i] == lb)
                        labels[i] = la;
                }
                CanonicalResult cr = canonicalize(info.size, labels);
                out.addState(cr.size, cr.numLabels, cr.key, info.poly, true);
            } else {
                out.addState(info.size, info.numLabels, key, info.poly, true);
            }
        }
        return out;
    }

    static StateMap removePosSwap(StateMap cur, int pos) {
        StateMap out = new StateMap();
        for (Map.Entry<Long, StateInfo> entry : cur.map.entrySet()) {
            long key = entry.getKey();
            StateInfo info = entry.getValue();
            byte[] labels = decodeLabels(key, info.size);
            byte label = labels[pos];
            boolean unique = true;
            for (int i = 0; i < info.size; i++) {
                if (i != pos && labels[i] == label) {
                    unique = false;
                    break;
                }
            }

            byte[] nsLabels = new byte[info.size - 1];
            for (int i = 0; i < info.size - 1; i++) {
                if (i == pos) {
                    nsLabels[i] = labels[info.size - 1];
                } else {
                    nsLabels[i] = labels[i];
                }
            }

            CanonicalResult cr = canonicalize(info.size - 1, nsLabels);
            long[] p = unique ? polyMulQ(info.poly) : info.poly;
            out.addState(cr.size, cr.numLabels, cr.key, p, false);
        }
        return out;
    }

    static StateMap removeLast(StateMap cur) {
        StateMap out = new StateMap();
        for (Map.Entry<Long, StateInfo> entry : cur.map.entrySet()) {
            long key = entry.getKey();
            StateInfo info = entry.getValue();
            byte[] labels = decodeLabels(key, info.size);
            int pos = info.size - 1;
            byte label = labels[pos];
            boolean unique = true;
            for (int i = 0; i < pos; i++) {
                if (labels[i] == label) {
                    unique = false;
                    break;
                }
            }

            byte[] nsLabels = Arrays.copyOf(labels, info.size - 1);
            CanonicalResult cr = canonicalize(info.size - 1, nsLabels);
            long[] p = unique ? polyMulQ(info.poly) : info.poly;
            out.addState(cr.size, cr.numLabels, cr.key, p, false);
        }
        return out;
    }

    static long[] chromaticPolyGrid(int rows, int cols) {
        if (rows > cols) {
            int temp = rows;
            rows = cols;
            cols = temp;
        }

        StateMap cur = new StateMap();
        long[] base = new long[MAX_DEG + 1];
        base[0] = 1;
        cur.addState(0, 0, encodeState(0, new byte[0]), base, false);

        int frontier = 0;
        for (int c = 0; c < cols; c++) {
            for (int r = 0; r < rows; r++) {
                cur = addVertexEnd(cur);
                frontier++;
                int newPos = frontier - 1;
                if (c > 0)
                    cur = addEdge(cur, r, newPos);
                if (r > 0)
                    cur = addEdge(cur, r - 1, newPos);
                if (c > 0) {
                    cur = removePosSwap(cur, r);
                    frontier--;
                }
            }
        }

        while (frontier > 0) {
            cur = removeLast(cur);
            frontier--;
        }

        for (StateInfo info : cur.map.values()) {
            if (info.size == 0)
                return info.poly;
        }
        return new long[MAX_DEG + 1];
    }

    static long lagrangeEval(long[] y, long x, long[] fact, long[] invFact) {
        int m = y.length - 1;
        if (x <= m)
            return y[(int) x];

        long[] prefix = new long[m + 2];
        long[] suffix = new long[m + 2];
        prefix[0] = 1;
        for (int i = 0; i <= m; i++) {
            long term = (x - i) % MOD;
            if (term < 0)
                term += MOD;
            prefix[i + 1] = (prefix[i] * term) % MOD;
        }
        suffix[m + 1] = 1;
        for (int i = m; i >= 0; i--) {
            long term = (x - i) % MOD;
            if (term < 0)
                term += MOD;
            suffix[i] = (suffix[i + 1] * term) % MOD;
        }

        long res = 0;
        for (int i = 0; i <= m; i++) {
            long num = (prefix[i] * suffix[i + 1]) % MOD;
            long denom = (invFact[i] * invFact[m - i]) % MOD;
            long term = (((y[i] * num) % MOD) * denom) % MOD;
            if ((m - i) % 2 != 0)
                term = MOD - term;
            res = (res + term) % MOD;
        }
        return res;
    }

    static long powMod(long a, long e) {
        long res = 1;
        a %= MOD;
        while (e > 0) {
            if ((e & 1) != 0)
                res = (res * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1;
        }
        return res;
    }

    static long sumPows(int d, long n, long[] fact, long[] invFact) {
        if (n <= 0)
            return 0;
        if (d == 0)
            return n % MOD;
        int m = d + 1;
        long[] y = new long[m + 1];
        for (int i = 1; i <= m; i++) {
            y[i] = (y[i - 1] + powMod(i, d)) % MOD;
        }
        return lagrangeEval(y, n % MOD, fact, invFact);
    }

    static long computeS(long[] coeffs, int maxDeg, long n, long[] fact, long[] invFact) {
        long res = 0;
        for (int d = 0; d <= maxDeg; d++) {
            if (coeffs[d] == 0)
                continue;
            long sumVal = sumPows(d, n, fact, invFact);
            res = (res + coeffs[d] * sumVal) % MOD;
        }
        return res;
    }

    public static String solve() {
        long[] fact = new long[MAX_DEG + 3];
        long[] invFact = new long[MAX_DEG + 3];
        fact[0] = 1;
        for (int i = 1; i < fact.length; i++)
            fact[i] = (fact[i - 1] * i) % MOD;
        invFact[invFact.length - 1] = powMod(fact[fact.length - 1], MOD - 2);
        for (int i = invFact.length - 1; i > 0; i--) {
            invFact[i - 1] = (invFact[i] * i) % MOD;
        }

        int R = 9;
        int C = 10;
        long N = 1112131415L;

        long[] poly = chromaticPolyGrid(R, C);
        return Long.toString(computeS(poly, R * C, N, fact, invFact));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}