x. Linear Recurrences

Given $f_n=c_1{f}_{n-1}+\cdots +c_k{f}_{n-k}$ and $f_0..f_{k-1}$, compute $f_n$ for $n\approx 10^{18}$; find the recurrence from terms (Berlekamp–Massey, and Reeds–Sloane for a composite modulus); and recurrences with polynomial coefficients (P-recursive sequences). Uses polyMul from NTT.

The unifying view. With the characteristic polynomial $C(x)=x^k-c_1{x}^{k-1}-\cdots -c_k$: if $x^n\equiv {\sum }_{i<k}{a}_i{x}^i(\mathrm{mod}C)$ then $f_n={\sum }_i{a}_i{f}_i$. All fast methods compute $x^n\mathrm{mod}C$; they differ in how they multiply.

Method	Multiplication	Time	When
iterate		$O(nk)$	small $n$
matrix power (Matrices)	matrices	$O(k^3\log n)$	$k\le 100$, or a general linear DP
Kitamasa	schoolbook mod $C$	$O(k^2\log n)$	$k\le 5000$, any modulus
Bostan–Mori	NTT on $P(x)/Q(x)$	$O(k\log k\log n)$	$k$ large
closed form $\sum {\alpha }_i{r}_i^n$	roots of $C$	$O(1)$	nice roots only; numerically dangerous

For $k=100$, $n=10^{18}$: matrices $6\cdot 10^7$ operations, Kitamasa $6\cdot 10^5$, Bostan–Mori about $10^5$. A root of multiplicity $m$ contributes $({\alpha }_0+{\alpha }_{1}n+\cdots +{\alpha }_{m-1}{n}^{m-1})r^n$ to the closed form.

The combination is the technique: brute-force about $2k+10$ terms, run Berlekamp–Massey, verify the recurrence predicts the terms it was not fitted to, then Bostan–Mori to $n=10^{18}$. It works because many sequences satisfy linear recurrences: walks in a fixed graph, tilings of a $k\times n$ board, any DP with a fixed finite state set and constant transitions, determinants of banded matrices, many counting sequences. A transition matrix gives one through its characteristic polynomial (Cayley–Hamilton, Characteristic Polynomial; Wiedemann for sparse ones).

Extension	How
$+$ constant / polynomial of degree $d$ / $r^n$ / $n^d{r}^n$	extend the state with $1$ / $n^d..1$ (binomial updates) / $r^n$ / all products; or just let Berlekamp–Massey find the longer recurrence
prefix sums $S_n={\sum }_{i\le n}{f}_i$	order $k+1$: multiply $Q$ by $1-x$
every $m$-th term $f_{mi+r}$	also linear recurrent (of order $\le k$): Graeffe-type step or BM on the subsequence
product $f_n{g}_n$, sum $f_n+g_n$	linear recurrent of order $\le k_1{k}_2$ / $k_1+k_2$

