Problem 393: Migrating Ants

Project Euler Problem 393 Solution

EulerSolve provides an optimized solution for Project Euler Problem 393, Migrating Ants, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $F(n)$ denote the number of valid migrations on an $n \times n$ grid. Every cell initially contains one ant, every ant must move to a neighboring cell, every cell must receive exactly one ant, and two ants are not allowed to traverse the same undirected edge in opposite directions. In graph terms, each vertex has out-degree $1$ and in-degree $1$, so a valid migration is a disjoint union of directed cycles, with 2-cycles forbidden. Mathematical Approach Step 1: Convert the migration into two perfect matchings Color the grid as a checkerboard, with black and white cells. Every legal move goes from one color to the other, because adjacent cells always have opposite parity. Split the directed edges of a migration into two sets: $$M_1=\{\text{moves from black to white}\},\qquad M_2=\{\text{moves from white to black}\}.$$ Now inspect $M_1$. Every black vertex has exactly one outgoing move, so every black cell is incident to exactly one edge of $M_1$. Every white vertex has exactly one incoming move, and the only possible source is a black neighbor, so every white cell is also incident to exactly one edge of $M_1$. Therefore $M_1$ is a perfect matching of the bipartite grid graph. By the same argument, $M_2$ is another perfect matching....

Detailed mathematical approach

Problem Summary

An $n\times n$ grid of squares holds $n^2$ ants, one per square. All move at the same moment to a square sharing an edge with the one they leave — usually four choices, fewer at an edge of the grid or in a corner. The move counts only if no two ants end on the same square and no two cross the same edge, and $f(n)$ is the number of moves that count. The statement gives $f(4)=88$ and asks for $f(10)$. ^[1]

The whole solution in one sentence: a migration amounts to splitting all the squares into neighbouring pairs twice over, once for the departures and once for the arrivals, never using the same pair twice — and these double pairings are counted one grid row at a time, remembering only which squares must still join the row below.

Mathematical Approach

Sections 1 to 4 replace the ants by two matchings, turning both prohibitions into one; Sections 5 to 8 build the recurrence.

1. A migration is a permutation without two-cycles

Let $\sigma(v)$ be the square the ant on $v$ moves to; the rules make $\sigma(v)$ a neighbour of $v$. No two ants end on the same square, so $\sigma$ is injective, and an injective map of a finite set into itself is a bijection — a permutation of the $n^2$ squares, which splits into disjoint cycles. ^[2] No cycle has length $1$, since no square neighbours itself, and two ants cross one edge exactly when $\sigma(v)=w$ and $\sigma(w)=v$, a cycle of length $2$. ^[2] So $f(n)$ counts the permutations sending every square to a neighbour with no cycle of length $1$ or $2$: the ants of each cycle march one step along it.

2. A chessboard colouring splits a migration into two matchings

Colour the square in row $r$, column $c$ black when $r+c$ is even and white otherwise. Neighbours differ by one in exactly one coordinate, so every edge joins black to white: the grid graph is bipartite. ^[3]

Let $M_1$ hold the undirected edges $\{b,\sigma(b)\}$ over the black squares $b$, and $M_2$ the edges $\{w,\sigma(w)\}$ over the white squares. Each black square lies on one edge of $M_1$, its own ant's. Each white square $w$ lies on one too: the ant reaching $w$ starts at $\sigma^{-1}(w)$, a neighbour of $w$ and so black, and that move is the only edge of $M_1$ at $w$. So $M_1$ meets every square exactly once, which is what a perfect matching means, ^[4] and likewise $M_2$. An edge $\{b,w\}$ lies in both exactly when $\sigma(b)=w$ and $\sigma(w)=b$, the forbidden two-cycle, so the ban says precisely that $M_1\cap M_2=\varnothing$.

Conversely, given such a pair, move each black square's ant along its $M_1$ edge and each white square's along its $M_2$ edge: every square then receives exactly one ant, and disjointness rules out two-cycles. The constructions undo each other, so $f(n)$ counts the ordered pairs of edge-disjoint perfect matchings of the $n\times n$ grid graph.

