Problem 996: Overtakes

Project Euler Problem 996 Solution

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

Problem Summary There are $n$ players arranged in their initial rank order. Each day a match is played between adjacent ranks. If the higher-ranked player wins, the order and all overtake counts are unchanged. If the lower-ranked player wins, the two adjacent players swap, and the winner's overtake count increases by $1$. After exactly $k$ days the final order must again be the initial order, and $F(n,k)$ counts how many overtake-count vectors can occur. The important observation is that the non-overtake matches are pure idle steps. They can be inserted anywhere without changing either the order or the count vector. Therefore only the sequence of actual overtakes matters, and because every overtake is an adjacent transposition, returning to the identity permutation requires an even number of overtakes. The implementation writes $E=\lfloor k/2\rfloor$ and counts all feasible count vectors with at most $E$ paired crossings. Mathematical Approach From overtakes to pair crossings Consider two fixed players $i$ and $j$, with $i$ initially above $j$. Their relative order can change only when they are adjacent and one overtakes the other. Since the final order is again the initial order, this pair must cross an even number of times....

Detailed mathematical approach

Problem Summary

There are $n$ players arranged in their initial rank order. Each day a match is played between adjacent ranks. If the higher-ranked player wins, the order and all overtake counts are unchanged. If the lower-ranked player wins, the two adjacent players swap, and the winner's overtake count increases by $1$. After exactly $k$ days the final order must again be the initial order, and $F(n,k)$ counts how many overtake-count vectors can occur.

The important observation is that the non-overtake matches are pure idle steps. They can be inserted anywhere without changing either the order or the count vector. Therefore only the sequence of actual overtakes matters, and because every overtake is an adjacent transposition, returning to the identity permutation requires an even number of overtakes. The implementation writes $E=\lfloor k/2\rfloor$ and counts all feasible count vectors with at most $E$ paired crossings.

Mathematical Approach

From overtakes to pair crossings

Consider two fixed players $i$ and $j$, with $i$ initially above $j$. Their relative order can change only when they are adjacent and one overtakes the other. Since the final order is again the initial order, this pair must cross an even number of times. If it crosses $2m_{ij}$ times, then $j$ overtakes $i$ exactly $m_{ij}$ times while moving upward, and $i$ overtakes $j$ exactly $m_{ij}$ times while moving back downward.

Thus every closed overtake history determines a loopless multigraph on the $n$ players: put $m_{ij}$ edges between players $i$ and $j$. The overtake count of player $i$ is exactly the degree of vertex $i$. If the total number of edges is $e$, then the total number of actual overtakes is $2e$, so a $k$-day schedule can realize any vector that is realizable with $e\le E=\lfloor k/2\rfloor$.

Why positive players form contiguous runs

Adjacent swaps impose a one-dimensional restriction. If two players interact, then every player initially between them must at some point lie inside the same moving block. Therefore the non-zero entries of a feasible degree vector decompose into disjoint contiguous runs. A run of length $1$ is impossible: a single isolated player cannot overtake anyone and return while all neighbours have zero count.

So a feasible vector is built from zero gaps and positive runs of lengths $\ell\ge 2$. Each run represents one connected component of the crossing multigraph on a consecutive block of players. This is exactly the decomposition used by the dynamic program: append a zero, or append a positive run of length $\ell$ after a zero-separated prefix.

Counting one positive run

For one run of length $\ell$ using $e$ graph edges, we need positive degrees

$$a_1,a_2,\ldots,a_\ell \ge 1,\qquad a_1+\cdots+a_\ell=2e.$$

A loopless multigraph with $e$ edges cannot have one vertex degree greater than $e$, because every incident edge contributes one degree to that vertex and one degree elsewhere. Conversely, for positive integer sequences the condition

$$\max_i a_i\le e$$

is sufficient for a loopless multigraph realization on the run. Hence the number of possible degree sequences for a run is the number of positive compositions of $2e$ into $\ell$ parts, excluding those with one part greater than $e$. Since at most one part can exceed $e$, this gives

