Problem 237: Tours on a $4 \times N$ Playing Board

Project Euler Problem 237 Solution

EulerSolve provides an optimized solution for Project Euler Problem 237, Tours on a $4 \times N$ Playing Board, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $T(N)$ be the number of tours on a $4\times N$ board that start at the upper-left square, end at the lower-left square, and visit every square exactly once. In graph-theoretic language, we are counting Hamiltonian paths on the $4\times N$ grid graph with those two prescribed endpoints. The Project Euler instance asks for $T(10^{12}) \bmod 10^8$. A direct search is hopeless, and even stepping through a recurrence one column at a time would still take $10^{12}$ iterations. The solution therefore needs two ingredients: a finite-width transfer argument to get a fixed recurrence, and fast matrix exponentiation to jump to the trillionth column length. Mathematical Approach The decisive simplification is that the board is extremely long but only 4 cells tall. Fixed width means a left-to-right sweep can leave only finitely many possible boundary patterns, so the counting problem becomes a small linear dynamical system. Frontier states on a strip of width 4 Sweep the board from left to right and imagine that the first $n$ columns have already been processed. Every square strictly inside that processed region must already have its final degree in the Hamiltonian path, namely degree 2, except for the two prescribed endpoints in the first column. Any unfinished information can only live on the cut between columns $n$ and $n+1$....

Detailed mathematical approach

Problem Summary

Number the rows of a $4\times n$ board $1$ to $4$ from the top and its columns $1$ to $n$ from the left. A tour starts in the top left square, moves one square up, down, left or right at a time, visits each of the $4n$ squares exactly once and ends in the bottom left square. $T(n)$ is the number of tours, that is, of Hamiltonian paths of the $4\times n$ grid graph between its two left corners. ^[2] The statement gives $T(10)=2329$ and asks for $T(10^{12})\bmod 10^8$. ^[1]

The whole solution in one sentence: cutting the first column off a tour leaves one or two paths of a few kinds on a smaller board, their counts obey five simple equations that reduce to a recurrence of order four, and a $4\times4$ matrix power computed by repeated squaring modulo $10^8$ jumps straight to $n=10^{12}$.

Mathematical Approach

We name a square by its row and column digits: 11 is the start, 41 the end and 23 the square in row 2 and column 3.

1. Small boards

Every square of a tour except the two ends is entered once and left once, so it has two neighbours along the tour; the ends have one. On one column the tour walks down, 11 → 21 → 31 → 41, so $T(1)=1$. On two columns the corner 12 has only the neighbours 11 and 22, so the tour runs 11 → 12 → 22, and likewise 32 → 42 → 41. The ends have now used their only step, so 21 is joined to 22 and 31, and the tour 11 → 12 → 22 → 21 → 31 → 32 → 42 → 41 is forced: $T(2)=1$. A search finds $T(3)=4$ and $T(4)=8$, but no search reaches $n=10^{12}$.

2. Covers by paths

Cutting the first column off a tour leaves pieces that are not tours, so we count more general objects. A cover of a board is a set of paths, each stepping between squares that share a side, that together visit every square exactly once. All covers below have their path ends in the first column. On a board with $m$ columns let

$T(m)$ count the covers by one path joining rows 1 and 4, the tours;
$B(m)$ count those by one path joining rows 1 and 2; turning the board upside down shows that the covers by one path joining rows 3 and 4 are just as many;
$C(m)$ count those by one path joining rows 2 and 3;
$D(m)$, $E(m)$ and $X(m)$ count those by two paths that join the rows in pairs: 1 with 2 and 3 with 4 for $D$, 1 with 4 and 2 with 3 for $E$, 1 with 3 and 2 with 4 for $X$.

For one column $T(1)=D(1)=1$, the cover for $D$ being the paths 11–21 and 31–41, and the other counts are $0$: a path through a whole column runs from row 1 to row 4, and for $E$ and $X$ the only neighbour of 11 is 21, an end of the other path.

3. Cutting off the first column