3. Cycle covers, and two consequences for free

The union $M_1\cup M_2$ gives every square degree two, so it is a disjoint union of cycles along which the matchings alternate. Cycles of a bipartite graph are even, ^[3] and length $2$ needs a repeated edge, so every cycle has length at least $4$. Conversely a cover of all squares by $c$ disjoint cycles splits into alternating halves in two ways per cycle, independently, giving $2^c$ ordered pairs. Hence

$$f(n)=\sum_{C}2^{c(C)}$$

over the covers $C$ by $c(C)$ disjoint cycles; the factor is the free direction of travel on each cycle.

Two consequences follow. A perfect matching needs an even number of vertices, ^[4] and $n^2$ is odd for odd $n$, so $f(n)=0$ there — the colouring agrees, since $\lceil n^2/2\rceil$ black and $\lfloor n^2/2\rfloor$ white squares cannot be paired off — and the published program returns $0$ for $n=1,3,5,7,9,11$. And $f(n)$ is even for $n\ge2$: exchanging the matchings reverses every arrow and has no fixed point, since disjoint non-empty matchings differ. The computed $f(2)=2$, $f(4)=88$, $f(6)=207408$ and $f(8)=22902801416$ are all even.

4. A worked example for $n=2$, and a check at $n=4$

For $n=2$ the four squares form one square cycle, the only cover, so $c=1$ and $f(2)=2$: the ants rotate clockwise or anticlockwise. Equivalently, the $2\times2$ grid has two perfect matchings, the horizontal pair of edges and the vertical pair; they are disjoint, and either may be $M_1$.

For $n=4$ the formula can be checked by exhaustive search. A cover gives each of the $16$ squares degree two, so it uses $16$ of the grid's $24$ edges, and the $\binom{24}{16}=735471$ sixteen-edge subsets can all be tested: $18$ of them make every square degree two. Following their cycles sorts them by length:

four $4$-cycles: $1$ cover, $2^4=16$;
two $4$-cycles and an $8$-cycle: $4$ covers, $8$ each, so $32$;
a $4$-cycle and a $12$-cycle: $5$ covers, $4$ each, so $20$;
two $8$-cycles: $2$ covers, $4$ each, so $8$;
one $16$-cycle: $6$ covers, $2$ each, so $12$.

The total $16+32+20+8+12=88$ is the statement's value, and the last line cross-checks the enumeration: a cycle through all $16$ squares is a Hamiltonian cycle of the $4\times4$ grid, and there are exactly $6$. ^[7]

5. Why the matchings cannot simply be listed

At $n=4$ brute force settles it: of the $36^2=1296$ ordered pairs of the grid's $36$ perfect matchings, exactly $88$ share no edge. The perfect matchings of the $n\times n$ grid graph are the domino tilings of the board, ^[5] numbering $2$, $36$, $6728$, $12988816$ and $258584046368$ for $n=2,4,6,8,10$, ^[6] in closed form, with $a_j=4\cos^2(\pi j/(n+1))$ and $n$ even,

$$\prod_{j=1}^{n/2}\prod_{k=1}^{n/2}(a_j+a_k),$$

found independently by Temperley and Fisher and by Kasteleyn. ^[5] Squaring it counts all ordered pairs, including those sharing an edge, so it does not impose the condition we need, and at $n=10$ listing the matchings is already out of reach. The published solutions therefore count the compatible pairs by a row recurrence, never holding a whole matching.

6. One row of one matching: the profile and the choices at a square

Work down the rows, both matchings at once. In one matching every edge touching row $r$ comes down from row $r-1$, lies inside row $r$, or goes down into row $r+1$. The rows above are finished and those below untouched, so row $r$ need only know which of its squares are matched from above, and the rows below only which columns send an edge down out of row $r$: two sets of columns, each written as a number whose bit $c$ has the value $2^c$ and is set when column $c$ belongs to the set. Passing such a profile from row to row is the standard broken-profile technique. ^[8] Counting pairs, the state is the pair of incoming masks, packed into one integer $s=i_1+2^{n}i_2$; grouping the fillings of the rows above by this state and adding their counts is dynamic programming, since all their continuations depend on the state alone.