$$R(\ell,e)=\binom{2e-1}{\ell-1}-\ell\binom{e-1}{\ell-1}.$$

This is the formula implemented in run_count_table. The first binomial counts all positive compositions; the second term subtracts the cases where a specified component is too large.

Dynamic programming over prefixes

Let $T(i,e)$ be the number of feasible vectors on the first $i$ players whose crossing multigraph has exactly $e$ edges. The code also keeps $Z(i,e)$, the number of such vectors whose last entry is zero. This makes it easy to enforce that two positive runs are separated by at least one zero.

The transition is:

$$Z(i,e)=T(i-1,e),$$

because appending a zero to a length $i-1$ vector gives a length $i$ vector ending in zero. Then a final positive run of length $\ell\ge 2$ can be appended after a zero-ending prefix:

$$T(i,e)\;{+}{=}\;Z(i-\ell,e-u)\,R(\ell,u),\qquad 2\le \ell\le i,\;1\le u\le e.$$

The base case is $T(0,0)=Z(0,0)=1$. After all $n$ players have been processed, the desired count for a day limit $k$ is cumulative:

$$F(n,k)=\sum_{e=0}^{\lfloor k/2\rfloor}T(n,e).$$

Polynomial tail for the large target

The direct DP is small in $n$, but $E=4567891/2$ is too large to iterate up to. The run formula $R(\ell,e)$ is a polynomial in $e$ of degree $\ell-1$, and the prefix recurrence is made from additions and convolutions of these polynomial pieces. For fixed $n$, the cumulative function

$$A_n(E)=\sum_{e=0}^{E}T(n,e)$$

agrees, from $E\ge n-1$ onward, with a polynomial of degree $n$. The program computes $A_n(E)$ only for the $n+1$ consecutive values

$$E=n-1,\;n,\;\ldots,\;2n-1,$$

takes finite differences, and evaluates the Newton forward expansion

$$A_n(E)=\sum_{r=0}^{n}\Delta^r A_n(n-1)\binom{E-(n-1)}{r}.$$

The target modulus is not prime, so the large binomial values are not computed by modular inverses. Instead, the implementation cancels the small denominator factors against the numerator factors over the integers, then multiplies the remaining factors modulo $1234567891$.

How the Code Works

The C++, Python, and Java programs first build binomial tables for the small range needed by the DP. They compute all $R(\ell,e)$, run the prefix dynamic program up to $E=2n-1$ for the target $n=123$, and form the cumulative values. Then they take finite differences beginning at $E=n-1$ and evaluate the degree-$n$ polynomial at $E=\lfloor4567891/2\rfloor$.

The validation block is deliberately stronger than just checking the two examples. For $2\le n\le5$, it compares the DP against a brute-force search of adjacent-swap states for small $E$. It also checks $F(3,4)=8$, $F(12,34)=2457178250$, and verifies the polynomial continuation against directly computed values for several small $n$.

Complexity Analysis

The DP only needs $E_{\max}=2n-1$ for the polynomial interpolation stage. The run table costs $O(n^2E_{\max})$ arithmetic after binomial precomputation, and the prefix DP costs $O(n^2E_{\max}^2)$ in this direct implementation. With $n=123$ and $E_{\max}=245$, this is practical.

The final large value of $k$ affects only $n+1$ modular binomial evaluations of order at most $n$. The algorithm never iterates through millions of days and never enumerates permutations for the target instance.

References

Problem page: Project Euler 996
Adjacent transposition: Wikipedia - Adjacent transposition
Degree sequence: Wikipedia - Degree in graph theory
Finite difference: Wikipedia - Finite difference
Dynamic programming: Wikipedia - Dynamic programming

Problem 996 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <set>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;