Let $m\ge2$. A square outside column 1 is not an end, so it has two neighbours along its path. Removing column 1 from a cover of one of the six kinds takes one neighbour from each square of column 2 joined to column 1, and none from the others. What remains is a cover of columns $2,\dots,m$ without cycles whose ends are exactly those joined squares: again a cover with all ends in its first column, on $m-1$ columns. Conversely, the piece in column 1 and such a smaller cover form a cover of the original kind exactly when the pieces link up into the required paths without closing a cycle.

A square of column 1 needs one neighbour if it is an end and two otherwise, among the squares above, below and to its right. Once we choose which of the vertical edges 11–21, 21–31 and 31–41 are used, a square takes the edge to its right exactly when it is one neighbour short; a choice fails if a square gets too many vertical edges or would need two edges to the right.

4. The tours

In a tour 11 and 41 need one neighbour, 21 and 31 two. The eight choices give:

all three edges: the column is the whole tour, possible only for $m=1$;
11–21 and 31–41: 21 and 31 go right, and the rest is a path joining rows 2 and 3, in $C(m-1)$ ways;
11–21 and 21–31: 31 and 41 go right, in $B(m-1)$ ways, and 21–31 with 31–41 likewise gives $B(m-1)$;
only 21–31: all four squares go right (below, ⇝ marks a path of the smaller cover). Smaller paths joining 1 with 4 and 2 with 3 would close a cycle through 21–31, while paths joining 1 with 2 and 3 with 4 link up into the tour 11 → 12 ⇝ 22 → 21 → 31 → 32 ⇝ 42 → 41, and paths joining 1 with 3 and 2 with 4 link up as well: $D(m-1)+X(m-1)$ ways;
only 11–21, only 31–41 or none: 31 or 21 would need two edges to the right.

Hence, for $m\ge2$,

$$\begin{aligned}T(m)={}&C(m-1)+2B(m-1)\\ &+D(m-1)+X(m-1).\end{aligned}$$

5. One path with other ends

For $B$ the ends are 11 and 21. The square 41 has only the neighbours 31 and 42 and uses both, and 11–21 is unused, since the path would be just 11–21. With 21–31, only 11 and 41 go right and the rest is a tour; without it, all four squares go right, and the smaller board must join 1 with 4 and 2 with 3, or 1 with 3 and 2 with 4, since joining 1 with 2 would close a cycle through 31–41. For $C$ the ends are 21 and 31, so 11 and 41 need two neighbours each; this forces 11–21, 31–41 and the edges to the right in rows 1 and 4, after which both ends are complete and the rest is a tour. Hence

$$\begin{aligned}B(m)={}&T(m-1)+E(m-1)\\ &+X(m-1),\\ C(m)={}&T(m-1).\end{aligned}$$

6. Two paths

For $D$, $E$ and $X$ all four squares of column 1 are ends, and no vertical edge may join two different paths. For $D$ this excludes 21–31. Using both 11–21 and 31–41 completes both paths, possible only for $m=1$; using one of them leaves the other path to cover the smaller board alone, in $B(m-1)$ ways each; using none sends all four squares right, and the smaller cover is again of kind $D$. For $E$, the edge 21–31 makes the inner path 21–31 and leaves a tour, and without it the smaller cover is again of kind $E$. For $X$ every vertical edge is excluded, so the smaller cover is again of kind $X$. Hence

$$\begin{aligned}D(m)&=2B(m-1)+D(m-1),\\ E(m)&=T(m-1)+E(m-1),\\ X(m)&=X(m-1).\end{aligned}$$

Since $X(1)=0$, induction gives $X(m)=0$ for all $m$: no cover joins rows 1 with 3 and 2 with 4, as such paths would have to cross.

7. Worked example: up to six columns

From the values for $m=1$, the equations give these values of $(T,B,C,D,E)$:

$m=1$: $(1,0,0,1,0)$; $m=2$: $(1,1,1,1,1)$; $m=3$: $(4,2,1,3,2)$;
$m=4$: $(8,6,4,7,6)$; $m=5$: $(23,14,8,19,14)$;
$m=6$: $(55,37,23,47,37)$.

