Problem 416: A Frog's Trip

Project Euler Problem 416 Solution

EulerSolve provides an optimized solution for Project Euler Problem 416, A Frog's Trip, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $p=2m$. We count ordered $p$-tuples of frog legs, where each leg starts at node $1$, ends at node $n$, and every jump has length $1$, $2$, or $3$. Among the internal nodes $2,3,\dots,n-1$, at most one node may be missed by every leg. The implementations return the count modulo $10^9$. Mathematical Approach Step 1: Encode One Leg as a Binary Word Look only at the internal nodes $2,\dots,n-1$. For one leg, write a binary word of length $n-2$: put $1$ when that internal node is visited and $0$ when it is skipped. This converts the path problem into a word problem. Because jumps are only $1$, $2$, or $3$, a leg can skip at most two consecutive internal nodes. Therefore valid legs are exactly the binary words with no substring $000$. So the problem is not really about geometric paths anymore. It is about $p=2m$ synchronized binary words of length $n-2$, each avoiding three consecutive zeros, and we want the union of their $1$-positions to miss at most one internal node. Step 2: A Three-State Automaton for One Leg Process the internal nodes from left to right. After some prefix has been processed, a single leg is determined by the number of consecutive processed internal nodes that were skipped since its last visit. That number can only be $0$, $1$, or $2$....

Detailed mathematical approach

Problem Summary

A row of $n$ squares holds a frog in the leftmost square. On the outward trip it jumps one, two or three squares to the right until it stands on the rightmost square, on the homeward trip the same to the left until it is back at the start, never jumping outside the row. The frog repeats this round trip $m$ times, and $F(m,n)$ counts the travels leaving at most one square unvisited. The statement gives $F(1,3)=4$, $F(1,4)=15$, $F(1,5)=46$, $F(2,3)=16$ and $F(2,100)\bmod 10^9$${}=429619151$, and asks for the last nine digits of $F(10,10^{12})$. ^[1]

The whole solution in one sentence: a round trip is really two independent one-way trips, so we sweep the squares from left to right with a matrix that remembers nothing about the $2m$ trips except how many have just skipped zero, one or two squares in a row, attach a marker $y$ with $y^2=0$ to every square that all of them skip, and cross a row of $10^{12}$ squares by repeated squaring.

Mathematical Approach

The whole method is one transfer matrix of modest size raised to a huge power, with a single extra coefficient to handle the words "at most one".

1. A round trip is two one-way trips

Read the homeward trip backwards: its jumps to the left become jumps of the same lengths to the right, its ends swap, and the squares it touches do not change, so this is a bijection between homeward and outward trips. The $m$ round trips are independent of one another, and so are the two directions inside one round trip, so a travel plan is exactly an ordered tuple of

$$p=2m$$

outward trips, each of which we call a leg. The two checkpoints with $n=3$ confirm this at once: for $n=3$ there are two legs, $1\to3$ and $1\to2\to3$, and only the middle square could stay unvisited, so every plan is acceptable and $F(1,3)=2^2=4$, $F(2,3)=2^4=16$.

Squares $1$ and $n$ are touched by every leg, so only the $n-2$ interior squares can remain unvisited, and such a square is unvisited precisely when all $p$ legs jump over it. We must count the $p$-tuples in which at most one interior square is skipped by every leg.

2. What one leg looks like

A leg is determined by the squares it lands on, the jumps being the differences of consecutive landings. Record it as a word of length $n-2$ over the letters "land" and "skip", one per interior square. Consecutive landings differ by at most three, so at most two interior squares lie between them, and the same bound holds before the first landing and after the last, each being a single jump. A leg is therefore legal exactly when its word has no three consecutive skips, and every such word comes from one leg.

Equivalently a leg is a composition of $n-1$ into parts $1$, $2$ and $3$, the parts being the jump lengths. ^[2] Either count gives, for n = 2, 3, ..., 11, the leg counts 1, 2, 4, 7, 13, 24, 44, 81, 149, 274: the tribonacci sequence, each term the sum of the previous three. ^[3]

So $n=5$ has $7$ legs and $m=1$ gives exactly $49$ ordered plans; $F(1,5)=46$ says three of them fail, each leaving two interior squares unvisited.

3. Three states per leg

Sweep the interior squares from left to right. The only thing about a leg that can influence its future is how many interior squares it has skipped since its last landing, counting square $1$ as a landing. Call that number the state; by Section 2 it is $0$, $1$ or $2$, starting at $0$. At the next interior square a leg in state $2$ must land, a third skip being illegal, and returns to state $0$; a leg in state $0$ or $1$ may land, returning to state $0$, or skip, moving up one.

Nothing has to be checked at the end: after the last interior square a leg in state $r$ still has to reach square $n$, a jump of length $r+1\le3$, which is always legal. Every state is accepting, with exactly one completion.

