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

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

Call these states $S_0,S_1,S_2$, where $S_r$ means “the current trailing run of skipped internal nodes has length $r$”. A visited node resets the state to $S_0$. A skipped node sends $S_0 \to S_1$ and $S_1 \to S_2$. The transition $S_2 \to S_3$ is impossible, because that would require a jump of length at least $4$.

This is the key compression: every leg has many possible detailed histories, but for the next step only the current zero-run length matters.

Step 3: Compress All $p=2m$ Legs into a Triple

Now process all $p$ legs simultaneously. After a prefix of internal nodes, let

$$F_t(a,b,c)$$

be the number of ordered $p$-tuples of partial legs such that exactly $a$ legs are in state $S_2$, exactly $b$ legs are in state $S_1$, and exactly $c$ legs are in state $S_0$, with

$$a+b+c=p.$$

The initial state is obvious: before any internal node is processed, no leg has skipped anything yet, so

$$F_0(0,0,p)=1,$$

and all other triples are zero.

The number of possible triples is

$$D=\#\{(a,b,c)\in \mathbb{Z}_{\ge 0}^3:a+b+c=p\}=\binom{p+2}{2}.$$

Since $p=2m$, this is quadratic in $m$. For the actual target $m=10$, we get $p=20$ and $D=231$.

Step 4: Transition for One More Internal Node

Suppose the old triple is $(A,B,C)$. At the next internal node:

If a leg is in $S_2$, it must visit the node, otherwise the path would contain three consecutive skips.

If a leg is in $S_1$, it may either visit the node and reset to $S_0$, or skip it and move to $S_2$.

If a leg is in $S_0$, it may either visit the node and stay effectively in $S_0$, or skip it and move to $S_1$.

Therefore, if the new triple is $(i,j,k)$, then $i$ must be the number of old $S_1$-legs that skip, and $j$ must be the number of old $S_0$-legs that skip. The rest visit. In the direct counting basis, the recurrence is

$$F_{t+1}(i,j,k)=\sum \binom{B}{i}\binom{C}{j}F_t(A,B,C),$$

where the sum runs over old triples satisfying

$$k=A+(B-i)+(C-j),\qquad A+B+C=p.$$

This already solves the combinatorics, but the implementations use a slightly different basis because it turns the coefficients into the multinomial form used by the matrix construction.

Step 5: Factorial-Scaled Basis and Matrix Entries

Define scaled coordinates

$$H_t(a,b,c)=a!\,b!\,c!\,F_t(a,b,c).$$

Write the old triple as $(u,v+i,j+w)$, so that $u+v+w=k$. After substituting the direct recurrence and cancelling factorials, we obtain

$$H_{t+1}(i,j,k)=\sum_{u+v+w=k}\binom{k}{u}\binom{k-u}{v}\,H_t(u,v+i,j+w).$$

This is exactly the transition rule used by the implementations. The coefficient

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

is a multinomial count: among the $k$ legs that do visit the current node, choose how many came from the old $S_2$, old $S_1$, and old $S_0$ groups.

So the dense matrix is indexed only by triples $(a,b,c)$, not by individual legs. That symmetry reduction is what makes the problem tractable.

Step 6: Mark Internal Nodes Missed by Every Leg

Let $g_q(n)$ denote the number of ordered $p$-tuples of full legs for which exactly $q$ internal nodes are missed by every leg. The target of the problem is

$$g_0(n)+g_1(n).$$

Introduce the generating polynomial

$$G_n(y)=\sum_{q\ge 0} g_q(n)(1+y)^q.$$

At one internal node there are two possibilities: either at least one leg visits it, contributing a factor $1$, or every leg skips it, contributing a factor $1+y$.

Let $T$ be the total transition matrix from Step 5, and let $Z$ be the submatrix corresponding to the special case “every leg skips the current node”. That special case is deterministic:

$$ (0,a,p-a)\longrightarrow (a,p-a,0). $$

Hence the weighted one-node transfer is

$$M(y)=T+yZ,$$

because $T$ already includes the all-skip case once, and the extra $yZ$ changes its weight from $1$ to $1+y$.

Step 7: Initial Vector, Final Closure, and the Derivative Filter

Let $e_0$ be the basis vector of the initial triple $(0,0,p)$. After the $n-2$ internal nodes are processed, we get

$$V_{n-2}(y)=M(y)^{n-2}e_0.$$

There is still a final jump to node $n$. A leg finishing in state $S_0$, $S_1$, or $S_2$ must use a final jump of length $1$, $2$, or $3$, respectively, so every terminal automaton state has exactly one legal completion. In the factorial-scaled basis, the corresponding closure row is obtained from the row indexed by $(0,0,p)$ in the total transition matrix. Denote this row by $L$.

Therefore

$$G_n(y)=L\,M(y)^{n-2}e_0.$$

Now evaluate at $y=-1$. Since $(1-1)^q$ vanishes for every $q\ge 1$, we get

$$G_n(-1)=g_0(n).$$

Also,

$$\left.\frac{d}{dy}(1+y)^q\right|_{y=-1}=\begin{cases}1,&q=1,\\0,&q\neq 1,\end{cases}$$

$$G_n'(-1)=g_1(n).$$

Hence the required count is

$$\boxed{g_0(n)+g_1(n)=G_n(-1)+G_n'(-1)\pmod{10^9}.}$$

Worked Example: $m=1,\ n=4$

Here $p=2$ and the internal nodes are $2$ and $3$. A single leg can be one of four paths:

$1\to4$, $1\to2\to4$, $1\to3\to4$, or $1\to2\to3\to4$.

So there are $4^2=16$ ordered pairs of legs. The only invalid pair is the pair in which both legs jump directly $1\to4$, because then both internal nodes are missed and we have $q=2$.

Therefore

$$g_0(4)=9,\qquad g_1(4)=6,\qquad g_2(4)=1,$$

and the desired value is

$$g_0(4)+g_1(4)=9+6=15.$$

This matches the small checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematics. They enumerate all triples $(a,b,c)$ with $a+b+c=2m$, precompute the needed binomial coefficients, build the dense total-transition matrix and the all-skip submatrix, and work entirely modulo $10^9$.

The hard part is the huge exponent $n-2$, so the implementation does not iterate over vertices one by one. Instead it evaluates $M(-1)^{n-2}$ and its derivative at $y=-1$ by binary matrix exponentiation. If $U(y)$ and $V(y)$ are matrix-valued factors, the derivative is updated by the product rule

$$\frac{d}{dy}(UV)=U'V+UV'.$$

After exponentiation, the implementation applies the terminal closure row, extracts $G_n(-1)$ and $G_n'(-1)$, adds them, and reduces the result modulo $10^9$. The Python version delegates to the same compiled algorithm, so all three language implementations compute the identical transfer-matrix formula.

Complexity Analysis

The state-space dimension is

$$D=\binom{2m+2}{2}=O(m^2).$$

Dense matrix multiplication dominates the runtime. Binary exponentiation needs $O(\log n)$ matrix products, so the overall complexity is

$$O(D^3\log n)$$

time and

$$O(D^2)$$

memory. The preliminary binomial table is only $O(m^2)$ and is negligible compared with the matrix work.

Footnotes and References

Problem page: https://projecteuler.net/problem=416
Generating function: Wikipedia — Generating function
Matrix exponentiation / binary exponentiation: cp-algorithms — Binary Exponentiation
Multinomial theorem: Wikipedia — Multinomial theorem
Transfer-matrix method: Wikipedia — Transfer-matrix method

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());
    }
}