A row is built by sweeping its columns left to right, keeping for each matching a mask of that row's squares already matched and starting from the incoming mask. At column $c$ a matching has three possibilities, named Filled, Horizontal and Vertical in the source:

the square is already matched, from above or by a horizontal edge at column $c-1$, and nothing is decided;
the square and column $c+1$ are both free and are joined horizontally, marking both;
the square is free and a row lies below, and an edge goes down, marking it and recording column $c$ in the outgoing mask.

A free square with neither of the last two kills the branch, and on the last row nothing may go down, forcing an empty outgoing mask.

7. The only two forbidden combinations, and the recurrence

The sweep creates each edge exactly once, at a moment we can name: the horizontal edge between columns $c$ and $c+1$ of row $r$ only from a Horizontal decision at column $c$ of row $r$, the vertical edge below column $c$ of row $r$ only from a Vertical decision there. So the matchings share an edge if and only if both make the same one of those decisions at the same column of the same row. Rejecting Horizontal with Horizontal and Vertical with Vertical at one column is therefore exactly edge-disjointness, and nothing else about the matchings need ever be compared.

Write $\mathbf i=(i_1,i_2)$ and $\mathbf o=(o_1,o_2)$ for a row's incoming and outgoing mask pairs, let $T_r(\mathbf i,\mathbf o)$ count the fillings of row $r$ from $\mathbf i$ to $\mathbf o$, and $D_r(\mathbf i)$ the ways to fill the rows above so that row $r$ receives $\mathbf i$. Row $r$ is independent of everything above it once $\mathbf i$ is fixed, so

$$D_{r+1}(\mathbf o)=\sum_{\mathbf i}D_r(\mathbf i)\,T_r(\mathbf i,\mathbf o).$$

Nothing enters the first row and nothing may leave the last, so $D_0(\mathbf 0)=1$ with $D_0$ zero elsewhere, and $f(n)=D_n(\mathbf 0)$. Only one bit of $r$ reaches $T_r$, whether $r$ is the last row, so two tables serve the whole grid and each entry is built once.

8. Two facts that keep the tables small

The masks $o_1$ and $o_2$ never share a column, since two downward edges in one column would be the same edge. So each column lies in $o_1$, in $o_2$ or in neither: at most $3^n$ states occur instead of the $4^n$ two free masks would allow, $59049$ rather than $1048576$ at $n=10$. Fewer are reached, the horizontal ban coupling the matchings too; counting gives $15$, $93$, $639$ and $4653$ states for $n=4,6,8,10$.

Every stored transition count is $1$. In one matching the squares of $i$ are matched from above and those of $o$ downwards, so the rest must pair horizontally — and a set of columns splits into pairs of neighbours in at most one way, since in each maximal run the first must take the second, the third the fourth, and an odd run cannot split. So $\mathbf i$ and $\mathbf o$ fix the row completely and $T_r(\mathbf i,\mathbf o)$ is $0$ or $1$; instrumenting the published sweep confirms this up to $n=12$, the completed sweeps always numbering the same as the distinct outgoing pairs. Run at $n=2$, the recurrence sends $D_0(\mathbf 0)=1$ to weight $1$ on each of $(0,3)$ and $(3,0)$, then both of those to $\mathbf 0$, giving $f(2)=2$ as in Section 4.

How the Code Works

All three programs implement Sections 6 to 8 with the same two layers and the same two bans, differing only in containers and integer types. ^[9]

In C++ a Solver holds $n$ as n_ and $2^n-1$ as bit_mask_. Its solve() keeps $D_r$ as an unordered_map<int, BigInt> keyed by the packed state $s$, starting from $0\mapsto1$; for each of the $n$ rows it unpacks every key into in1 and in2, asks get_transitions for the reachable outgoing states, accumulates into a fresh map and swaps it in. Afterwards it returns the entry for state $0$, or zero when absent, as happens for odd $n$.