4. Only the head count matters

The $p$ legs obey identical rules, never interact, and have futures depending on nothing but their state. After $t$ interior squares the whole plan therefore compresses to the triple $(a,b,c)$ with $a+b+c=p$, counting the legs in state $2$, state $1$ and state $0$. The number of triples is

$$D=\binom{p+2}{2}=\frac{(p+1)(p+2)}{2},$$

the number of ways to place two separators among $p+2$ positions, namely $6$ for $m=1$, $15$ for $m=2$ and $231$ for the required $m=10$.

5. The transition for one interior square

Let $F_t(a,b,c)$ be the number of ordered $p$-tuples of legal partial legs over the first $t$ interior squares leaving $a$ legs in state $2$, $b$ in state $1$ and $c$ in state $0$. At the next square all $a$ legs in state $2$ must land; choose freely which $i$ of the $b$ legs in state $1$ and which $j$ of the $c$ legs in state $0$ skip, and every other leg lands. The new triple is $(i,j,k)$ with

$$k=a+(b-i)+(c-j)=p-i-j,$$

so the new count of legs in state $0$ is forced, and

$$F_{t+1}(i,j,k)=\sum\binom{b}{i}\binom{c}{j}F_t(a,b,c),$$

the sum running over the old triples with $a+b+c=p$, $b\ge i$ and $c\ge j$. Before the first interior square every leg sits on square $1$ in state $0$, so $F_0(0,0,p)=1$. Since every final state has exactly one completion, the sum of $F_{n-2}$ over all triples counts the plans with no condition on unvisited squares.

6. The scaled basis the programs use

Fix which labelled legs sit in which group: legs $1,\dots,a$ in state $2$, the next $b$ in state $1$, the rest in state $0$. Relabelling turns the tuples of one such assignment into those of any other, so all $p!/(a!\,b!\,c!)$ assignments carry equally many, namely

$$H_t(a,b,c)=\frac{a!\,b!\,c!}{p!}\,F_t(a,b,c),$$

an integer, with $H_0(0,0,p)=1$. That is what the programs store. Write the old triple as $(u,\,v+i,\,j+w)$: $u$ legs in state $2$, $v$ state-$1$ legs that land, $w$ state-$0$ legs that land, so $u+v+w=k$. After the factorials cancel, the two binomials merge into one multinomial coefficient

$$\frac{k!}{u!\,v!\,w!}=\binom{k}{u}\binom{k-u}{v},$$

which splits the $k$ legs that land here into the three groups they came from, ^[4] and the recurrence becomes

$$H_{t+1}(i,j,k)=\sum\frac{k!}{u!v!w!}H_t(u,v{+}i,j{+}w)$$

with the sum over all $u,v,w\ge0$ adding up to $k$.

Two things are gained. The entries are now the multinomials the programs store, and since $H_0=1$ and every entry is an integer the computation needs only integer addition and multiplication: it may run modulo $10^9=2^9\cdot5^9$ although that modulus is composite, because congruence classes form a ring for any modulus and the factorial quotient only derives the recurrence, never divides inside it. ^[5] And undoing the scaling is free, since

$$F_t(a,b,c)=\frac{p!}{a!\,b!\,c!}\,H_t(a,b,c),$$

so the total over all states is $H_t$ contracted with the multinomial row $p!/(a!\,b!\,c!)$ — which is already stored, since the recurrence with $i=j=0$ and $k=p$ produces exactly that as the matrix row indexed by $(0,0,p)$. No coincidence: one more interior square on which every leg lands cannot change the count and collapses every state to $(0,0,p)$, so applying that row is the same as summing over all states. Write $T$ for the matrix and $L$ for that closure row. A row of $T$ is a new state and a column an old one, so a column of counts advances one square when $T$ multiplies it from the left — the opposite of the usual adjacency-matrix convention. ^[8]

7. Marking the squares nobody visits

Let $g_q(n)$ be the number of plans in which exactly $q$ interior squares are skipped by every leg, so the answer is $g_0(n)+g_1(n)$. Give each such square a marker $y$ and every other square the weight $1$; a plan then carries the weight $y^q$, and the weighted total is

$$P_n(y)=\sum_{q\ge0}g_q(n)\,y^q.$$

For $m=1$ and $n=4$ the $16$ ordered plans split into $9$, $6$ and $1$ by the number of squares they miss, so $P_4(y)=9+6y+y^2$ and we want $9+6$.

All legs skip a square only if none is in state $2$, and then every state-$1$ leg moves to state $2$ and every state-$0$ leg to state $1$. The all-skip event is thus deterministic: one way only, sending $(0,b,c)$ to $(b,c,0)$. Let $Z$ be the matrix of that single transition; its nonzero entries are $1$ in the scaled basis too, since the scaling factor cancels between $(0,b,c)$ and $(b,c,0)$. As $T$ already counts the all-skip option once with weight $1$, the weighted matrix of one interior square is