For instance $T(3)=C(2)+2B(2)+D(2)=4$. In the order of the cases of Section 4, the four tours on three columns are

11 → 21 → 22 → 12 → 13 → 23 → 33 → 43 → 42 → 32 → 31 → 41;
11 → 21 → 31 → 32 → 22 → 12 → 13 → 23 → 33 → 43 → 42 → 41;
11 → 12 → 13 → 23 → 33 → 43 → 42 → 32 → 22 → 21 → 31 → 41;
11 → 12 → 13 → 23 → 22 → 21 → 31 → 32 → 33 → 43 → 42 → 41.

8. The recurrence

With $X=0$ the equations read, for $m\ge2$,

$$\begin{aligned}T(m)&=C(m-1)+2B(m-1)\\ &\quad+D(m-1),\\ B(m)&=T(m-1)+E(m-1),\\ C(m)&=T(m-1),\\ D(m)&=2B(m-1)+D(m-1),\\ E(m)&=T(m-1)+E(m-1).\end{aligned}$$

We eliminate the helper counts in four steps. First, as $E(1)=0$, the last equation gives $E(m)=T(1)+\dots+T(m-1)$. Then $B(m)=T(m-1)+E(m-1)$${}=E(m)$ for $m\ge2$, and $B(1)=E(1)$, so $B=E$ and $B(m)-B(m-1)=T(m-1)$ for $m\ge2$. Second, by the equation for $D$ the first equation reads $T(m)=C(m-1)+D(m)$, and with $C(m-1)=T(m-2)$ this gives $T(m)-T(m-2)=D(m)$ for $m\ge3$. Third, subtracting this identity at $m-1$ from the one at $m$ and using the equation for $D$ gives, for $m\ge4$,

$$\begin{aligned}T(m)&{}-T(m-1)-T(m-2)\\ &{}+T(m-3)=2B(m-1).\end{aligned}$$

Fourth, doing the same with this identity and using the first step gives, for $m\ge5$,

$$\begin{aligned}T(m)&{}-2T(m-1)+2T(m-3)\\ &{}-T(m-4)=2T(m-2).\end{aligned}$$

Rearranged, this is the recurrence of the solution: for $m\ge5$,

$$\begin{aligned}T(m)={}&2T(m-1)+2T(m-2)\\ &-2T(m-3)+T(m-4).\end{aligned}$$

From $T(1),\dots,T(4)=1,1,4,8$ it gives $T(5)=2\cdot8+2\cdot4$${}-2\cdot1+1=23$, then $55$, $144$, $360$, $921$ and $T(10)=2329$, as in the statement; the coefficient $-2$ comes from taking differences. The OEIS lists the same recurrence and values. ^[3]

A by-product: adding the edge 11–21 to a path counted by $B(m)$ closes a Hamiltonian cycle, and every Hamiltonian cycle of the board uses both edges at the corner 11, so $B(m)=T(1)+\dots+T(m-1)$ also counts these cycles, $1,2,6,14,37,\dots$ for $m\ge2$, as the OEIS records. ^[4]

9. Transfer matrices and growth

The equations say $v(m)=M\,v(m-1)$ for the column vector $v(m)$ of $T,B,C,D,E$, with $v(1)=(1,0,0,1,0)$ and

$$M=\begin{pmatrix}0&2&1&1&0\\1&0&0&0&1\\1&0&0&0&0\\0&2&0&1&0\\1&0&0&0&1\end{pmatrix}.$$

For example, the first row of $M$ turns the counts $(4,2,1,3,2)$ for three columns into $T(4)=2\cdot2+1+3=8$. This is the transfer-matrix method, in which entries of matrix powers count walks, here walks through the five kinds, one column per step; it makes the generating function $G(x)=\sum_{n\ge1}T(n)\,x^n$ rational with denominator $\det(I-xM)$, ^[5] here $1-2x-2x^2+2x^3-x^4$. Indeed, multiplying $G(x)$ by this polynomial cancels every coefficient from $x^5$ on, by the recurrence, and the first four leave $x-x^2$:

$$G(x)=\frac{x(1-x)}{1-2x-2x^2+2x^3-x^4}.$$

