Problem 654: Neighbourly Constraints

Project Euler Problem 654 Solution

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

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

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

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

Problem 654 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;

constexpr u64 kMod = 1'000'000'007ULL;

u64 mod_add(u64 a, u64 b) {
    a += b;
    if (a >= kMod) a -= kMod;
    return a;
}

u64 mod_sub(u64 a, u64 b) {
    return (a >= b) ? (a - b) : (a + kMod - b);
}

u64 mod_mul(u64 a, u64 b) {
    return (a * b) % kMod;
}

u64 mod_pow(u64 base, u64 exp) {
    u64 result = 1;
    base %= kMod;
    while (exp > 0) {
        if (exp & 1ULL) result = mod_mul(result, base);
        base = mod_mul(base, base);
        exp >>= 1ULL;
    }
    return result;
}

std::vector<u64> build_sequence(int n, int terms) {
    const int states = n - 1;
    std::vector<u64> vec(static_cast<std::size_t>(states), 1);
    std::vector<u64> prefix(static_cast<std::size_t>(states), 0);
    std::vector<u64> next_vec(static_cast<std::size_t>(states), 0);

    std::vector<u64> seq;
    seq.reserve(static_cast<std::size_t>(terms));

    for (int step = 1; step <= terms; ++step) {
        u64 total = 0;
        for (u64 x : vec) {
            total += x;
            if (total >= kMod) total -= kMod;
        }
        seq.push_back(total);

        u64 run = 0;
        for (int i = 0; i < states; ++i) {
            run += vec[static_cast<std::size_t>(i)];
            if (run >= kMod) run -= kMod;
            prefix[static_cast<std::size_t>(i)] = run;
        }

        for (int j = 0; j < states; ++j) {
            next_vec[static_cast<std::size_t>(j)] =
                prefix[static_cast<std::size_t>(states - 1 - j)];
        }
        vec.swap(next_vec);
    }

    return seq;
}

std::vector<u64> berlekamp_massey(const std::vector<u64>& s) {
    std::vector<u64> C{1}, B{1};
    int L = 0;
    int m = 1;
    u64 b = 1;

    for (int n = 0; n < static_cast<int>(s.size()); ++n) {
        u64 d = s[static_cast<std::size_t>(n)];
        for (int i = 1; i <= L; ++i) {
            d = (d + mod_mul(C[static_cast<std::size_t>(i)], s[static_cast<std::size_t>(n - i)])) % kMod;
        }

        if (d == 0) {
            ++m;
            continue;
        }

        const std::vector<u64> T = C;
        const u64 coef = mod_mul(d, mod_pow(b, kMod - 2));

        if (C.size() < B.size() + static_cast<std::size_t>(m)) {
            C.resize(B.size() + static_cast<std::size_t>(m), 0);
        }

        for (std::size_t i = 0; i < B.size(); ++i) {
            const std::size_t idx = i + static_cast<std::size_t>(m);
            C[idx] = mod_sub(C[idx], mod_mul(coef, B[i]));
        }

        if (2 * L <= n) {
            L = n + 1 - L;
            B = T;
            b = d;
            m = 1;
        } else {
            ++m;
        }
    }

    std::vector<u64> rec(static_cast<std::size_t>(L), 0);
    for (int i = 1; i <= L; ++i) {
        rec[static_cast<std::size_t>(i - 1)] = mod_sub(0, C[static_cast<std::size_t>(i)]);
    }
    return rec;
}