vector<mint> berlekampMassey(const vector<mint>& s) {      // shortest c with s[i] = sum_j c[j] s[i-1-j] for all i >= |c|
    int n = sz(s), L = 0, m = 0; vector<mint> C(n + 1, 0), B(n + 1, 0), T; C[0] = B[0] = 1; mint b = 1;
    for (int i = 0; i < n; i++) {
        m++; mint d = s[i]; for (int j = 1; j <= L; j++) d += C[j] * s[i - j];
        if (d == mint(0)) continue;
        T = C; mint coef = d / b;
        for (int j = m; j <= n; j++) C[j] -= coef * B[j - m];
        if (2 * L > i) continue;
        L = i + 1 - L; B = T; b = d; m = 0;
    }
    vector<mint> c(L); for (int j = 0; j < L; j++) c[j] = -C[j + 1];
    return c;
}
mint kitamasa(const vector<mint>& f, const vector<mint>& c, ll n) {   // f_i = sum_j c[j] f_{i-1-j}; O(k^2 log n)
    int k = sz(c); if (n < sz(f)) return f[n];
    if (k == 0) return 0;
    auto mulMod = [&](const vector<mint>& a, const vector<mint>& b) {
        vector<mint> r(2 * k, 0);
        for (int i = 0; i < k; i++) if (a[i] != mint(0)) for (int j = 0; j < k; j++) r[i + j] += a[i] * b[j];
        for (int i = 2 * k - 1; i >= k; i--) if (r[i] != mint(0)) for (int j = 1; j <= k; j++) r[i - j] += r[i] * c[j - 1];   // x^k = sum c_j x^{k-j}
        r.resize(k); return r;
    };
    vector<mint> res(k, 0), base(k, 0); res[0] = 1;
    if (k == 1) base[0] = c[0]; else base[1] = 1;
    for (; n; n >>= 1, base = mulMod(base, base)) if (n & 1) res = mulMod(res, base);
    mint ans = 0; for (int i = 0; i < k; i++) ans += res[i] * f[i];
    return ans;
}
mint bostanMori(Poly P, Poly Q, ll n) {                     // [x^n] P / Q, Q[0] != 0, O(k log k log n)
    while (n) {
        Poly Qn = Q; for (int i = 1; i < sz(Qn); i += 2) Qn[i] = -Qn[i];
        Poly U = polyMul(P, Qn), V = polyMul(Q, Qn);
        P.assign((sz(U) + 1) / 2 + 1, 0); Q.assign((sz(V) + 1) / 2, 0);
        for (int i = n & 1; i < sz(U); i += 2) P[i / 2] = U[i];
        for (int i = 0; i < sz(V); i += 2) Q[i / 2] = V[i];
        n >>= 1;
    }
    return P.empty() ? mint(0) : P[0] / Q[0];
}
mint linearRecurrenceNth(const vector<mint>& f, const vector<mint>& c, ll n) {   // via Bostan-Mori
    int k = sz(c); if (n < sz(f)) return f[n];
    Poly Q(k + 1); Q[0] = 1; for (int j = 0; j < k; j++) Q[j + 1] = -c[j];
    Poly P = polyMul(Poly(f.begin(), f.begin() + k), Q); P.resize(k);   // P = F Q mod x^k
    return bostanMori(P, Q, n);
}
mint guessNth(const vector<mint>& terms, ll n) { return linearRecurrenceNth(terms, berlekampMassey(terms), n); }

Proof: $x^n\mathrm{mod}C$, Kitamasa and Bostan–Mori

The key identity. Define the linear functional $\lambda (x^i)=f_i$ on polynomials. The recurrence says $\lambda (x^{j}C(x))=f_{j+k}-{\sum }_i{c}_i{f}_{j+k-i}=0$ for all $j\ge 0$, so $\lambda$ vanishes on the ideal $(C)$, and $f_n=\lambda (x^n)=\lambda (x^n\mathrm{mod}C)={\sum }_i{a}_i{f}_i$. Kitamasa computes $x^n\mathrm{mod}C$ by binary exponentiation with schoolbook products reduced by $x^k\equiv {\sum }_j{c}_j{x}^{k-j}$. Bostan–Mori. $F=\sum f_i{x}^i$ satisfies $FQ=P$ with $Q=1-{\sum }_j{c}_j{x}^j$ and $\mathrm{deg}P<k$: the coefficient of $x^m$ in $FQ$ for $m\ge k$ is $f_m-{\sum }_j{c}_j{f}_{m-j}=0$. Then $\frac{P}{Q}=\frac{P(x)Q(-x)}{Q(x)Q(-x)}$; $Q(x)Q(-x)=V(x^2)$ is even because it is invariant under $x\to -x$; splitting $P(x)Q(-x)=U_e(x^2)+x{U}_o(x^2)$ gives $[x^n]\frac{P}{Q}=[x^{\lfloor n/2\rfloor }]\frac{U_{n\mathrm{mod}2}}{V}$. The degrees stay $\le k$, each round is two multiplications and halves $n$, and at $n=0$ the coefficient is $P(0)/Q(0)$. Closure properties. A sequence is linear recurrent of order $\le k$ iff $f_n=u^T{A}^{n}v$ for some $k\times k$ matrix $A$ (companion matrix one way, Cayley–Hamilton the other). Then $f_{mi+r}=u^T(A^m{)}^i{A}^{r}v$ (order $\le k$), $f_n+g_n$ uses the block-diagonal $A\oplus B$ (order $\le k_1+k_2$), $f_n{g}_n=(u\otimes w{)}^T(A\otimes B{)}^n(v\otimes z)$ (order $\le k_1{k}_2$), and prefix sums satisfy the recurrence with characteristic polynomial $(x-1)C(x)$, i.e. $Q$ multiplied by $1-x$. Checked against direct iteration for random recurrences.