constexpr int TARGET_N = 123;
constexpr u64 TARGET_K = 4'567'891ULL;
constexpr u64 TARGET_E = TARGET_K / 2;
constexpr u64 MOD = 1'234'567'891ULL;

std::vector<std::vector<u64>> binomial_table(const int n, const int r, const u64 mod) {
    std::vector<std::vector<u64>> c(static_cast<std::size_t>(n + 1),
                                    std::vector<u64>(static_cast<std::size_t>(r + 1), 0));
    c[0][0] = 1;
    for (int i = 1; i <= n; ++i) {
        c[static_cast<std::size_t>(i)][0] = 1;
        const int hi = std::min(i, r);
        for (int j = 1; j <= hi; ++j) {
            const u64 value = c[static_cast<std::size_t>(i - 1)][static_cast<std::size_t>(j - 1)] +
                              c[static_cast<std::size_t>(i - 1)][static_cast<std::size_t>(j)];
            c[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] =
                mod == 0 ? value : value % mod;
        }
    }
    return c;
}

u64 choose_from_table(const std::vector<std::vector<u64>>& c, const int n, const int r) {
    if (n < 0 || r < 0 || n < r) {
        return 0;
    }
    if (n >= static_cast<int>(c.size()) || r >= static_cast<int>(c[0].size())) {
        return 0;
    }
    return c[static_cast<std::size_t>(n)][static_cast<std::size_t>(r)];
}

std::vector<std::vector<u64>> run_count_table(const int n, const int emax, const u64 mod) {
    const auto comb = binomial_table(2 * emax + 1, n, mod);
    std::vector<std::vector<u64>> runs(static_cast<std::size_t>(n + 1),
                                       std::vector<u64>(static_cast<std::size_t>(emax + 1), 0));

    for (int length = 2; length <= n; ++length) {
        for (int e = 1; e <= emax; ++e) {
            const u64 all_positive = choose_from_table(comb, 2 * e - 1, length - 1);
            const u64 too_large = choose_from_table(comb, e - 1, length - 1);
            if (mod == 0) {
                runs[static_cast<std::size_t>(length)][static_cast<std::size_t>(e)] =
                    all_positive - static_cast<u64>(length) * too_large;
            } else {
                const u64 bad = static_cast<u64>(length) * too_large % mod;
                runs[static_cast<std::size_t>(length)][static_cast<std::size_t>(e)] =
                    (all_positive + mod - bad) % mod;
            }
        }
    }

    return runs;
}

std::vector<u64> cumulative_counts(const int n, const int emax, const u64 mod) {
    const auto runs = run_count_table(n, emax, mod);
    std::vector<std::vector<u64>> total(static_cast<std::size_t>(n + 1),
                                        std::vector<u64>(static_cast<std::size_t>(emax + 1), 0));
    std::vector<std::vector<u64>> zero(static_cast<std::size_t>(n + 1),
                                       std::vector<u64>(static_cast<std::size_t>(emax + 1), 0));
    total[0][0] = 1;
    zero[0][0] = 1;

    for (int i = 1; i <= n; ++i) {
        zero[static_cast<std::size_t>(i)] = total[static_cast<std::size_t>(i - 1)];
        std::vector<u64> row = zero[static_cast<std::size_t>(i)];

        for (int length = 2; length <= i; ++length) {
            const auto& source = zero[static_cast<std::size_t>(i - length)];
            const auto& run = runs[static_cast<std::size_t>(length)];
            for (int before = 0; before <= emax; ++before) {
                const u64 source_value = source[static_cast<std::size_t>(before)];
                if (source_value == 0) {
                    continue;
                }
                for (int used = 1; before + used <= emax; ++used) {
                    const u64 run_value = run[static_cast<std::size_t>(used)];
                    if (run_value == 0) {
                        continue;
                    }
                    u64& cell = row[static_cast<std::size_t>(before + used)];
                    if (mod == 0) {
                        cell += source_value * run_value;
                    } else {
                        cell = (cell + source_value * run_value) % mod;
                    }
                }
            }
        }

        total[static_cast<std::size_t>(i)] = std::move(row);
    }

    std::vector<u64> cumulative(static_cast<std::size_t>(emax + 1), 0);
    u64 sum = 0;
    for (int e = 0; e <= emax; ++e) {
        if (mod == 0) {
            sum += total[static_cast<std::size_t>(n)][static_cast<std::size_t>(e)];
        } else {
            sum = (sum + total[static_cast<std::size_t>(n)][static_cast<std::size_t>(e)]) % mod;
        }
        cumulative[static_cast<std::size_t>(e)] = sum;
    }
    return cumulative;
}