get_transitions(in1, in2, last_row) is the table $T_r$, cached under in1 | (in2 << n_) | (last_row << 24) so that both tables share one map. On a miss it runs dfs(c, occ1, occ2, out1, out2) across the columns, calling options_for_cell once per matching for the Move records of Section 6; the double loop over the two lists skips the forbidden combinations, and at c == n_ increments the counter for the outgoing state.

Arguments and self-checks are C++ only. parse_arguments accepts --n= and --skip-checkpoints and checks $1\le n\le 12$, though its helper parse_int_after_prefix accumulates the digits in an unchecked int, so a long enough digit string overflows before that test: --n=4294967298 wraps round to $2$ and is accepted. Unless skipped, run_checkpoints verifies the hand-derived $f(2)=2$ and the statement's $f(4)=88$, then compares the recurrence at $n=3$ with brute_force_small(3), which tries every neighbour for every square and filters with is_valid_assignment: Manhattan distance $1$ for each move, exactly one ant arriving per square, and to[to[u]] == u never true, the two-cycle ban at its bluntest. A failure exits with status $2$ before printing.

Python is the same algorithm in about sixty lines, with dictionaries, tuples in place of Move and unbounded integers; Java mirrors the C++ with HashMap, an enum MoveType and BigInteger. Only the C++ takes $n$ on the command line: Python's solve accepts it as a parameter but its entry point uses the default, and Java fixes static int n = 10. The arbitrary-precision types earn their place at the top of the accepted range, $f(12)$ having $26$ decimal digits against the $19$ a signed 64-bit integer holds; the multiplicities, always $1$ by Section 8, would fit anywhere.

Complexity and Verification

A sweep has at most two surviving branches per column, since where both matchings are free only two of the four combinations survive, so it visits $O(2^n)$ nodes. With at most $3^n$ incoming states by Section 8 and two values of the last-row flag, filling the cache costs $O(6^n)$ operations and stores $O(6^n)$ outgoing entries; each of the $n$ rows then walks the list of every state it holds, again $O(6^n)$. The total is $O(n\cdot 6^n)$ arithmetic operations — exponential in the side length, linear in the number of rows, and linear in it even without the cache, since rebuilding a list once per row costs the same order. Memory is $O(6^n)$ for the cache, $O(3^n)$ for the two layers of $D_r$ and $O(n)$ for the sweep's stack. Both bounds count operations, not bit work, and assume constant-time table access.

The measured figures fall far below them. At $n=10$ the recurrence reaches $4653$ states and stores $9306$ transition lists — each reachable state is asked once with the last-row flag clear and once with it set — holding $177830$ outgoing entries built from $788908$ sweep nodes; at $n=12$ these become $35169$, $70338$, $2840330$ and $12792262$, so the states and lists grow by about $7.6$ and the entries and nodes by about $16$. On the machine used here, compiled at -O2 with the checkpoints skipped, the C++ needs about $0.07$ seconds at $n=10$ and $0.8$ at $n=12$, Java about $0.19$ and Python about $0.33$ at $n=10$ — single untuned runs, not machine-independent costs. All three print the same value.

Besides the C++ checkpoints, these checks were run for this article.

Direct enumeration of the ant moves, rejecting a repeated target or a two-cycle at once, gives $0,2,0,88,0$ for $n=1,\dots,5$.
The brute force of Section 5 reaches $n=6$, where the $6728$ matchings yield $207408$ edge-disjoint ordered pairs, so Section 2 is tested with no recurrence at all.
The degree-two search of Section 4 gives its $18$ covers and $\sum_C 2^{c(C)}=88$, the $6$ single-cycle covers matching the published Hamiltonian count; ^[7] the matching counts, recomputed by a one-layer recurrence, agree with the product formula and the sequence. ^[6]
A second transfer matrix, written independently, agrees with all three programs for every $n$ from $1$ to $12$, confirming $f(6)=207408$ and $f(8)=22902801416$. It enumerates all $3^n$ decision vectors of a row and validates each afterwards instead of pruning, and compares the matchings by intersecting their horizontal and vertical masks.