Proof: Berlekamp–Massey finds the shortest recurrence

Massey’s lemma. If a recurrence of length $L$ generates $s_0..s_{N-1}$ but not $s_N$, then every recurrence generating $s_0..s_N$ has length $L^{\prime }\ge N+1-L$. Suppose $s_m={\sum }_{i=1}^L{a}_i{s}_{m-i}$ for $L\le m\le N-1$, $s_m={\sum }_{j=1}^{L^{\prime }}{b}_j{s}_{m-j}$ for $L^{\prime }\le m\le N$, and $L+L^{\prime }\le N$. Then ${\sum }_i{a}_i{s}_{N-i}={\sum }_i{a}_i{\sum }_j{b}_j{s}_{N-i-j}$ (since $L^{\prime }\le N-i\le N$) $={\sum }_j{b}_j{\sum }_i{a}_i{s}_{N-j-i}={\sum }_j{b}_j{s}_{N-j}$ (since $L\le N-j\le N-1$) $=s_N$, a contradiction. The algorithm keeps the current connection polynomial $C$ (length $L$) and the polynomial $B$ from just before the last length change, whose discrepancy was $b$ at the position $m$ steps ago. If the new discrepancy $d\ne 0$, $C-\frac{d}{b}{x}^{m}B$ cancels it: the shifted $B$ produces zero discrepancy at all earlier positions and exactly $b$ at the current one. When $2L\le i$, the new length is $i+1-L$, which by the lemma is the minimum possible; otherwise the length stays $L$, still $\ge i+1-L$. By induction $L$ is always the minimal length for the prefix. $2k$ terms determine an order-$k$ recurrence uniquely (two recurrences of length $\le k$ agreeing on $2k$ terms generate the same sequence). $O(N^2)$. Verified on sequences generated by random recurrences (the found recurrence is no longer and reproduces all terms).

Reeds–Sloane: composite modulus

Berlekamp–Massey divides by discrepancies, which fails modulo a composite $m$. Reeds–Sloane (1985) factors $m=\prod p^e$; for each prime power it runs $e$ coupled Berlekamp–Massey instances (one per power of $p$ dividing the discrepancy, using only units for division), and the connection polynomials are glued with CRT. Give at least $2d$ terms for a recurrence of order $d$; cost $O(N^2\sum e)$. For a prime modulus use plain Berlekamp–Massey.