u64 binomial_large_mod(const u64 n, int r) {
    if (r < 0 || n < static_cast<u64>(r)) {
        return 0;
    }
    if (static_cast<u64>(r) > n - static_cast<u64>(r)) {
        r = static_cast<int>(n - static_cast<u64>(r));
    }

    std::vector<u64> numerator(static_cast<std::size_t>(r));
    for (int i = 0; i < r; ++i) {
        numerator[static_cast<std::size_t>(i)] = n - static_cast<u64>(r) + 1 + static_cast<u64>(i);
    }

    for (int d = 2; d <= r; ++d) {
        u64 remaining = static_cast<u64>(d);
        for (u64& value : numerator) {
            const u64 g = std::gcd(value, remaining);
            if (g == 1) {
                continue;
            }
            value /= g;
            remaining /= g;
            if (remaining == 1) {
                break;
            }
        }
        assert(remaining == 1);
    }

    u64 result = 1;
    for (const u64 value : numerator) {
        result = result * (value % MOD) % MOD;
    }
    return result;
}

u64 polynomial_tail_value(const int n, const u64 e) {
    const int start = n - 1;
    const int degree = n;
    const int emax = start + degree;
    const std::vector<u64> cumulative = cumulative_counts(n, emax, MOD);

    std::vector<u64> current;
    current.reserve(static_cast<std::size_t>(degree + 1));
    for (int i = 0; i <= degree; ++i) {
        current.push_back(cumulative[static_cast<std::size_t>(start + i)]);
    }

    std::vector<u64> differences;
    differences.reserve(static_cast<std::size_t>(degree + 1));
    for (int r = 0; r <= degree; ++r) {
        differences.push_back(current[0]);
        std::vector<u64> next;
        next.reserve(current.size() - 1);
        for (std::size_t i = 0; i + 1 < current.size(); ++i) {
            next.push_back((current[i + 1] + MOD - current[i]) % MOD);
        }
        current = std::move(next);
    }

    const u64 x = e - static_cast<u64>(start);
    u64 result = 0;
    for (int r = 0; r <= degree; ++r) {
        result = (result + differences[static_cast<std::size_t>(r)] * binomial_large_mod(x, r)) % MOD;
    }
    return result;
}

u64 brute_count(const int n, const int emax) {
    using State = std::pair<std::vector<int>, std::vector<int>>;

    std::vector<int> identity(static_cast<std::size_t>(n));
    std::iota(identity.begin(), identity.end(), 0);

    std::set<State> states;
    states.insert({identity, std::vector<int>(static_cast<std::size_t>(n), 0)});

    std::set<std::vector<int>> results;
    results.insert(std::vector<int>(static_cast<std::size_t>(n), 0));

    for (int step = 1; step <= 2 * emax; ++step) {
        std::set<State> next;
        for (const auto& [perm, count] : states) {
            for (int p = 0; p + 1 < n; ++p) {
                std::vector<int> next_perm = perm;
                std::vector<int> next_count = count;
                ++next_count[static_cast<std::size_t>(next_perm[static_cast<std::size_t>(p + 1)])];
                std::swap(next_perm[static_cast<std::size_t>(p)], next_perm[static_cast<std::size_t>(p + 1)]);
                if (next_perm == identity) {
                    results.insert(next_count);
                }
                next.insert({std::move(next_perm), std::move(next_count)});
            }
        }
        states = std::move(next);
    }

    return static_cast<u64>(results.size());
}