By partial fractions, $T(n)$ is a combination of the $n$-th powers of the four distinct roots of $x^4-2x^3-2x^2+2x-1$. ^[6] Let $r=2.5386\ldots$ be the largest; the others have absolute values below $1.28$, so $T(n)=c\,r^n+O(1.28^n)$ with a constant $c=0.2098\ldots$, and $T(10^{12})$ has about $4.05\cdot10^{11}$ digits. We compute only its remainder.

10. Jumping ahead with a matrix power

Let $X_n$ be the column vector of $T(n)$, $T(n-1)$, $T(n-2)$, $T(n-3)$. Then $X_n=A\,X_{n-1}$ for $n\ge5$, with

$$A=\begin{pmatrix}2&2&-2&1\\1&0&0&0\\0&1&0&0\\0&0&1&0\end{pmatrix},$$

whose first row applies the recurrence while the other rows shift the values down, ^[6] so $X_n=A^{n-4}X_4$ with $X_4=(8,4,1,1)$. Exponentiation by squaring ^[7] computes $A^e$ for $e=n-4$: a product $R$ starts as the identity matrix and a power $P$ as $A$; while $e\gt0$, $R$ is multiplied by $P$ if the last binary digit of $e$ is $1$, then $P$ is squared and that digit dropped. As $P$ runs through $A,A^2,A^4,\dots$, $R$ ends as the product of the $A^{2^j}$ over the digits $1$ of $e$, which is $A^e$.

Worked example, $n=10$: here $e=6$, which is 110 in binary. The digit $0$ only squares $P$ to $A^2$; the next digit gives $R=A^2$ and $P=A^4$; the last gives $R=A^6$, whose first row is $(236,124,-147,92)$. Applied to $X_4=(8,4,1,1)$ it gives

$$\begin{aligned}T(10)&=236\cdot8+124\cdot4\\ &\quad-147+92=2329.\end{aligned}$$

For $n=10^{12}$ the exponent $e=999\,999\,999\,996$ has $40$ binary digits, $22$ of them equal to $1$, so $40$ squarings and $22$ products replace about $10^{12}$ steps of the recurrence.

11. Arithmetic modulo $10^8$

Congruences modulo $10^8$ can be added and multiplied, ^[8] so every entry may be reduced as soon as it is formed, and the first entry at the end is $T(n)\bmod 10^8$. The entry $-2$ is stored as $10^8-2$, so all stored numbers lie in $0,\dots,10^8-1$. No division occurs, so it does not matter that $10^8$ is not a prime. A product of two stored numbers plus a stored number is below $10^{16}$, far below the signed 64-bit limit $2^{63}\approx9.2\cdot10^{18}$.

How the Code Works

The C++ program ^[9] starts in main with an Options record holding $n=10^{12}$, the modulus $10^8$ and run_checkpoints = true. parse_arguments accepts --n=, --modulo= and --skip-checkpoints; it returns false on any other argument, on $n\lt1$ or on a modulus below $2$, and main then exits with code 1. Its digit parsers do not detect overflow. Unless --skip-checkpoints is given, run_checkpoints compares solve(n, 1000000000) with brute_tours(n) for $n=1,\dots,4$ and requires solve(10, 1000000000) to be $2329$; if a check fails, main exits with code 2. brute_tours does not use Sections 3 to 8: its table dp[mask][v] counts the paths from the start that visit exactly the squares of the bit mask mask and end at v, so the full mask at the bottom left square counts the tours; it returns $0$ for boards with more than $20$ squares.

solve(n, mod) returns stored values for $n\le4$. Otherwise it fills step with $A$, storing $-2$ as mod_norm(-2, mod), sets acc to the identity and runs the loop of Section 10 on power = n - 4: acc = mat_mul(acc, step, mod) when the lowest bit is set, then step is squared and power shifted right. mat_vec_mul applies acc to base, which is $X_4$, and main prints the first entry. mat_mul multiplies $4\times4$ arrays of long long, skips zero entries of the left factor and reduces after every multiply-add with mod_norm, which returns a remainder from $0$ to the modulus minus $1$ even for negative input.