One latent defect: the cache key reserves bit $24$ for the last-row flag, but $i_2$ occupies the bits $n$ to $2n-1$, which reach bit $24$ once $n\ge 13$. The flag is OR-ed in, so it is then lost whenever that bit of $i_2$ is already set, and the last row silently reuses the ordinary-row list. The C++ range check and the default $n=10$ of the other two keep every published run clear of this, though Python's solve(n) does not itself enforce the limit; a flag field above both masks would. At $n=14$ the reuse happens $90945$ times without changing the printed value, since only the empty outgoing pair is read at the end and a filling reaching it places no downward edge.

Footnotes and References

The references supply the background facts used above; the matching reduction, the cycle-cover tally, the two forbidden combinations and the multiplicity-one claim are derived here.

Project Euler 393 — Migrating Ants. The official statement moves every ant of an n × n grid at once to a neighbouring square, counts only the moves with no two ants on one square and no two crossing one edge, gives f(4) = 88 and asks for f(10).
Wikipedia — Permutation. A permutation of a finite set decomposes into disjoint cycles, and one exchanging two elements while fixing the rest is a transposition, a single 2-cycle.
Wikipedia — Bipartite graph. Colouring the two parts differently gives every edge endpoints of differing colours, equivalently there is no odd cycle; grid graphs are named as bipartite.
Wikipedia — Perfect matching. A perfect matching is a set of edges with every vertex incident to exactly one of them, and it can exist only when the vertex count is even.
Wikipedia — Domino tiling. A domino tiling is a perfect matching of the grid graph, and the tilings of an m × n rectangle number the product over j and k of 4cos²(πj/(m + 1)) + 4cos²(πk/(n + 1)), found independently by Temperley and Fisher and by Kasteleyn in 1961.
OEIS A004003 — domino tilings of a 2n × 2n square. The sequence begins 1, 2, 36, 6728, 12988816, 258584046368, the matching counts used in Section 5.
OEIS A003763 — Hamiltonian cycles on a 2n × 2n grid. The count begins 1, 6, 1072, 4638576, so the 4 × 4 grid carries exactly 6, cross-checking the tally in Section 4.
cp-algorithms — Dynamic programming on broken profile. Domino tilings are counted with a mask recording the frontier row's profile, a row filled cell by cell: a covered cell is skipped, a vertical domino marks the cell in the next row's mask, a horizontal domino covers two free neighbours. The page gives no complexity bound, so the bounds above are our own.
C++ — Euler393.cpp. The linked immutable revision holds the Solver class, brute_force_small, the checkpoints, the range check 1 ≤ n ≤ 12 and the boost::multiprecision::cpp_int accumulator; the extra checks above are not packaged in it.

Problem 393 source code

C++

#include <cstdint>
#include <functional>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using BigInt = boost::multiprecision::cpp_int;

struct Options {
    int n = 10;
    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;
    }
    int parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(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, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1 && options.n <= 12;
}

enum class MoveType : std::uint8_t { Filled, Horizontal, Vertical };

struct Move {
    MoveType type;
    int occ_add;
    int out_add;
};

class Solver {
   public:
    explicit Solver(const int n) : n_(n), bit_mask_((1 << n) - 1) {}

    BigInt solve() {
        std::unordered_map<int, BigInt> dp;
        dp.reserve(1 << 12);
        dp.emplace(0, BigInt(1));

        for (int row = 0; row < n_; ++row) {
            const bool last_row = (row == n_ - 1);
            std::unordered_map<int, BigInt> next;
            next.reserve(dp.size() * 4 + 16);

            for (const auto& entry : dp) {
                const int state = entry.first;
                const BigInt& ways = entry.second;
                const int in1 = state & bit_mask_;
                const int in2 = (state >> n_) & bit_mask_;

                const auto& transitions = get_transitions(in1, in2, last_row);
                for (const auto& tr : transitions) {
                    next[tr.first] += ways * tr.second;
                }
            }

            dp.swap(next);
        }

        const auto it = dp.find(0);
        if (it == dp.end()) {
            return BigInt(0);
        }
        return it->second;
    }