struct ReedsSloane {
    using vec = vector<ll>;
    static void extend(vec& a, size_t d, ll value = 0) { if (d > a.size()) a.resize(d, value); }
    static ll exgcd(ll a, ll b, ll& x, ll& y) { if (!b) { x = 1, y = 0; return a; } ll g = exgcd(b, a % b, y, x); y -= x * (a / b); return g; }
    static ll inverse(ll a, ll m) { ll x, y; return exgcd(a, m, x, y) == 1 ? (x % m + m) % m : -1; }
    static int L(const vec& a, const vec& b) {
        int da = (sz(a) > 1 || (sz(a) == 1 && a[0])) ? sz(a) - 1 : -40000, db = (sz(b) > 1 || (sz(b) == 1 && b[0])) ? sz(b) - 1 : -40000;
        return max(da, db + 1);
    }
    static pair<vec, vec> primePower(const vec& s, ll mod, ll p, int e) {   // mod = p^e
        vector<vec> a(e), b(e), an(e), bn(e), ao(e), bo(e);
        vec t(e), u(e), r(e), to(e, 1), uo(e), pw(e + 1, 1);
        for (int i = 1; i <= e; i++) pw[i] = pw[i - 1] * p;
        for (int i = 0; i < e; i++) {
            a[i] = {pw[i]}, an[i] = {pw[i]}, b[i] = {0}, bn[i] = {s[0] * pw[i] % mod}, t[i] = s[0] * pw[i] % mod;
            if (t[i] == 0) t[i] = 1, u[i] = e; else for (u[i] = 0; t[i] % p == 0; t[i] /= p, u[i]++) {}
        }
        for (size_t k = 1; k < s.size(); k++) {
            for (int g = 0; g < e; g++) if (L(an[g], bn[g]) > L(a[g], b[g])) { ao[g] = a[e - 1 - u[g]], bo[g] = b[e - 1 - u[g]], to[g] = t[e - 1 - u[g]], uo[g] = u[e - 1 - u[g]], r[g] = k - 1; }
            a = an, b = bn;
            for (int o = 0; o < e; o++) {
                ll d = 0; for (size_t i = 0; i < a[o].size() && i <= k; i++) d = (d + (i128)a[o][i] * s[k - i]) % mod;
                if (d == 0) { t[o] = 1, u[o] = e; continue; }
                for (u[o] = 0, t[o] = d; t[o] % p == 0; t[o] /= p, u[o]++) {}
                int g = e - 1 - u[o];
                if (L(a[g], b[g]) == 0) { extend(bn[o], k + 1); bn[o][k] = (bn[o][k] + d) % mod; }
                else {
                    ll coef = (i128)t[o] * inverse(to[g], mod) % mod * pw[u[o] - uo[g]] % mod; int m = k - r[g];
                    extend(an[o], ao[g].size() + m); extend(bn[o], bo[g].size() + m);
                    for (size_t i = 0; i < ao[g].size(); i++) an[o][i + m] = ((an[o][i + m] - (i128)coef * ao[g][i]) % mod + mod) % mod;
                    while (!an[o].empty() && an[o].back() == 0) an[o].pop_back();
                    for (size_t i = 0; i < bo[g].size(); i++) bn[o][i + m] = ((bn[o][i + m] - (i128)coef * bo[g][i]) % mod + mod) % mod;
                    while (!bn[o].empty() && bn[o].back() == 0) bn[o].pop_back();
                }
            }
        }
        return {an[0], bn[0]};
    }
    static pair<vec, int> connection(vec s, ll mod) {      // A (A[0] = 1) and L: sum_j A[j] s[i-j] = 0 for i >= L
        vector<tuple<ll, ll, int>> fac; ll m = mod;
        for (ll p = 2; p * p <= m; p++) if (m % p == 0) { ll pe = 1; int e = 0; while (m % p == 0) m /= p, pe *= p, e++; fac.emplace_back(pe, p, e); }
        if (m > 1) fac.emplace_back(m, m, 1);
        vector<vec> as; size_t n = 0; int len = 0;
        for (auto [pe, p, e] : fac) {
            vec ss = s; for (auto& x : ss) x = (x % pe + pe) % pe;
            auto [a, b] = primePower(ss, pe, p, e);
            len = max(len, L(a, b)); as.push_back(a); n = max(n, as.back().size());
        }
        vec res(n, 0);
        for (size_t i = 0; i < n; i++) {                    // CRT of the coefficients
            ll x = 0, M = 1;
            for (size_t j = 0; j < fac.size(); j++) { ll pe = get<0>(fac[j]), c = i < as[j].size() ? as[j][i] : 0; ll k = (c - x % pe + pe) % pe * inverse(M % pe, pe) % pe; x += M * k, M *= pe; }
            res[i] = x % mod;
        }
        len = max(len, 0); res.resize(max<size_t>(res.size(), len + 1), 0);
        return {res, len};
    }
};
vector<ll> reedsSloane(const vector<ll>& s, ll mod) {      // s[i] = sum_{j=1}^{d} c[j-1] s[i-j] (mod m) for i >= d = |c|
    if (mod == 1 || s.empty()) return {};
    auto [A, len] = ReedsSloane::connection(s, mod);
    vector<ll> c(len); for (int j = 1; j <= len; j++) c[j - 1] = (mod - A[j] % mod) % mod;
    return c;
}
ll linearRecurrenceMod(const vector<ll>& first, const vector<ll>& c, ll n, ll mod) {   // Kitamasa for any modulus
    int d = sz(c); if (n < sz(first)) return first[n] % mod;
    if (d == 0) return 0;
    auto mulMod = [&](const vector<ll>& a, const vector<ll>& b) {
        vector<ll> r(2 * d, 0);
        for (int i = 0; i < d; i++) if (a[i]) for (int j = 0; j < d; j++) r[i + j] = (r[i + j] + (i128)a[i] * b[j]) % mod;
        for (int i = 2 * d - 1; i >= d; i--) if (r[i]) for (int j = 1; j <= d; j++) r[i - j] = (r[i - j] + (i128)r[i] * c[j - 1]) % mod;
        r.resize(d); return r;
    };
    vector<ll> res(d, 0), base(d, 0); res[0] = 1 % mod;
    if (d == 1) base[0] = c[0] % mod; else base[1] = 1;
    for (ll e = n; e > 0; e >>= 1, base = mulMod(base, base)) if (e & 1) res = mulMod(res, base);
    ll ans = 0; for (int i = 0; i < d; i++) ans = (ans + (i128)res[i] * first[i]) % mod;
    return ans;
}