$$N(y)=(T-Z)+yZ.$$

8. Cutting the polynomial off at $y^2$

Only the coefficients of $y^0$ and $y^1$ are wanted, so we may work modulo $y^2$. There every element is a pair $A+By$ with $y^2=0$ — a dual number — and

$$(A+By)(C+Dy)=AC+(AD+BC)y.$$

A product never lowers a power of $y$, so a term of degree two or more can never reach the constant or the linear coefficient: dropping those terms after every product is exact. For a matrix polynomial $Q(y)$ the coefficient of $y$ is $Q'(0)$, the entrywise derivative in the scalar $y$ at $0$, and that is why the programs call the second component a derivative. ^[6]

Here $N(y)=M+yZ$ with $M=T-Z$, so the pair is $(M,Z)$ and the pair of $N(y)^{\,n-2}$ is $(M^{\,n-2},R)$, with $R$ assembled by the product rule above. The order of its two products matters: at the exponent $2$ the linear coefficient is $MZ+ZM$, and these two matrices differ. With $e_0$ the unit vector at $(0,0,p)$,

$$\begin{aligned}g_0(n)&=L\,M^{\,n-2}e_0,\\ g_1(n)&=L\,R\,e_0.\end{aligned}$$

As $Z$ is zero in the row $(0,0,p)$, the matrices $T$ and $M$ agree there, so one stored matrix serves as both the transition and the closure row $L$.

9. Reaching the far end by squaring

Put $e=n-2$. Matrix multiplication is associative, so the powers with exponents $1,2,4,8,\dots$ can be built by squaring and the ones picked out by the binary digits of $e$ multiplied together: $\lfloor\log_2 e\rfloor$ squarings plus at most as many further products in place of $e$ single steps. ^[7] For $n=10^{12}$ the exponent has 40 binary digits of which 23 are one, so 39 squarings and 23 further products replace a trillion interior squares — each tripled, because every step carries a pair. Nothing in Sections 1 to 9 is an estimate, so the count is exact before the modulus is applied.

10. Worked example: one round trip

Take $m=1$, so $p=2$ and $D=6$. Number the six states $1$ to $6$ in the order $(0,0,2)$, $(0,1,1)$, $(0,2,0)$, $(1,0,1)$, $(1,1,0)$, $(2,0,0)$. The closure row of Section 6 is the list of multinomials $2!/(a!\,b!\,c!)$, that is $L=(1,2,1,2,2,1)$. The matrix $Z$ copies $V_1$ to position $3$, $V_2$ to $5$ and $V_3$ to $6$, while $M=T-Z$ has three nonzero rows only: row $1$ is $L$ itself, by the identity of Section 6, row $2$ returns $V_1+V_2+V_4$ and row $4$ returns $V_2+V_3+V_5$. The empty rows are the states with no leg in state $0$, reachable only by an all-skip square.

Working on pairs as in Section 8 means iterating two vectors together,

$$V_{t+1}=MV_t,\quad R_{t+1}=ZV_t+MR_t,$$

from $R_0=0$ and $V_0$ the unit vector at state $1$. Three steps give

$t=1$: $V=(1,1,0,0,0,0)$ and $R=(0,0,1,0,0,0)$, so $g_0=3$ and $g_1=1$: $F(1,3)=4$;
$t=2$: $V=(3,2,0,1,0,0)$ and $R=(1,0,1,1,1,0)$, so $g_0=9$ and $g_1=6$: $F(1,4)=15$;
$t=3$: $V=(9,6,0,2,0,0)$ and $R=(6,2,3,2,2,0)$, so $g_0=25$ and $g_1=21$: $F(1,5)=46$.

The last line is worth doing by hand: contracting with $L$ gives $9+12+4=25$ for $V$ and $6+4+3+4+4$${}=21$ for $R$.

One run thus yields all three checkpoints for $m=1$. The middle line is the distribution used in Section 7: for $n=4$ the legs are $1\to4$, $1\to2\to4$, $1\to3\to4$ and $1\to2\to3\to4$, and of the $16$ ordered pairs only the one in which both legs jump straight from $1$ to $4$ misses both interior squares.

How the Code Works

The C++ program builds the matrices of Sections 6 and 7 once, then runs Sections 8 and 9. Matrix is a flat std::vector of 32-bit words with an operator() for row and column, plus identity, add, subtract, transpose and multiply.