   private:
    int n_;
    int bit_mask_;
    std::unordered_map<std::uint32_t, std::vector<std::pair<int, std::uint32_t>>> transition_cache_;

    std::vector<Move> options_for_cell(const int c, const int occ, const bool last_row) const {
        std::vector<Move> options;
        const int bit = 1 << c;
        if ((occ & bit) != 0) {
            options.push_back(Move{MoveType::Filled, 0, 0});
            return options;
        }

        if (c + 1 < n_ && (occ & (1 << (c + 1))) == 0) {
            options.push_back(Move{MoveType::Horizontal, bit | (1 << (c + 1)), 0});
        }
        if (!last_row) {
            options.push_back(Move{MoveType::Vertical, bit, bit});
        }

        return options;
    }

    const std::vector<std::pair<int, std::uint32_t>>& get_transitions(
        const int in1, const int in2, const bool last_row) {
        const std::uint32_t key = static_cast<std::uint32_t>(in1 | (in2 << n_) | (static_cast<int>(last_row) << 24));
        const auto cache_it = transition_cache_.find(key);
        if (cache_it != transition_cache_.end()) {
            return cache_it->second;
        }

        std::unordered_map<int, std::uint32_t> accum;
        accum.reserve(128);

        std::function<void(int, int, int, int, int)> dfs =
            [&](const int c, const int occ1, const int occ2, const int out1, const int out2) {
                if (c == n_) {
                    const int out_state = out1 | (out2 << n_);
                    ++accum[out_state];
                    return;
                }

                const std::vector<Move> ops1 = options_for_cell(c, occ1, last_row);
                const std::vector<Move> ops2 = options_for_cell(c, occ2, last_row);

                for (const Move& m1 : ops1) {
                    for (const Move& m2 : ops2) {
                        if (m1.type == MoveType::Horizontal && m2.type == MoveType::Horizontal) {
                            continue;  // Same undirected horizontal edge used by both matchings.
                        }
                        if (m1.type == MoveType::Vertical && m2.type == MoveType::Vertical) {
                            continue;  // Same undirected vertical edge used by both matchings.
                        }

                        const int next_occ1 = occ1 | m1.occ_add;
                        const int next_occ2 = occ2 | m2.occ_add;
                        const int next_out1 = out1 | m1.out_add;
                        const int next_out2 = out2 | m2.out_add;
                        dfs(c + 1, next_occ1, next_occ2, next_out1, next_out2);
                    }
                }
            };

        dfs(0, in1, in2, 0, 0);

        std::vector<std::pair<int, std::uint32_t>> packed;
        packed.reserve(accum.size());
        for (const auto& entry : accum) {
            packed.push_back(entry);
        }

        auto [it, inserted] = transition_cache_.emplace(key, std::move(packed));
        (void)inserted;
        return it->second;
    }
};

bool is_valid_assignment(const int n, const std::vector<int>& to) {
    const int total = n * n;
    std::vector<int> indeg(static_cast<std::size_t>(total), 0);
    auto rc = [n](const int id) {
        return std::pair<int, int>{id / n, id % n};
    };

    for (int u = 0; u < total; ++u) {
        const int v = to[static_cast<std::size_t>(u)];
        if (v < 0 || v >= total) {
            return false;
        }
        const auto [r1, c1] = rc(u);
        const auto [r2, c2] = rc(v);
        const int manhattan = std::abs(r1 - r2) + std::abs(c1 - c2);
        if (manhattan != 1) {
            return false;
        }
        ++indeg[static_cast<std::size_t>(v)];
        if (indeg[static_cast<std::size_t>(v)] > 1) {
            return false;
        }
    }

    for (int u = 0; u < total; ++u) {
        const int v = to[static_cast<std::size_t>(u)];
        if (to[static_cast<std::size_t>(v)] == u) {
            return false;  // Crossing the same edge in opposite directions.
        }
    }

    for (int v = 0; v < total; ++v) {
        if (indeg[static_cast<std::size_t>(v)] != 1) {
            return false;
        }
    }

    return true;
}