std::vector<u64> combine_poly(const std::vector<u64>& a,
                              const std::vector<u64>& b,
                              const std::vector<u64>& rec) {
    const int k = static_cast<int>(rec.size());
    std::vector<u64> tmp(static_cast<std::size_t>(2 * k - 1), 0);

    for (int i = 0; i < k; ++i) {
        if (a[static_cast<std::size_t>(i)] == 0) continue;
        for (int j = 0; j < k; ++j) {
            if (b[static_cast<std::size_t>(j)] == 0) continue;
            const u64 add =
                mod_mul(a[static_cast<std::size_t>(i)], b[static_cast<std::size_t>(j)]);
            u64& cell = tmp[static_cast<std::size_t>(i + j)];
            cell += add;
            if (cell >= kMod) cell -= kMod;
        }
    }

    for (int i = 2 * k - 2; i >= k; --i) {
        const u64 x = tmp[static_cast<std::size_t>(i)];
        if (x == 0) continue;
        for (int t = 1; t <= k; ++t) {
            const int idx = i - t;
            const u64 add = mod_mul(x, rec[static_cast<std::size_t>(t - 1)]);
            u64& cell = tmp[static_cast<std::size_t>(idx)];
            cell += add;
            if (cell >= kMod) cell -= kMod;
        }
    }

    tmp.resize(static_cast<std::size_t>(k));
    return tmp;
}

u64 linear_recurrence_nth(const std::vector<u64>& init,
                          const std::vector<u64>& rec,
                          u64 n_one_based) {
    const int k = static_cast<int>(rec.size());
    if (n_one_based <= static_cast<u64>(k)) return init[static_cast<std::size_t>(n_one_based - 1)];

    u64 exp = n_one_based - 1;
    std::vector<u64> result(static_cast<std::size_t>(k), 0);
    result[0] = 1;

    std::vector<u64> base(static_cast<std::size_t>(k), 0);
    if (k == 1) {
        base[0] = rec[0];
    } else {
        base[1] = 1;
    }

    while (exp > 0) {
        if (exp & 1ULL) result = combine_poly(result, base, rec);
        exp >>= 1ULL;
        if (exp > 0) base = combine_poly(base, base, rec);
    }

    u64 ans = 0;
    for (int i = 0; i < k; ++i) {
        ans = (ans + mod_mul(result[static_cast<std::size_t>(i)], init[static_cast<std::size_t>(i)])) % kMod;
    }
    return ans;
}

u64 solve_case(int n, u64 m) {
    const int need_terms = 2 * (n - 1) + 8;
    const std::vector<u64> seq = build_sequence(n, need_terms);
    if (m <= seq.size()) return seq[static_cast<std::size_t>(m - 1)];

    const std::vector<u64> rec = berlekamp_massey(seq);
    const int k = static_cast<int>(rec.size());
    std::vector<u64> init(seq.begin(), seq.begin() + k);
    return linear_recurrence_nth(init, rec, m);
}

}  // namespace

