Detailed mathematical approach

Problem Summary

For integers \(n \ge 2\) and \(m \ge 1\), let \(T(n,m)\) be the number of length-\(m\) sequences \((a_1,\dots,a_m)\) such that

$1 \le a_i \le n-1 \quad (1 \le i \le m),\qquad a_i+a_{i+1} \le n \quad (1 \le i \le m-1).$

The task is to evaluate \(T(5000,10^{12})\) modulo \(10^9+7\). A direct dynamic program over all lengths up to \(10^{12}\) is impossible, so the solution first studies the finite transition system behind the neighbour rule, then reconstructs a short linear recurrence and jumps directly to the required term.

Mathematical Approach

The admissible sequences form a finite-state process on the values \(1,2,\dots,n-1\). Once that process is written in linear-algebra form, the distant term can be recovered from a much shorter recurrence.

Step 1: Encode the Neighbour Rule as a Transition Matrix

Use the last entry of the sequence as the state. If the current value is \(x\) and the next value is \(y\), the pair is valid exactly when

$x+y \le n.$

This gives a matrix \(A\in\{0,1\}^{(n-1)\times(n-1)}\) defined by

$A_{y,x}=\begin{cases} 1,& x+y\le n,\\ 0,& x+y\gt n. \end{cases}$

Each row is a consecutive block of ones followed by zeros. The first row contains \(n-1\) ones, the next row contains \(n-2\), and so on down to a single one. This triangular pattern is what makes the warm-up computation efficient.

Step 2: Write the Dynamic Programming Recurrence

Let \(f_t(y)\) be the number of valid sequences of length \(t\) ending at \(y\). For length \(1\), every state is allowed, so

$f_1(y)=1 \qquad (1\le y\le n-1).$

For \(t\ge 1\), the previous value can be any \(x\) with \(x\le n-y\), hence

$f_{t+1}(y)=\sum_{x=1}^{n-y} f_t(x).$

If we define prefix sums

$P_t(k)=\sum_{x=1}^{k} f_t(x),$

then the transition becomes

$f_{t+1}(y)=P_t(n-y).$

That identity explains the reverse-prefix update used by the implementation: one full layer can be produced in linear time. The total number of valid sequences of length \(t\) is

$T(n,t)=\sum_{y=1}^{n-1} f_t(y).$

In vector form, with \(F_t=(f_t(1),\dots,f_t(n-1))^T\) and \(\mathbf{1}\) the all-ones vector,

$F_t=A^{t-1}\mathbf{1},\qquad T(n,t)=\mathbf{1}^T A^{t-1}\mathbf{1}.$

Step 3: Why a Linear Recurrence Must Exist

The matrix \(A\) has size \(n-1\), so its minimal polynomial has degree at most \(n-1\). By the Cayley-Hamilton principle, every scalar sequence extracted from powers of \(A\), including \(T(n,t)\), must satisfy a linear recurrence of some order \(r\le n-1\):

$T(n,k+r)=c_1T(n,k+r-1)+c_2T(n,k+r-2)+\cdots+c_rT(n,k)\pmod{10^9+7}.$

This is the key structural fact. Instead of iterating the dynamic program all the way to \(m\), we only need enough initial terms to recover the recurrence and then evaluate the distant term quickly.

Step 4: Recover the Shortest Recurrence from Initial Terms

The implementation generates an initial segment

$T(n,1),T(n,2),T(n,3),\dots$

modulo \(10^9+7\), and then applies the Berlekamp-Massey algorithm to find the shortest linear recurrence that matches those values. If the recovered order is \(r\), its characteristic relation can be written as

$x^r-c_1x^{r-1}-c_2x^{r-2}-\cdots-c_r=0.$

Because \(r\le n-1\), any sample length comfortably above \(2(n-1)\) is enough for recovery; the implementations use \(2(n-1)+8\) warm-up terms.

Step 5: Compute the Distant Term by Polynomial Reduction

Once the recurrence is known, the problem becomes a standard fast nth-term computation. Reduce \(x^{m-1}\) modulo the recurrence polynomial:

$x^{m-1}\equiv q_0+q_1x+\cdots+q_{r-1}x^{r-1}\pmod{x^r-c_1x^{r-1}-\cdots-c_r}.$

Then the desired term is

$T(n,m)\equiv q_0T(n,1)+q_1T(n,2)+\cdots+q_{r-1}T(n,r)\pmod{10^9+7}.$

The coefficients \(q_i\) are obtained by binary exponentiation inside the quotient ring defined by the recurrence, so the huge exponent \(m=10^{12}\) contributes only a \(\log m\) factor.

Worked Example: \(n=3\)

Now the allowed values are \(1\) and \(2\). The neighbour rule forbids only the pair \((2,2)\), so

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

The totals by length are

$T(3,1)=2,\qquad T(3,2)=3,\qquad T(3,3)=5,\qquad T(3,4)=8.$

In this case the recurrence is simply

$T(3,t)=T(3,t-1)+T(3,t-2),$

so the checkpoint \(T(3,4)=8\) follows immediately.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. First they generate an initial block of totals using the end-state dynamic program. Instead of recomputing every state transition naively, they build prefix sums of the current layer and read those sums in reverse order, which turns one transition step into an \(O(n)\) operation.

Next they run Berlekamp-Massey modulo \(10^9+7\) on the warm-up totals and obtain the shortest recurrence for \(T(n,m)\). If the requested length already lies inside the precomputed prefix, the answer is returned directly. Otherwise the implementation performs binary exponentiation on reduced polynomials, producing the coefficients \(q_0,\dots,q_{r-1}\) and therefore the required distant term.

Complexity Analysis

Each warm-up dynamic-programming step costs \(O(n)\) time and \(O(n)\) memory. Since only \(O(n)\) initial terms are generated, the warm-up phase costs \(O(n^2)\) time overall. Berlekamp-Massey on a recurrence of order \(r\le n-1\) costs \(O(r^2)\), and the final nth-term stage costs \(O(r^2\log m)\). Therefore the total complexity is

$O(n^2+r^2+r^2\log m)=O(n^2\log m),$

with \(O(n)\) memory.

Footnotes and References

1. Project Euler problem page: Problem 654: Neighbourly Constraints
2. Berlekamp-Massey algorithm: Wikipedia - Berlekamp-Massey algorithm
3. Fast computation of linear recurrences: cp-algorithms - Linear Recurrence
4. Cayley-Hamilton theorem: Wikipedia - Cayley-Hamilton theorem
5. Prefix sums: Wikipedia - Prefix sum