Proof sketch: Reeds–Sloane

Over $\mathbb{Z}/p^e$ the discrepancy $d$ may be a non-unit $p^{u}t$ with $t$ a unit. The algorithm keeps, for every $u<e$, the most recent connection pair $(a,b)$ whose discrepancy had exactly $p$-valuation $u$; a new discrepancy of valuation $u$ is cancelled with the stored pair of valuation $e-1-u$ scaled by $t{t}_{\text{old}}^{-1}{p}^{u-u_{\text{old}}}$, which only inverts units. Running the $e$ instances started from $p^i$ (i.e. $p^i$ times the trivial connection) handles all valuations at once; Reeds and Sloane prove the result is a shortest linear recurrence modulo $p^e$ (the $\mathbb{Z}/p^e$ analogue of Massey’s lemma). Different prime powers are independent and their connection polynomials combine coefficientwise by CRT, taking the maximal length. We verified the output on random recurrences modulo composite $m$ (e.g. $2^{10}$, $360$, $10^9+9\cdot 2$): the returned recurrence reproduces every given term and is never longer than the generating one, and linearRecurrenceMod matches iteration.

P-recursive sequences

${\sum }_{i=0}^r{P}_i(n)a_{n-i}=0$ with polynomial coefficients of degree $\le d$ (“holonomic”): factorials ($a_n-n{a}_{n-1}=0$), Catalan ($(n+1)a_n=(4n-2)a_{n-1}$), Motzkin, Schröder, central binomials, derangements, any algebraic generating function, sums ${\sum }_k$ of hypergeometric terms. Guess the relation from about $(r+1)(d+2)$ terms by solving for the $(r+1)(d+1)$ unknown coefficients (each $n\ge r$ gives one linear equation), then extend in $O(nrd)$; the $n$-th term for huge $n$ uses polynomial matrix products with sample shifting ($\tilde{O}(\sqrt{n})$, as for $n!\mathrm{mod}p$ in Evaluation and Interpolation).