build_basis enumerates the triples in the order $a=0,\dots,p$ and $b=0,\dots,p-a$, fills the lookup table idx, and records in idx_c_pow the position of $(0,0,p)$ — at once the start vector $e_0$ and the index of the closure row. build_combinations fills Pascal's triangle up to row $p$. build_allowed walks each new triple $(i,j,k)$, loops $u$ to $k$ and $v$ to $k-u$, looks up the old triple $(u,\,v+i)$ and adds comb[k][u] * comb[k-u][v]; this is $T$. Its largest entry for $p=20$ is 133024320, below the modulus, so the stored numbers are the true multinomial coefficients. build_forbidden writes $Z$, and the constructor of Solver forms combined${}=M=T-Z$.

solve sets internal${}=n-2$ and calls pow_with_derivative(combined, forbidden, internal), which walks the bits of the exponent keeping two matrix pairs: cur with cur_deriv for the current repeated square, result with result_deriv for the accumulated product. Each update of a second component is the two-term product rule of Section 8, in the correct order. It then takes the column idx_c_pow of the two returned matrices, which are $M^{\,n-2}e_0$ and $Re_0$, contracts each with the row idx_c_pow of allowed — that is with $L$ — and returns their sum modulo $10^9$, which main prints. multiply transposes its right operand first, so both inner loops walk along rows of memory, and accumulates a whole dot product in an unsigned __int128, a compiler extension, ^[9] since $231$ products of numbers below $10^9$ reach about $2.3\cdot10^{20}$.

The program carries its own brute force too: build_path_masks lists every leg as a bitmask of visited squares, and brute_count convolves $p$ copies of that list over all $2^n$ subsets, keeping those that miss at most one interior square. run_validation checks the five statement values, then compares the matrix method with brute_count for $m=1,2$ and $n=3,\dots,7$; -m and -n select other parameters, --no-validate skips the checks. One warning about -m: build_allowed multiplies two 64-bit binomial coefficients before reducing, which holds only while every multinomial stays below $2^{64}$, that is up to $m=22$. At $m=23$ the largest coefficient wraps silently; the required $m=10$ is far inside the safe range.

Java runs the same algorithm on flat long[] arrays. Its matMul, used for the value component, skips zero entries of the left operand — which pays at the start, since $T$ is mostly zero — and reduces after every multiply-and-add, while the second component is built by explicit loops in powDeriv that reduce each product before adding it to the sums s1 and s2. Nothing then exceeds $2D\cdot10^9$, so a signed 64-bit long suffices. Here $m$ and $n$ are hard-coded and no checks run. The Python file is a bridge, not a port: it looks for solutionsCpp/Euler416.cpp beside its own parent directory, compiles it with clang++ or g++ into a cached binary when that binary is missing or older than the source, and runs it with the checks still on. parse_output tries two patterns for a labelled result, but this program prints a bare number, so the last nonempty line comes back. The Python result is the C++ result by construction.

Complexity and Verification

The state space has $D=(p+1)(p+2)/2$ elements, of order $m^2$ and equal to $231$ for $m=10$, so one dense matrix product performs $D^3=12326391$ multiply-and-add steps and, in C++, $D^2=53361$ reductions. The exponent $n-2$ is 999999999998, with 40 binary digits of which 23 are one, and pow_with_derivative spends three products per squaring and three per set bit, two of them for the second component: $3\cdot(23+39)=186$ matrix products, about $2.3\cdot10^9$ multiply-and-add steps. The time is $O(m^6\log n)$ and the memory $O(m^4)$: one matrix holds $231^2$ entries, 213444 bytes in C++ and twice that in Java. Measured on one current desktop, the C++ program finishes in well under a second and the Java program in about two.

Dense products are right even though the input is sparse. For $p=20$ the matrix $T$ has 10626 nonzero entries out of 53361, about 19.9 per cent, and $M$ has 21 empty rows, exactly the triples with $c=0$; those rows stay empty in every power, but the rest fills in at once: $M^2$ has 26565 nonzero entries, and $M^3$ and $M^4$ have 48510 each, about 90.9 per cent. The brute force is a checking aid only: its $2^n$ unreduced 64-bit counters would overflow at $n=6$ already if called with $m=10$.

Beyond the checks packaged with the C++ program, the following were verified independently for this article.

An enumeration listing each leg as its set of visited squares and grouping all $p$-tuples by the union reproduces $F(1,3)=4$, $F(1,4)=15$, $F(1,5)=46$ and $F(2,3)=16$, with the splits by unvisited squares 3 + 1, 9 + 6 + 1, 25 + 21 + 3 and 15 + 1.
A second transfer-matrix program, written from the derivation above and not from the published source, agrees with that enumeration for all $m=1,2,3$ and $n=2,3,\dots,9$ and reproduces $F(2,100)\bmod 10^9$${}=429619151$.
The leg counts match the tribonacci sequence up to $n=11$, the row $(0,0,p)$ of $T$ equals $p!/(a!\,b!\,c!)$ for $p=20$, and C++ and Java print the same value at $m=10$, $n=10^{12}$.

Footnotes and References