void run_checkpoints() {
    for (int n = 2; n <= 5; ++n) {
        const std::vector<u64> values = cumulative_counts(n, 4, 0);
        for (int e = 0; e <= 4; ++e) {
            assert(values[static_cast<std::size_t>(e)] == brute_count(n, e));
        }
    }

    assert(cumulative_counts(3, 2, 0)[2] == 8);
    assert(cumulative_counts(12, 17, 0)[17] == 2'457'178'250ULL);

    for (int n = 2; n <= 8; ++n) {
        const int start = n - 1;
        const int degree = n;
        const int emax = start + degree + 5;
        const std::vector<u64> values = cumulative_counts(n, emax, 0);

        std::vector<u64> current;
        for (int i = 0; i <= degree; ++i) {
            current.push_back(values[static_cast<std::size_t>(start + i)]);
        }

        std::vector<u64> differences;
        for (int r = 0; r <= degree; ++r) {
            differences.push_back(current[0]);
            std::vector<u64> next;
            for (std::size_t i = 0; i + 1 < current.size(); ++i) {
                next.push_back(current[i + 1] - current[i]);
            }
            current = std::move(next);
        }

        for (int e = start; e <= emax; ++e) {
            u64 predicted = 0;
            u64 choose = 1;
            for (int r = 0; r <= degree; ++r) {
                if (r > 0) {
                    choose = choose * static_cast<u64>(e - start - r + 1) / static_cast<u64>(r);
                }
                predicted += differences[static_cast<std::size_t>(r)] * choose;
            }
            assert(predicted == values[static_cast<std::size_t>(e)]);
        }
    }
}

}  // namespace

int main() {
    run_checkpoints();
    std::cout << polynomial_tail_value(TARGET_N, TARGET_E) << '\n';
    return 0;
}

Python

from math import gcd


TARGET_N = 123
TARGET_K = 4_567_891
TARGET_E = TARGET_K // 2
MOD = 1_234_567_891


def binomial_table(n, r, mod):
    c = [[0] * (r + 1) for _ in range(n + 1)]
    c[0][0] = 1
    for i in range(1, n + 1):
        c[i][0] = 1
        hi = min(i, r)
        for j in range(1, hi + 1):
            value = c[i - 1][j - 1] + c[i - 1][j]
            c[i][j] = value if mod == 0 else value % mod
    return c


def choose_from_table(c, n, r):
    if n < 0 or r < 0 or n < r:
        return 0
    if n >= len(c) or r >= len(c[0]):
        return 0
    return c[n][r]


def run_count_table(n, emax, mod):
    comb = binomial_table(2 * emax + 1, n, mod)
    runs = [[0] * (emax + 1) for _ in range(n + 1)]

    for length in range(2, n + 1):
        for e in range(1, emax + 1):
            all_positive = choose_from_table(comb, 2 * e - 1, length - 1)
            too_large = choose_from_table(comb, e - 1, length - 1)
            if mod == 0:
                runs[length][e] = all_positive - length * too_large
            else:
                runs[length][e] = (all_positive - length * too_large) % mod
    return runs


def cumulative_counts(n, emax, mod):
    runs = run_count_table(n, emax, mod)
    total = [[0] * (emax + 1) for _ in range(n + 1)]
    zero = [[0] * (emax + 1) for _ in range(n + 1)]
    total[0][0] = 1
    zero[0][0] = 1

    for i in range(1, n + 1):
        zero[i] = total[i - 1][:]
        row = zero[i][:]

        for length in range(2, i + 1):
            source = zero[i - length]
            run = runs[length]
            for before, source_value in enumerate(source):
                if source_value == 0:
                    continue
                for used in range(1, emax - before + 1):
                    run_value = run[used]
                    if run_value == 0:
                        continue
                    if mod == 0:
                        row[before + used] += source_value * run_value
                    else:
                        row[before + used] = (row[before + used] + source_value * run_value) % mod

        total[i] = row

    cumulative = []
    running = 0
    for e in range(emax + 1):
        running += total[n][e]
        if mod:
            running %= mod
        cumulative.append(running)
    return cumulative