The Python version has the same functions on lists of unbounded integers, without argument parsing or checks; its % is never negative for a positive modulus, so the correction in mod_norm never fires, and solve returns a string for $n\ge5$ but an integer for $n\le4$. The Java version has modNorm, matMul and matVecMul on long arrays and a solve() with $n$ and the modulus fixed, also without checks.

Complexity and Verification

For $n=10^{12}$ the loop runs $40$ times and calls mat_mul $62$ times, with at most $64$ multiply-adds each (the last squaring is not needed). Counting operations on numbers of machine size, the time is $O(\log n)$ and the extra memory $O(1)$. The C++ checks cost more than the computation, as the table of brute_tours for $n=4$ has $2^{16}\cdot16=1\,048\,576$ entries; still, in our runs each program finished within a few hundredths of a second. Since solve returns stored values for $n\le4$, the brute-force comparison checks those values; the matrix power is tested by $T(10)=2329$. For this article we also verified:

a depth-first search over all covers, independent of Sections 3 to 8, gives the six counts for $m\le16$; they satisfy all equations of Sections 4 to 6, with $X=0$ and $T(10)=2329$;
the recurrence reproduces the $30$ values listed in ^[3], and the sums $T(1)+\dots+T(m-1)$ the $32$ cycle counts listed in ^[4];
with --n= and --modulo=, the C++ program agrees with step-by-step evaluation of the recurrence in $13$ cases, among them $n=10^6$ with the moduli $10^8$ and $10^9+7$;
powers of the matrix $M$ of Section 9 modulo $10^8$, computed separately, agree with step-by-step evaluation for all $n\le3000$ and give the same residue for $n=10^{12}$ as the three programs.

Footnotes and References

The references supply background facts; the cover counts, their equations and the elimination are derived in this article.

Project Euler 237 — Tours on a 4 × N Playing Board. The official statement defines T(n) for tours from the top left to the bottom left corner of a 4 × n board, gives T(10) = 2329 and asks for T(10¹²) modulo 10⁸.
Wikipedia — Hamiltonian path. A Hamiltonian path visits each vertex of a graph exactly once, and a Hamiltonian cycle is a cycle that does so.
OEIS A181688. The numbers T(n), called there maximal self-avoiding walks from the NW to the SW corner of a 4-by-n grid: 30 terms, the recurrence of Section 8 for n > 4 and the generating function of Section 9.
OEIS A006864. The numbers of Hamiltonian cycles of the 4 × n grid, 0, 1, 2, 6, 14, 37, …, with the formula a(n) = A181688(1) + … + A181688(n − 1) of Section 8.
R. P. Stanley, Enumerative Combinatorics, Volume 1, second edition, Section 4.7, The Transfer-matrix Method: powers of the adjacency matrix A of a weighted directed graph count its walks (Theorem 4.7.1), and their generating functions are rational with denominator det(I − λA) (Theorem 4.7.2).
Wikipedia — Linear recurrence with constant coefficients. A recurrence of order n becomes a first-order matrix equation, the matrix having the coefficients in its first row and a single 1 in each other row; its generating function is rational with denominator 1 − a₁x − … − aₙxⁿ, and partial fractions turn it into a combination of n-th powers of the roots.
Wikipedia — Exponentiation by squaring. Each step removes the last binary digit of the exponent, so about log₂ n squarings and at most as many multiplications suffice; the method is often used for matrix powers and modular powers.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition and multiplication: congruent summands give congruent sums, and congruent factors congruent products.
C++ — Euler237.cpp. The linked immutable C++ revision contains mat_mul, solve, brute_tours and run_checkpoints; it is the implementation record for this article.

Problem 237 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;

struct Options {
    u64 n = 1000000000000ULL;
    int modulo = 100000000;
    bool run_checkpoints = 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 = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
    }
    value = parsed;
    return 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 c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '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_u64_after_prefix(arg, "--n=", options.n) ||
            parse_int_after_prefix(arg, "--modulo=", options.modulo)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1ULL && options.modulo >= 2;
}