The sources below supply the background facts used above. The two one-way trips, the three states, the scaled basis with its closure row and the marker $y$ are derived in this article.

Project Euler 416 — A Frog's Trip. The official statement sends a frog from the leftmost of n squares to the rightmost and back, m times, in jumps of one, two or three squares, and defines F(m, n) as the number of such travels leaving at most one square unvisited, with the five checkpoints quoted in the summary.
Wikipedia — Composition (combinatorics). A composition writes n as an ordered sum of positive integers; A-restricted compositions take their parts from a fixed set A. A leg in Section 2 is a composition of n − 1 with parts in {1, 2, 3}.
OEIS A000073 — Tribonacci numbers. Defined by a(n) = a(n−1) + a(n−2) + a(n−3) with a(0) = a(1) = 0 and a(2) = 1. Two comments state that a(n) counts the compositions of n − 2 with no part above 3 and the binary sequences of length n − 3 without three consecutive zeros: the legs of Section 2 in both readings.
Wikipedia — Multinomial theorem. It defines the multinomial coefficient n!/(k₁! ⋯ k_m!), counting the ways of depositing n distinct objects into m distinct bins with kᵢ in bin i, and gives the product form C(n, k₁)·C(n − k₁, k₂)⋯.
Wikipedia — Modular arithmetic. The congruence classes modulo m form a commutative ring, a field exactly when m is prime, and a class is invertible exactly when it is coprime to m — Section 6 uses the ring operations only.
Wikipedia — Dual number. Dual numbers are the expressions a + bε with ε² = 0, multiplying by (a + bε)(c + dε) = ac + (ad + bc)ε. The article also derives P(a + bε) = P(a) + bP′(a)ε, the identity behind Section 8.
Wikipedia — Exponentiation by squaring. The method needs ⌊log₂ n⌋ squarings and at most that many further multiplications, one less than the number of ones in the binary expansion of n. It works in non-commutative semigroups, hence for matrices.
Wikipedia — Adjacency matrix. Its section on matrix powers states that entry (i, j) of the n-th power is the number of walks of length n from vertex i to vertex j — the orientation opposite to the one used here.
GCC manual — 128-bit integers. The types __int128 and unsigned __int128 are a GCC extension, supported where an integer mode is wide enough for 128 bits.

Problem 416 source code

C++

#include <algorithm>
#include <cstdint>
#include <functional>
#include <cstdlib>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

namespace {

using u32 = uint32_t;
using u64 = uint64_t;

constexpr u32 kMod = 1000000000u;
constexpr int kDefaultM = 10;
constexpr long long kDefaultN = 1000000000000LL;

struct Matrix {
    int n = 0;
    std::vector<u32> a;

    Matrix() = default;
    explicit Matrix(int n_) : n(n_), a(static_cast<size_t>(n_) * n_, 0) {}