vector<vector<mint>> findPRecursive(const vector<mint>& terms, int deg) {   // rel[i][e]: coefficient of n^e in P_i, at shift n - i
    int n = sz(terms), B = (n + 2) / (deg + 2), C = B * (deg + 1), R = n - (B - 1);
    if (B < 2 || R < C - 1) return {};
    vector<vector<mint>> mat(R, vector<mint>(C));
    for (int y = 0; y < R; y++) for (int b = 0; b < B; b++) { mint v = terms[y + b]; for (int e = 0; e <= deg; e++) mat[y][b * (deg + 1) + e] = v, v *= mint(y + b); }
    int rank = 0;
    for (int x = 0; x < C; x++) {
        int pivot = -1; for (int y = rank; y < R; y++) if (mat[y][x] != mint(0)) { pivot = y; break; }
        if (pivot < 0) break;
        swap(mat[rank], mat[pivot]); mint inv = mat[rank][x].inv();
        for (int x2 = x; x2 < C; x2++) mat[rank][x2] *= inv;
        for (int y = rank + 1; y < R; y++) if (mat[y][x] != mint(0)) { mint c = mat[y][x]; for (int x2 = x; x2 < C; x2++) mat[y][x2] -= c * mat[rank][x2]; }
        rank++;
    }
    if (rank == C) return {};
    for (int y = rank - 1; y >= 0; y--) if (mat[y][rank] != mint(0)) { mint c = mat[y][rank]; for (int y2 = 0; y2 < y; y2++) mat[y2][rank] -= c * mat[y2][y]; }
    int order = rank / (deg + 1); vector<vector<mint>> rel(order + 1, vector<mint>(deg + 1, 0));
    rel[0][rank % (deg + 1)] = 1;
    for (int y = rank - 1; y >= 0; y--) { int k = order - y / (deg + 1), e = y % (deg + 1); rel[k][e] = -mat[y][rank]; }
    return rel;
}
vector<mint> extendPRecursive(const vector<vector<mint>>& rel, vector<mint> terms, int n) {   // terms up to index n
    int order = sz(rel) - 1, deg = sz(rel[0]) - 1;
    while (sz(terms) <= n) {
        int m = sz(terms); mint s = 0;
        for (int i = 1; i <= order; i++) { int k = m - i; mint t = terms[k]; for (int e = 0; e <= deg; e++) s += t * rel[i][e], t *= mint(k); }
        mint den = 0, pw = 1; for (int e = 0; e <= deg; e++) den += pw * rel[0][e], pw *= mint(m);
        terms.push_back(-s / den);
    }
    terms.resize(n + 1); return terms;
}

Proof: guessing and extending P-recursive relations

Column block $b$ of the matrix holds $a_{y+b}(y+b{)}^e$, so a kernel vector gives polynomials $Q_b$ with ${\sum }_b{Q}_b(y+b)a_{y+b}=0$ for all rows $y$, i.e. a relation of order $B-1$ written with forward shifts; the code returns it as ${\sum }_i{P}_i(m-i)a_{m-i}=0$ with $P_0$ the leading block (normalised by its first free variable from the row echelon form, which the back substitution expresses through the pivots). With $R\ge C-1$ equations, a non-trivial kernel means the data satisfies some relation of the stated shape; it is the true relation if enough terms were given (verify on extra terms, as with Berlekamp–Massey). Extension: isolate $a_m=-\frac{{\sum }_{i\ge 1}{P}_i(m-i)a_{m-i}}{P_0(m)}$, valid while $P_0(m)\ne 0$. Checked on factorials, Catalan, Motzkin and central binomial coefficients (found from 20–30 terms, predictions verified to index 500).

Related notes

• Matrix Exponentiation DP
• Wiedemann (Berlekamp–Massey on Krylov sequences)
• Composition (Bostan–Mori in two variables)