i64 mod_norm(const i64 x, const int mod) {
    i64 v = x % mod;
    if (v < 0) {
        v += mod;
    }
    return v;
}

using Mat4 = std::array<std::array<i64, 4>, 4>;
using Vec4 = std::array<i64, 4>;

Mat4 mat_mul(const Mat4& a, const Mat4& b, const int mod) {
    Mat4 c{};
    for (int i = 0; i < 4; ++i) {
        for (int k = 0; k < 4; ++k) {
            if (a[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)] == 0) {
                continue;
            }
            for (int j = 0; j < 4; ++j) {
                c[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] =
                    mod_norm(c[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] +
                                 a[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)] *
                                     b[static_cast<std::size_t>(k)][static_cast<std::size_t>(j)],
                             mod);
            }
        }
    }
    return c;
}

Vec4 mat_vec_mul(const Mat4& a, const Vec4& x, const int mod) {
    Vec4 y{};
    for (int i = 0; i < 4; ++i) {
        i64 sum = 0;
        for (int j = 0; j < 4; ++j) {
            sum += a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] *
                   x[static_cast<std::size_t>(j)];
        }
        y[static_cast<std::size_t>(i)] = mod_norm(sum, mod);
    }
    return y;
}

i64 solve(const u64 n, const int mod) {
    // T(1)=1, T(2)=1, T(3)=4, T(4)=8
    if (n == 1ULL || n == 2ULL) {
        return 1 % mod;
    }
    if (n == 3ULL) {
        return 4 % mod;
    }
    if (n == 4ULL) {
        return 8 % mod;
    }

    Mat4 step{};
    step[0] = {2, 2, mod_norm(-2, mod), 1};
    step[1] = {1, 0, 0, 0};
    step[2] = {0, 1, 0, 0};
    step[3] = {0, 0, 1, 0};

    Mat4 acc{};
    for (int i = 0; i < 4; ++i) {
        acc[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)] = 1;
    }

    u64 power = n - 4ULL;
    while (power > 0ULL) {
        if ((power & 1ULL) != 0ULL) {
            acc = mat_mul(acc, step, mod);
        }
        step = mat_mul(step, step, mod);
        power >>= 1U;
    }

    const Vec4 base{8, 4, 1, 1};  // [T4, T3, T2, T1]
    const Vec4 out = mat_vec_mul(acc, base, mod);
    return out[0];
}

u64 brute_tours(int n) {
    const int rows = 4;
    const int cols = n;
    const int total = rows * cols;
    if (total > 20) {
        return 0;
    }

    const auto vid = [cols](const int r, const int c) { return r * cols + c; };
    std::vector<std::vector<int>> g(static_cast<std::size_t>(total));
    for (int r = 0; r < rows; ++r) {
        for (int c = 0; c < cols; ++c) {
            const int v = vid(r, c);
            if (r > 0) {
                g[static_cast<std::size_t>(v)].push_back(vid(r - 1, c));
            }
            if (r + 1 < rows) {
                g[static_cast<std::size_t>(v)].push_back(vid(r + 1, c));
            }
            if (c > 0) {
                g[static_cast<std::size_t>(v)].push_back(vid(r, c - 1));
            }
            if (c + 1 < cols) {
                g[static_cast<std::size_t>(v)].push_back(vid(r, c + 1));
            }
        }
    }

    const int start = vid(0, 0);
    const int target = vid(3, 0);
    const int max_mask = 1 << total;
    std::vector<std::vector<u64>> dp(static_cast<std::size_t>(max_mask),
                                     std::vector<u64>(static_cast<std::size_t>(total), 0ULL));
    dp[1 << start][static_cast<std::size_t>(start)] = 1ULL;

    for (int mask = 0; mask < max_mask; ++mask) {
        for (int v = 0; v < total; ++v) {
            const u64 ways = dp[static_cast<std::size_t>(mask)][static_cast<std::size_t>(v)];
            if (ways == 0) {
                continue;
            }
            for (int to : g[static_cast<std::size_t>(v)]) {
                if ((mask & (1 << to)) != 0) {
                    continue;
                }
                const int nmask = mask | (1 << to);
                dp[static_cast<std::size_t>(nmask)][static_cast<std::size_t>(to)] += ways;
            }
        }
    }
    return dp[static_cast<std::size_t>(max_mask - 1)][static_cast<std::size_t>(target)];
}