def binomial_large_mod(n, r):
    if r < 0 or n < r:
        return 0
    if r > n - r:
        r = n - r

    numerator = [n - r + 1 + i for i in range(r)]
    for d in range(2, r + 1):
        remaining = d
        for i, value in enumerate(numerator):
            g = gcd(value, remaining)
            if g == 1:
                continue
            numerator[i] = value // g
            remaining //= g
            if remaining == 1:
                break
        assert remaining == 1

    result = 1
    for value in numerator:
        result = result * (value % MOD) % MOD
    return result


def polynomial_tail_value(n, e):
    start = n - 1
    degree = n
    emax = start + degree
    cumulative = cumulative_counts(n, emax, MOD)

    current = [cumulative[start + i] for i in range(degree + 1)]
    differences = []
    for _ in range(degree + 1):
        differences.append(current[0])
        current = [(current[i + 1] - current[i]) % MOD for i in range(len(current) - 1)]

    x = e - start
    result = 0
    for r, diff in enumerate(differences):
        result = (result + diff * binomial_large_mod(x, r)) % MOD
    return result


def brute_count(n, emax):
    identity = tuple(range(n))
    states = {(identity, (0,) * n)}
    results = {(0,) * n}

    for _ in range(1, 2 * emax + 1):
        nxt = set()
        for perm, count in states:
            for p in range(n - 1):
                next_perm = list(perm)
                next_count = list(count)
                next_count[next_perm[p + 1]] += 1
                next_perm[p], next_perm[p + 1] = next_perm[p + 1], next_perm[p]
                next_perm_t = tuple(next_perm)
                next_count_t = tuple(next_count)
                if next_perm_t == identity:
                    results.add(next_count_t)
                nxt.add((next_perm_t, next_count_t))
        states = nxt

    return len(results)


def run_checkpoints():
    for n in range(2, 6):
        values = cumulative_counts(n, 4, 0)
        for e in range(5):
            assert values[e] == brute_count(n, e)

    assert cumulative_counts(3, 2, 0)[2] == 8
    assert cumulative_counts(12, 17, 0)[17] == 2_457_178_250

    for n in range(2, 9):
        start = n - 1
        degree = n
        emax = start + degree + 5
        values = cumulative_counts(n, emax, 0)

        current = [values[start + i] for i in range(degree + 1)]
        differences = []
        for _ in range(degree + 1):
            differences.append(current[0])
            current = [current[i + 1] - current[i] for i in range(len(current) - 1)]

        for e in range(start, emax + 1):
            predicted = 0
            choose = 1
            for r, diff in enumerate(differences):
                if r > 0:
                    choose = choose * (e - start - r + 1) // r
                predicted += diff * choose
            assert predicted == values[e]


def main():
    run_checkpoints()
    print(polynomial_tail_value(TARGET_N, TARGET_E))


if __name__ == "__main__":
    main()

Java

import java.util.Arrays;
import java.util.HashSet;
import java.util.Objects;
import java.util.Set;

public class Euler996 {
    static final int TARGET_N = 123;
    static final long TARGET_K = 4_567_891L;
    static final long TARGET_E = TARGET_K / 2;
    static final long MOD = 1_234_567_891L;

    static long[][] binomialTable(int n, int r, long mod) {
        long[][] c = new long[n + 1][r + 1];
        c[0][0] = 1;
        for (int i = 1; i <= n; i++) {
            c[i][0] = 1;
            int hi = Math.min(i, r);
            for (int j = 1; j <= hi; j++) {
                long value = c[i - 1][j - 1] + c[i - 1][j];
                c[i][j] = mod == 0 ? value : value % mod;
            }
        }
        return c;
    }