int main() {
    assert(solve_case(3, 4) == 8);
    assert(solve_case(5, 5) == 246);
    assert(solve_case(10, 100) == 862820094);
    assert(solve_case(100, 10) == 782136797);

    std::cout << solve_case(5000, 1'000'000'000'000ULL) << "\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

import java.util.ArrayList;

public class Euler654 {

    static final long kMod = 1000000007L;

    static long modAdd(long a, long b) {
        a += b;
        return (a >= kMod) ? (a - kMod) : a;
    }

    static long modSub(long a, long b) {
        return (a >= b) ? (a - b) : (a + kMod - b);
    }

    static long modMul(long a, long b) {
        return (a * b) % kMod;
    }

    static long modPow(long base, long exp) {
        long result = 1;
        base %= kMod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                result = modMul(result, base);
            base = modMul(base, base);
            exp >>= 1;
        }
        return result;
    }

    static ArrayList<Long> buildSequence(int n, int terms) {
        int states = n - 1;
        long[] vec = new long[states];
        long[] prefix = new long[states];
        long[] nextVec = new long[states];
        for (int i = 0; i < states; i++)
            vec[i] = 1;

        ArrayList<Long> seq = new ArrayList<>(terms);

        for (int step = 1; step <= terms; ++step) {
            long total = 0;
            for (int i = 0; i < states; ++i) {
                total += vec[i];
            }
            total %= kMod;
            seq.add(total);

            long run = 0;
            for (int i = 0; i < states; ++i) {
                run += vec[i];
                if (run >= kMod)
                    run -= kMod;
                prefix[i] = run;
            }

            for (int j = 0; j < states; ++j) {
                nextVec[j] = prefix[states - 1 - j];
            }
            long[] tmp = vec;
            vec = nextVec;
            nextVec = tmp;
        }

        return seq;
    }

    static ArrayList<Long> berlekampMassey(ArrayList<Long> s) {
        ArrayList<Long> C = new ArrayList<>();
        ArrayList<Long> B = new ArrayList<>();
        C.add(1L);
        B.add(1L);
        int L = 0;
        int m = 1;
        long b = 1;

        for (int n = 0; n < s.size(); ++n) {
            long d = s.get(n);
            for (int i = 1; i <= L; ++i) {
                d = (d + modMul(C.get(i), s.get(n - i))) % kMod;
            }

            if (d == 0) {
                ++m;
                continue;
            }

            ArrayList<Long> T = new ArrayList<>(C);
            long coef = modMul(d, modPow(b, kMod - 2));

            while (C.size() < B.size() + m) {
                C.add(0L);
            }

            for (int i = 0; i < B.size(); ++i) {
                int idx = i + m;
                C.set(idx, modSub(C.get(idx), modMul(coef, B.get(i))));
            }

            if (2 * L <= n) {
                L = n + 1 - L;
                B = T;
                b = d;
                m = 1;
            } else {
                ++m;
            }
        }

        ArrayList<Long> rec = new ArrayList<>(L);
        for (int i = 1; i <= L; ++i) {
            rec.add(modSub(0, C.get(i)));
        }
        return rec;
    }

    static long[] combinePoly(long[] a, long[] b, long[] rec) {
        int k = rec.length;
        long[] tmp = new long[2 * k - 1];

        for (int i = 0; i < a.length; ++i) {
            if (a[i] == 0)
                continue;
            for (int j = 0; j < b.length; ++j) {
                if (b[j] == 0)
                    continue;
                tmp[i + j] = (tmp[i + j] + modMul(a[i], b[j])) % kMod;
            }
        }

        for (int i = 2 * k - 2; i >= k; --i) {
            long x = tmp[i];
            if (x == 0)
                continue;
            for (int t = 1; t <= k; ++t) {
                int idx = i - t;
                tmp[idx] = (tmp[idx] + modMul(x, rec[t - 1])) % kMod;
            }
        }

        long[] ret = new long[k];
        System.arraycopy(tmp, 0, ret, 0, k);
        return ret;
    }

    static long linearRecurrenceNth(ArrayList<Long> initList, ArrayList<Long> recList, long nOneBased) {
        int k = recList.size();
        if (nOneBased <= k)
            return initList.get((int) (nOneBased - 1));

        long[] rec = new long[k];
        long[] init = new long[k];
        for (int i = 0; i < k; i++) {
            rec[i] = recList.get(i);
            init[i] = initList.get(i);
        }

        long exp = nOneBased - 1;
        long[] result = new long[k];
        result[0] = 1;

        long[] base = new long[k];
        if (k == 1) {
            base[0] = rec[0];
        } else {
            base[1] = 1;
        }

        while (exp > 0) {
            if ((exp & 1) != 0) {
                result = combinePoly(result, base, rec);
            }
            exp >>= 1;
            if (exp > 0) {
                base = combinePoly(base, base, rec);
            }
        }

        long ans = 0;
        for (int i = 0; i < k; ++i) {
            ans = (ans + modMul(result[i], init[i])) % kMod;
        }
        return ans;
    }

    static long solveCase(int n, long m) {
        int needTerms = 2 * (n - 1) + 8;
        ArrayList<Long> seq = buildSequence(n, needTerms);
        if (m <= seq.size())
            return seq.get((int) (m - 1));

        ArrayList<Long> rec = berlekampMassey(seq);
        int k = rec.size();
        ArrayList<Long> init = new ArrayList<>(seq.subList(0, k));
        return linearRecurrenceNth(init, rec, m);
    }

    public static String solve() {
        long ans = solveCase(5000, 1000000000000L);
        return Long.toString(ans);
    }

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