bool run_checkpoints() {
    for (int n = 1; n <= 4; ++n) {
        if (solve(static_cast<u64>(n), 1000000000) != static_cast<i64>(brute_tours(n))) {
            std::cerr << "Checkpoint failed for brute cross-check at n=" << n << '\n';
            return false;
        }
    }
    if (solve(10, 1000000000) != 2329) {
        std::cerr << "Checkpoint failed for T(10)=2329 sample" << '\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;
    }
    std::cout << solve(options.n, options.modulo) << '\n';
    return 0;
}

Python

def mod_norm(x, mod):
    v = x % mod
    if v < 0:
        v += mod
    return v

def mat_mul(a, b, mod):
    c = [[0] * 4 for _ in range(4)]
    for i in range(4):
        for k in range(4):
            if a[i][k] == 0:
                continue
            for j in range(4):
                c[i][j] = mod_norm(c[i][j] + a[i][k] * b[k][j], mod)
    return c

def mat_vec_mul(a, x, mod):
    y = [0] * 4
    for i in range(4):
        sum_val = 0
        for j in range(4):
            sum_val += a[i][j] * x[j]
        y[i] = mod_norm(sum_val, mod)
    return y

def solve(n=10**12, mod=10**8):
    if n == 1 or n == 2: return 1 % mod
    if n == 3: return 4 % mod
    if n == 4: return 8 % mod

    step = [
        [2, 2, mod_norm(-2, mod), 1],
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 1, 0]
    ]

    acc = [
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ]

    power = n - 4
    while power > 0:
        if power % 2 == 1:
            acc = mat_mul(acc, step, mod)
        step = mat_mul(step, step, mod)
        power //= 2

    base = [8, 4, 1, 1]
    out = mat_vec_mul(acc, base, mod)
    return str(out[0])

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

Java

public class Euler237 {
    static long modNorm(long x, int mod) {
        long v = x % mod;
        if (v < 0) {
            v += mod;
        }
        return v;
    }

    static long[][] matMul(long[][] a, long[][] b, int mod) {
        long[][] c = new long[4][4];
        for (int i = 0; i < 4; ++i) {
            for (int k = 0; k < 4; ++k) {
                if (a[i][k] == 0)
                    continue;
                for (int j = 0; j < 4; ++j) {
                    c[i][j] = modNorm(c[i][j] + a[i][k] * b[k][j], mod);
                }
            }
        }
        return c;
    }

    static long[] matVecMul(long[][] a, long[] x, int mod) {
        long[] y = new long[4];
        for (int i = 0; i < 4; ++i) {
            long sum = 0;
            for (int j = 0; j < 4; ++j) {
                sum += a[i][j] * x[j];
            }
            y[i] = modNorm(sum, mod);
        }
        return y;
    }

    public static String solve() {
        long n = 1_000_000_000_000L;
        int mod = 100_000_000;

        if (n == 1 || n == 2)
            return String.valueOf(1 % mod);
        if (n == 3)
            return String.valueOf(4 % mod);
        if (n == 4)
            return String.valueOf(8 % mod);

        long[][] step = {
                { 2, 2, modNorm(-2, mod), 1 },
                { 1, 0, 0, 0 },
                { 0, 1, 0, 0 },
                { 0, 0, 1, 0 }
        };

        long[][] acc = {
                { 1, 0, 0, 0 },
                { 0, 1, 0, 0 },
                { 0, 0, 1, 0 },
                { 0, 0, 0, 1 }
        };

        long power = n - 4;
        while (power > 0) {
            if (power % 2 == 1) {
                acc = matMul(acc, step, mod);
            }
            step = matMul(step, step, mod);
            power /= 2;
        }

        long[] base = { 8, 4, 1, 1 };
        long[] out = matVecMul(acc, base, mod);
        return String.valueOf(out[0]);
    }

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