    static long chooseFromTable(long[][] c, int n, int r) {
        if (n < 0 || r < 0 || n < r) {
            return 0;
        }
        if (n >= c.length || r >= c[0].length) {
            return 0;
        }
        return c[n][r];
    }

    static long[][] runCountTable(int n, int emax, long mod) {
        long[][] comb = binomialTable(2 * emax + 1, n, mod);
        long[][] runs = new long[n + 1][emax + 1];

        for (int length = 2; length <= n; length++) {
            for (int e = 1; e <= emax; e++) {
                long allPositive = chooseFromTable(comb, 2 * e - 1, length - 1);
                long tooLarge = chooseFromTable(comb, e - 1, length - 1);
                if (mod == 0) {
                    runs[length][e] = allPositive - (long) length * tooLarge;
                } else {
                    runs[length][e] = Math.floorMod(allPositive - (long) length * tooLarge, mod);
                }
            }
        }
        return runs;
    }

    static long[] cumulativeCounts(int n, int emax, long mod) {
        long[][] runs = runCountTable(n, emax, mod);
        long[][] total = new long[n + 1][emax + 1];
        long[][] zero = new long[n + 1][emax + 1];
        total[0][0] = 1;
        zero[0][0] = 1;

        for (int i = 1; i <= n; i++) {
            zero[i] = Arrays.copyOf(total[i - 1], emax + 1);
            long[] row = Arrays.copyOf(zero[i], emax + 1);

            for (int length = 2; length <= i; length++) {
                long[] source = zero[i - length];
                long[] run = runs[length];
                for (int before = 0; before <= emax; before++) {
                    long sourceValue = source[before];
                    if (sourceValue == 0) {
                        continue;
                    }
                    for (int used = 1; before + used <= emax; used++) {
                        long runValue = run[used];
                        if (runValue == 0) {
                            continue;
                        }
                        if (mod == 0) {
                            row[before + used] += sourceValue * runValue;
                        } else {
                            row[before + used] = (row[before + used] + sourceValue * runValue) % mod;
                        }
                    }
                }
            }

            total[i] = row;
        }

        long[] cumulative = new long[emax + 1];
        long sum = 0;
        for (int e = 0; e <= emax; e++) {
            sum += total[n][e];
            if (mod != 0) {
                sum %= mod;
            }
            cumulative[e] = sum;
        }
        return cumulative;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = a % b;
            a = b;
            b = t;
        }
        return Math.abs(a);
    }

    static long binomialLargeMod(long n, int r) {
        if (r < 0 || n < r) {
            return 0;
        }
        if ((long) r > n - r) {
            r = (int) (n - r);
        }

        long[] numerator = new long[r];
        for (int i = 0; i < r; i++) {
            numerator[i] = n - r + 1L + i;
        }

        for (int d = 2; d <= r; d++) {
            long remaining = d;
            for (int i = 0; i < numerator.length; i++) {
                long g = gcd(numerator[i], remaining);
                if (g == 1) {
                    continue;
                }
                numerator[i] /= g;
                remaining /= g;
                if (remaining == 1) {
                    break;
                }
            }
            assert remaining == 1;
        }

        long result = 1;
        for (long value : numerator) {
            result = result * (value % MOD) % MOD;
        }
        return result;
    }

    static long polynomialTailValue(int n, long e) {
        int start = n - 1;
        int degree = n;
        int emax = start + degree;
        long[] cumulative = cumulativeCounts(n, emax, MOD);

        long[] current = new long[degree + 1];
        for (int i = 0; i <= degree; i++) {
            current[i] = cumulative[start + i];
        }

        long[] differences = new long[degree + 1];
        for (int r = 0; r <= degree; r++) {
            differences[r] = current[0];
            long[] next = new long[current.length - 1];
            for (int i = 0; i + 1 < current.length; i++) {
                next[i] = Math.floorMod(current[i + 1] - current[i], MOD);
            }
            current = next;
        }

        long x = e - start;
        long result = 0;
        for (int r = 0; r <= degree; r++) {
            result = (result + differences[r] * binomialLargeMod(x, r)) % MOD;
        }
        return result;
    }

    static long bruteCount(int n, int emax) {
        int[] identityArray = new int[n];
        for (int i = 0; i < n; i++) {
            identityArray[i] = i;
        }
        IntVector identity = new IntVector(identityArray);

        Set<State> states = new HashSet<>();
        states.add(new State(identity, new IntVector(new int[n])));

        Set<IntVector> results = new HashSet<>();
        results.add(new IntVector(new int[n]));

        for (int step = 1; step <= 2 * emax; step++) {
            Set<State> next = new HashSet<>();
            for (State state : states) {
                for (int p = 0; p + 1 < n; p++) {
                    int[] nextPerm = state.perm.values.clone();
                    int[] nextCount = state.count.values.clone();
                    nextCount[nextPerm[p + 1]]++;
                    int temp = nextPerm[p];
                    nextPerm[p] = nextPerm[p + 1];
                    nextPerm[p + 1] = temp;

                    IntVector permVector = new IntVector(nextPerm);
                    IntVector countVector = new IntVector(nextCount);
                    if (permVector.equals(identity)) {
                        results.add(countVector);
                    }
                    next.add(new State(permVector, countVector));
                }
            }
            states = next;
        }
        return results.size();
    }

    static void runCheckpoints() {
        for (int n = 2; n <= 5; n++) {
            long[] values = cumulativeCounts(n, 4, 0);
            for (int e = 0; e <= 4; e++) {
                assert values[e] == bruteCount(n, e);
            }
        }

        assert cumulativeCounts(3, 2, 0)[2] == 8;
        assert cumulativeCounts(12, 17, 0)[17] == 2_457_178_250L;

        for (int n = 2; n <= 8; n++) {
            int start = n - 1;
            int degree = n;
            int emax = start + degree + 5;
            long[] values = cumulativeCounts(n, emax, 0);

            long[] current = new long[degree + 1];
            for (int i = 0; i <= degree; i++) {
                current[i] = values[start + i];
            }

            long[] differences = new long[degree + 1];
            for (int r = 0; r <= degree; r++) {
                differences[r] = current[0];
                long[] next = new long[current.length - 1];
                for (int i = 0; i + 1 < current.length; i++) {
                    next[i] = current[i + 1] - current[i];
                }
                current = next;
            }

            for (int e = start; e <= emax; e++) {
                long predicted = 0;
                long choose = 1;
                for (int r = 0; r <= degree; r++) {
                    if (r > 0) {
                        choose = choose * (e - start - r + 1L) / r;
                    }
                    predicted += differences[r] * choose;
                }
                assert predicted == values[e];
            }
        }
    }

    public static void main(String[] args) {
        runCheckpoints();
        System.out.println(polynomialTailValue(TARGET_N, TARGET_E));
    }

    static final class IntVector {
        final int[] values;

        IntVector(int[] values) {
            this.values = values;
        }

        @Override
        public boolean equals(Object other) {
            return other instanceof IntVector && Arrays.equals(values, ((IntVector) other).values);
        }

        @Override
        public int hashCode() {
            return Arrays.hashCode(values);
        }
    }

    static final class State {
        final IntVector perm;
        final IntVector count;

        State(IntVector perm, IntVector count) {
            this.perm = perm;
            this.count = count;
        }

        @Override
        public boolean equals(Object other) {
            if (!(other instanceof State)) {
                return false;
            }
            State that = (State) other;
            return perm.equals(that.perm) && count.equals(that.count);
        }

        @Override
        public int hashCode() {
            return Objects.hash(perm, count);
        }
    }
}