    u32 &operator()(int r, int c) { return a[static_cast<size_t>(r) * n + c]; }
    const u32 &operator()(int r, int c) const { return a[static_cast<size_t>(r) * n + c]; }
};

Matrix identity(int n) {
    Matrix m(n);
    for (int i = 0; i < n; ++i) m(i, i) = 1;
    return m;
}

Matrix transpose(const Matrix &m) {
    Matrix t(m.n);
    for (int i = 0; i < m.n; ++i) {
        for (int j = 0; j < m.n; ++j) {
            t(j, i) = m(i, j);
        }
    }
    return t;
}

Matrix add(const Matrix &a, const Matrix &b) {
    Matrix c(a.n);
    const size_t size = a.a.size();
    for (size_t i = 0; i < size; ++i) {
        u32 val = a.a[i] + b.a[i];
        if (val >= kMod) val -= kMod;
        c.a[i] = val;
    }
    return c;
}

Matrix subtract(const Matrix &a, const Matrix &b) {
    Matrix c(a.n);
    const size_t size = a.a.size();
    for (size_t i = 0; i < size; ++i) {
        u32 val = (a.a[i] >= b.a[i]) ? (a.a[i] - b.a[i]) : (a.a[i] + kMod - b.a[i]);
        c.a[i] = val;
    }
    return c;
}

Matrix multiply(const Matrix &a, const Matrix &b, int threads) {
    const int n = a.n;
    Matrix c(n);
    Matrix bt = transpose(b);

    auto worker = [&](int row_start, int row_end) {
        for (int i = row_start; i < row_end; ++i) {
            const u32 *row_a = &a.a[static_cast<size_t>(i) * n];
            for (int j = 0; j < n; ++j) {
                const u32 *row_b = &bt.a[static_cast<size_t>(j) * n];
                unsigned __int128 sum = 0;
                for (int k = 0; k < n; ++k) {
                    sum += static_cast<unsigned __int128>(row_a[k]) * row_b[k];
                }
                c(i, j) = static_cast<u32>(sum % kMod);
            }
        }
    };

    if (threads <= 1 || n < 64) {
        worker(0, n);
        return c;
    }

    int use_threads = std::min(threads, n);
    int chunk = (n + use_threads - 1) / use_threads;
    std::vector<std::thread> pool;
    pool.reserve(use_threads);
    for (int t = 0; t < use_threads; ++t) {
        int start = t * chunk;
        int end = std::min(n, start + chunk);
        if (start >= end) break;
        pool.emplace_back(worker, start, end);
    }
    for (auto &th : pool) th.join();
    return c;
}

struct Triple {
    int a = 0;
    int b = 0;
    int c = 0;
};

struct Basis {
    int p = 0;
    int dim = 0;
    int idx_c_pow = 0;
    std::vector<Triple> triples;
    std::vector<std::vector<int>> idx;
};

Basis build_basis(int p) {
    Basis basis;
    basis.p = p;
    basis.idx.assign(p + 1, std::vector<int>(p + 1, -1));
    for (int a = 0; a <= p; ++a) {
        for (int b = 0; b <= p - a; ++b) {
            int c = p - a - b;
            basis.idx[a][b] = static_cast<int>(basis.triples.size());
            basis.triples.push_back({a, b, c});
        }
    }
    basis.dim = static_cast<int>(basis.triples.size());
    basis.idx_c_pow = basis.idx[0][0];
    return basis;
}

std::vector<std::vector<u64>> build_combinations(int p) {
    std::vector<std::vector<u64>> comb(p + 1, std::vector<u64>(p + 1, 0));
    comb[0][0] = 1;
    for (int n = 1; n <= p; ++n) {
        comb[n][0] = comb[n][n] = 1;
        for (int k = 1; k < n; ++k) {
            comb[n][k] = comb[n - 1][k - 1] + comb[n - 1][k];
        }
    }
    return comb;
}

Matrix build_allowed(const Basis &basis, const std::vector<std::vector<u64>> &comb) {
    Matrix allowed(basis.dim);
    for (int idx_new = 0; idx_new < basis.dim; ++idx_new) {
        const Triple &t = basis.triples[idx_new];
        int i = t.a;
        int j = t.b;
        int k = t.c;
        for (int u = 0; u <= k; ++u) {
            for (int v = 0; v <= k - u; ++v) {
                int w = k - u - v;
                int a_old = u;
                int b_old = v + i;
                int idx_old = basis.idx[a_old][b_old];
                u64 coeff = comb[k][u] * comb[k - u][v];
                u32 add = static_cast<u32>(coeff % kMod);
                u32 &cell = allowed(idx_new, idx_old);
                cell = static_cast<u32>((cell + add) % kMod);
            }
        }
    }
    return allowed;
}

Matrix build_forbidden(const Basis &basis) {
    Matrix forbidden(basis.dim);
    for (int idx_new = 0; idx_new < basis.dim; ++idx_new) {
        const Triple &t = basis.triples[idx_new];
        if (t.c != 0) continue;
        int idx_old = basis.idx[0][t.a];
        forbidden(idx_new, idx_old) = 1;
    }
    return forbidden;
}

struct PowResult {
    Matrix value;
    Matrix deriv;
};

PowResult pow_with_derivative(const Matrix &base, const Matrix &base_deriv, u64 exp, int threads) {
    const int n = base.n;
    Matrix result = identity(n);
    Matrix result_deriv(n);
    Matrix cur = base;
    Matrix cur_deriv = base_deriv;

    u64 e = exp;
    while (e > 0) {
        if (e & 1ULL) {
            Matrix tmp1 = multiply(result_deriv, cur, threads);
            Matrix tmp2 = multiply(result, cur_deriv, threads);
            Matrix next_deriv = add(tmp1, tmp2);
            Matrix next = multiply(result, cur, threads);
            result = std::move(next);
            result_deriv = std::move(next_deriv);
        }
        e >>= 1ULL;
        if (e == 0) break;
        Matrix tmp1 = multiply(cur_deriv, cur, threads);
        Matrix tmp2 = multiply(cur, cur_deriv, threads);
        Matrix next_deriv = add(tmp1, tmp2);
        Matrix next = multiply(cur, cur, threads);
        cur = std::move(next);
        cur_deriv = std::move(next_deriv);
    }
    return {std::move(result), std::move(result_deriv)};
}

struct Solver {
    int m = 0;
    int p = 0;
    Basis basis;
    Matrix allowed;
    Matrix forbidden;
    Matrix combined;

    explicit Solver(int m_) : m(m_), p(2 * m_), basis(build_basis(2 * m_)) {
        auto comb = build_combinations(p);
        allowed = build_allowed(basis, comb);
        forbidden = build_forbidden(basis);
        combined = subtract(allowed, forbidden);  // y = -1
    }