std::uint64_t brute_force_small(const int n) {
    const int total = n * n;
    std::vector<std::vector<int>> neighbors(static_cast<std::size_t>(total));
    auto id = [n](const int r, const int c) {
        return r * n + c;
    };

    for (int r = 0; r < n; ++r) {
        for (int c = 0; c < n; ++c) {
            const int u = id(r, c);
            if (r > 0) neighbors[static_cast<std::size_t>(u)].push_back(id(r - 1, c));
            if (r + 1 < n) neighbors[static_cast<std::size_t>(u)].push_back(id(r + 1, c));
            if (c > 0) neighbors[static_cast<std::size_t>(u)].push_back(id(r, c - 1));
            if (c + 1 < n) neighbors[static_cast<std::size_t>(u)].push_back(id(r, c + 1));
        }
    }

    std::vector<int> to(static_cast<std::size_t>(total), -1);
    std::uint64_t count = 0;

    std::function<void(int)> dfs = [&](const int u) {
        if (u == total) {
            if (is_valid_assignment(n, to)) {
                ++count;
            }
            return;
        }
        for (const int v : neighbors[static_cast<std::size_t>(u)]) {
            to[static_cast<std::size_t>(u)] = v;
            dfs(u + 1);
        }
    };

    dfs(0);
    return count;
}