    u32 solve(long long n, int threads) const {
        if (n <= 1) return 1;
        u64 internal = (n >= 2) ? static_cast<u64>(n - 2) : 0;
        PowResult pow = (internal == 0)
                            ? PowResult{identity(basis.dim), Matrix(basis.dim)}
                            : pow_with_derivative(combined, forbidden, internal, threads);

        const int idx0 = basis.idx_c_pow;
        std::vector<u32> col(pow.value.n);
        std::vector<u32> col_deriv(pow.value.n);
        for (int i = 0; i < pow.value.n; ++i) {
            col[i] = pow.value(i, idx0);
            col_deriv[i] = pow.deriv(i, idx0);
        }

        unsigned __int128 sum_val = 0;
        unsigned __int128 sum_deriv = 0;
        const u32 *row = &allowed.a[static_cast<size_t>(idx0) * allowed.n];
        for (int j = 0; j < allowed.n; ++j) {
            sum_val += static_cast<unsigned __int128>(row[j]) * col[j];
            sum_deriv += static_cast<unsigned __int128>(row[j]) * col_deriv[j];
        }
        u32 f_val = static_cast<u32>(sum_val % kMod);
        u32 f_deriv = static_cast<u32>(sum_deriv % kMod);
        u32 ans = f_val + f_deriv;
        if (ans >= kMod) ans -= kMod;
        return ans;
    }
};

std::vector<int> build_path_masks(int n) {
    std::vector<int> masks;
    if (n <= 0) return masks;

    std::function<void(int, int)> dfs = [&](int pos, int mask) {
        if (pos == n) {
            masks.push_back(mask);
            return;
        }
        for (int step = 1; step <= 3; ++step) {
            int next = pos + step;
            if (next > n) continue;
            int next_mask = mask | (1 << (next - 1));
            dfs(next, next_mask);
        }
    };

    dfs(1, 1);
    return masks;
}

u64 brute_count(int m, int n) {
    if (n <= 1) return 1;
    std::vector<int> masks = build_path_masks(n);
    int legs = 2 * m;
    int max_mask = 1 << n;
    std::vector<u64> dp(max_mask, 0);
    dp[0] = 1;
    for (int step = 0; step < legs; ++step) {
        std::vector<u64> next(max_mask, 0);
        for (int mask = 0; mask < max_mask; ++mask) {
            u64 val = dp[mask];
            if (val == 0) continue;
            for (int pmask : masks) {
                next[mask | pmask] += val;
            }
        }
        dp.swap(next);
    }

    int internal = std::max(0, n - 2);
    int internal_mask = ((1 << n) - 1) ^ 1 ^ (1 << (n - 1));
    u64 total = 0;
    for (int mask = 0; mask < max_mask; ++mask) {
        u64 val = dp[mask];
        if (val == 0) continue;
        int visited_internal = __builtin_popcount(mask & internal_mask);
        int unvisited = internal - visited_internal;
        if (unvisited <= 1) total += val;
    }
    return total % kMod;
}

void run_validation() {
    struct Check {
        int m;
        long long n;
        u32 expected;
    };
    const Check checks[] = {
        {1, 3, 4},
        {1, 4, 15},
        {1, 5, 46},
        {2, 3, 16},
        {2, 100, 429619151},
    };
    for (const auto &chk : checks) {
        Solver solver(chk.m);
        u32 got = solver.solve(chk.n, 1);
        if (got != chk.expected) {
            std::cerr << "Validation failed for m=" << chk.m << ", n=" << chk.n
                      << ": got " << got << ", expected " << chk.expected << '\n';
            std::exit(1);
        }
    }

    for (int m = 1; m <= 2; ++m) {
        Solver solver(m);
        for (int n = 3; n <= 7; ++n) {
            u32 fast = solver.solve(n, 1);
            u32 brute = static_cast<u32>(brute_count(m, n));
            if (fast != brute) {
                std::cerr << "Brute validation failed for m=" << m << ", n=" << n
                          << ": got " << fast << ", expected " << brute << '\n';
                std::exit(1);
            }
        }
    }
    std::cerr << "Validation checkpoints passed.\n";
}

void usage(const char *argv0) {
    std::cerr << "Usage: " << argv0 << " [-m M] [-n N] [-t THREADS] [--no-validate]\n";
}

}  // namespace

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

    int m = kDefaultM;
    long long n = kDefaultN;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;
    bool validate = true;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "-m" && i + 1 < argc) {
            m = std::stoi(argv[++i]);
        } else if (arg == "-n" && i + 1 < argc) {
            n = std::stoll(argv[++i]);
        } else if (arg == "-t" && i + 1 < argc) {
            threads = std::max(1, std::stoi(argv[++i]));
        } else if (arg == "--no-validate") {
            validate = false;
        } else if (arg == "--validate") {
            validate = true;
        } else {
            usage(argv[0]);
            return 1;
        }
    }

    if (validate) run_validation();

    Solver solver(m);
    u32 answer = solver.solve(n, threads);
    std::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