bool run_checkpoints() {
    Solver s2(2);
    const BigInt f2 = s2.solve();
    if (f2 != 2) {
        std::cerr << "Checkpoint failed: f(2)\n";
        return false;
    }

    Solver s4(4);
    const BigInt f4 = s4.solve();
    if (f4 != 88) {
        std::cerr << "Checkpoint failed: f(4)\n";
        return false;
    }

    const std::uint64_t brute3 = brute_force_small(3);
    Solver s3(3);
    const BigInt f3 = s3.solve();
    if (f3 != brute3) {
        std::cerr << "Checkpoint failed: brute-force mismatch for n=3\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    Solver solver(options.n);
    std::cout << solver.solve() << '\n';
    return 0;
}

Python

def solve(n=10):
    bit_mask = (1 << n) - 1
    dp = {0: 1}
    transition_cache = {}
    
    def options_for_cell(c, occ, last_row):
        options = []
        bit = 1 << c
        if occ & bit:
            options.append((0, 0, 0))
            return options
            
        if c + 1 < n and (occ & (1 << (c + 1))) == 0:
            options.append((1, bit | (1 << (c + 1)), 0))
            
        if not last_row:
            options.append((2, bit, bit))
            
        return options
        
    def get_transitions(in1, in2, last_row):
        key = in1 | (in2 << n) | (int(last_row) << 24)
        if key in transition_cache:
            return transition_cache[key]
            
        accum = {}
        
        def dfs(c, occ1, occ2, out1, out2):
            if c == n:
                out_state = out1 | (out2 << n)
                accum[out_state] = accum.get(out_state, 0) + 1
                return
                
            ops1 = options_for_cell(c, occ1, last_row)
            ops2 = options_for_cell(c, occ2, last_row)
            
            for mt1, m1_occ, m1_out in ops1:
                for mt2, m2_occ, m2_out in ops2:
                    if mt1 == 1 and mt2 == 1: continue
                    if mt1 == 2 and mt2 == 2: continue
                    
                    dfs(c + 1, occ1 | m1_occ, occ2 | m2_occ, out1 | m1_out, out2 | m2_out)
                    
        dfs(0, in1, in2, 0, 0)
        items = list(accum.items())
        transition_cache[key] = items
        return items

    for row in range(n):
        last_row = (row == n - 1)
        next_dp = {}
        for state, ways in dp.items():
            in1 = state & bit_mask
            in2 = (state >> n) & bit_mask
            transitions = get_transitions(in1, in2, last_row)
            
            for out_state, ways_add in transitions:
                next_dp[out_state] = next_dp.get(out_state, 0) + ways * ways_add
                
        dp = next_dp
        
    return str(dp.get(0, 0))

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler393 {

    enum MoveType {
        FILLED, HORIZONTAL, VERTICAL
    }

    static class Move {
        MoveType type;
        int occAdd;
        int outAdd;

        Move(MoveType t, int occ, int out) {
            this.type = t;
            this.occAdd = occ;
            this.outAdd = out;
        }
    }

    static class Transition {
        int outState;
        int ways;

        Transition(int s, int w) {
            outState = s;
            ways = w;
        }
    }

    static int n = 10;
    static int bitMask = (1 << n) - 1;
    static Map<Integer, List<Transition>> transitionCache = new HashMap<>();

    static List<Move> optionsForCell(int c, int occ, boolean lastRow) {
        List<Move> options = new ArrayList<>(3);
        int bit = 1 << c;
        if ((occ & bit) != 0) {
            options.add(new Move(MoveType.FILLED, 0, 0));
            return options;
        }

        if (c + 1 < n && (occ & (1 << (c + 1))) == 0) {
            options.add(new Move(MoveType.HORIZONTAL, bit | (1 << (c + 1)), 0));
        }
        if (!lastRow) {
            options.add(new Move(MoveType.VERTICAL, bit, bit));
        }

        return options;
    }

    static void dfs(int c, int occ1, int occ2, int out1, int out2, boolean lastRow, Map<Integer, Integer> accum) {
        if (c == n) {
            int outState = out1 | (out2 << n);
            accum.put(outState, accum.getOrDefault(outState, 0) + 1);
            return;
        }

        List<Move> ops1 = optionsForCell(c, occ1, lastRow);
        List<Move> ops2 = optionsForCell(c, occ2, lastRow);

        for (Move m1 : ops1) {
            for (Move m2 : ops2) {
                if (m1.type == MoveType.HORIZONTAL && m2.type == MoveType.HORIZONTAL)
                    continue;
                if (m1.type == MoveType.VERTICAL && m2.type == MoveType.VERTICAL)
                    continue;

                dfs(c + 1, occ1 | m1.occAdd, occ2 | m2.occAdd, out1 | m1.outAdd, out2 | m2.outAdd, lastRow, accum);
            }
        }
    }

    static List<Transition> getTransitions(int in1, int in2, boolean lastRow) {
        int key = in1 | (in2 << n) | ((lastRow ? 1 : 0) << 24);
        List<Transition> cached = transitionCache.get(key);
        if (cached != null)
            return cached;

        Map<Integer, Integer> accum = new HashMap<>();
        dfs(0, in1, in2, 0, 0, lastRow, accum);

        List<Transition> packed = new ArrayList<>(accum.size());
        for (Map.Entry<Integer, Integer> e : accum.entrySet()) {
            packed.add(new Transition(e.getKey(), e.getValue()));
        }

        transitionCache.put(key, packed);
        return packed;
    }

    static String solve() {
        Map<Integer, BigInteger> dp = new HashMap<>();
        dp.put(0, BigInteger.ONE);

        for (int row = 0; row < n; ++row) {
            boolean lastRow = (row == n - 1);
            Map<Integer, BigInteger> nextDp = new HashMap<>();

            for (Map.Entry<Integer, BigInteger> entry : dp.entrySet()) {
                int state = entry.getKey();
                BigInteger ways = entry.getValue();

                int in1 = state & bitMask;
                int in2 = (state >> n) & bitMask;

                List<Transition> transitions = getTransitions(in1, in2, lastRow);
                for (Transition tr : transitions) {
                    BigInteger add = ways.multiply(BigInteger.valueOf(tr.ways));
                    BigInteger current = nextDp.getOrDefault(tr.outState, BigInteger.ZERO);
                    nextDp.put(tr.outState, current.add(add));
                }
            }

            dp = nextDp;
        }

        BigInteger ans = dp.getOrDefault(0, BigInteger.ZERO);
        return ans.toString();
    }

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