public class Euler416 {
    static final long MOD = 1000000000;

    public static String solve() {
        int m = 10;
        long n = 1000000000000L;
        int p = 2 * m;
        int dim = 0;
        int[][] triples;
        int[][] idx = new int[p + 1][p + 1];
        java.util.List<int[]> tl = new java.util.ArrayList<>();
        for (int a = 0; a <= p; a++)
            for (int b = 0; b <= p - a; b++) {
                idx[a][b] = tl.size();
                tl.add(new int[] { a, b, p - a - b });
            }
        dim = tl.size();
        triples = tl.toArray(new int[0][]);
        int idxCP = idx[0][0];
        long[][] C = new long[p + 1][p + 1];
        C[0][0] = 1;
        for (int nn = 1; nn <= p; nn++) {
            C[nn][0] = C[nn][nn] = 1;
            for (int k = 1; k < nn; k++)
                C[nn][k] = C[nn - 1][k - 1] + C[nn - 1][k];
        }
        long[] allowed = new long[dim * dim], forbidden = new long[dim * dim];
        for (int inew = 0; inew < dim; inew++) {
            int at = triples[inew][0], bt = triples[inew][1], ct = triples[inew][2];
            for (int u = 0; u <= ct; u++)
                for (int v = 0; v <= ct - u; v++) {
                    int ao = u, bo = v + at, iold = idx[ao][bo];
                    long coeff = C[ct][u] * C[ct - u][v] % MOD;
                    allowed[inew * dim + iold] = (allowed[inew * dim + iold] + coeff) % MOD;
                }
        }
        for (int inew = 0; inew < dim; inew++) {
            int ct = triples[inew][2];
            if (ct != 0)
                continue;
            int iold = idx[0][triples[inew][0]];
            forbidden[inew * dim + iold] = 1;
        }
        long[] combined = new long[dim * dim];
        for (int i = 0; i < dim * dim; i++)
            combined[i] = (allowed[i] - forbidden[i] + MOD) % MOD;
        long internal = n - 2;
        long[] pv, pd;
        if (internal == 0) {
            pv = matId(dim);
            pd = new long[dim * dim];
        } else {
            long[][] res = powDeriv(combined, forbidden, internal, dim);
            pv = res[0];
            pd = res[1];
        }
        long sv = 0, sd = 0;
        for (int j = 0; j < dim; j++) {
            sv += allowed[idxCP * dim + j] * pv[j * dim + idxCP] % MOD;
            sd += allowed[idxCP * dim + j] * pd[j * dim + idxCP] % MOD;
        }
        return String.valueOf((sv % MOD + sd % MOD) % MOD);
    }

    static long[] matMul(long[] x, long[] y, int d) {
        long[] c = new long[d * d];
        for (int i = 0; i < d; i++)
            for (int k = 0; k < d; k++) {
                long xik = x[i * d + k];
                if (xik == 0)
                    continue;
                for (int j = 0; j < d; j++)
                    c[i * d + j] = (c[i * d + j] + xik * y[k * d + j]) % MOD;
            }
        return c;
    }

    static long[] matId(int d) {
        long[] m = new long[d * d];
        for (int i = 0; i < d; i++)
            m[i * d + i] = 1;
        return m;
    }

    static long[][] powDeriv(long[] base, long[] deriv, long exp, int d) {
        long[] rr = matId(d), rd = new long[d * d];
        while (exp > 0) {
            if ((exp & 1) != 0) {
                long[] nrd = new long[d * d];
                for (int i = 0; i < d; i++)
                    for (int j = 0; j < d; j++) {
                        long s1 = 0, s2 = 0;
                        for (int k = 0; k < d; k++) {
                            s1 += rd[i * d + k] * base[k * d + j] % MOD;
                            s2 += rr[i * d + k] * deriv[k * d + j] % MOD;
                        }
                        nrd[i * d + j] = (s1 + s2) % MOD;
                    }
                rd = nrd;
                rr = matMul(rr, base, d);
            }
            exp >>= 1;
            if (exp == 0)
                break;
            long[] nd = new long[d * d];
            for (int i = 0; i < d; i++)
                for (int j = 0; j < d; j++) {
                    long s1 = 0, s2 = 0;
                    for (int k = 0; k < d; k++) {
                        s1 += deriv[i * d + k] * base[k * d + j] % MOD;
                        s2 += base[i * d + k] * deriv[k * d + j] % MOD;
                    }
                    nd[i * d + j] = (s1 + s2) % MOD;
                }
            deriv = nd;
            base = matMul(base, base, d);
        }
        return new long[][] { rr, rd };
    }

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