Solving Linear Recurrences

Given $f_n=c_1{f}_{n-1}+\cdots +c_k{f}_{n-k}$ and the first $k$ terms, compute $f_n$ for astronomically large $n$.

The unifying view

All fast methods compute the same object:

$x^n\mathrm{mod}C(x),C(x)=x^k-c_1{x}^{k-1}-\cdots -c_k,$

the characteristic polynomial. If $x^n\equiv {\sum }_{i<k}{a}_i{x}^i$, then

$f_n={\sum }_{i<k}{a}_i{f}_i.$

They differ only in how the polynomial multiplication modulo $C$ is performed.

Method	Multiplication	Time
Matrix exponentiation	matrices	$O(k^3\log n)$
Kitamasa	schoolbook polynomials	$O(k^2\log n)$
Bostan-Mori	NTT	$O(k\log k\log n)$

For $k=100$, $n=10^{18}$: matrices are $6\cdot 10^7$ operations, Kitamasa $6\cdot 10^5$, Bostan-Mori $\approx 10^5$.

Bostan-Mori — the fastest

The generating function of a linear recurrence is rational: $F(x)=\frac{P(x)}{Q(x)}$ with $\mathrm{deg}P<\mathrm{deg}Q=k$. To extract $[x^n]\frac{P}{Q}$:

$\frac{P(x)}{Q(x)}=\frac{P(x)Q(-x)}{Q(x)Q(-x)},$

and $Q(x)Q(-x)$ is a polynomial in $x^2$. So writing $V(x^2)=Q(x)Q(-x)$ and splitting $P(x)Q(-x)$ into even and odd parts $U_e(x^2)+x{U}_o(x^2)$:

$[x^n]\frac{P}{Q}=\{\begin{array}{ll}[x^{n/2}]\frac{U_e}{V}& n\text{ even}\\ [x^{(n-1)/2}]\frac{U_o}{V}& n\text{ odd}\end{array}$

Each step halves $n$ and costs two polynomial multiplications.

long long bostanMori(Poly P, Poly Q, long long n, long long MOD) {
    while (n) {
        Poly Qneg = Q;
        for (size_t i = 1; i < Qneg.size(); i += 2) Qneg[i] = (MOD - Qneg[i]) % MOD;
        Poly U = mul(P, Qneg), V = mul(Q, Qneg);
        P.assign((U.size() + 1 - (n & 1)) / 2, 0);
        for (size_t i = (n & 1); i < U.size(); i += 2) P[i / 2] = U[i];
        Q.assign((V.size() + 1) / 2, 0);
        for (size_t i = 0; i < V.size(); i += 2) Q[i / 2] = V[i];
        n >>= 1;
    }
    return P.empty() ? 0 : P[0] * powmod(Q[0], MOD - 2, MOD) % MOD;
}

$O(k\log k\log n)$ with NTT, $O(k^2\log n)$ with schoolbook (still competitive with Kitamasa and shorter to write).

Finding the recurrence: Berlekamp-Massey

Often you can compute the first few terms but do not know the recurrence. Berlekamp-Massey recovers the shortest linear recurrence from $2k$ terms in $O(k^2)$.

The combination is the technique: brute-force ~200 terms, run Berlekamp-Massey, then Bostan-Mori to reach $n=10^{18}$.

This works surprisingly often, because these sequences satisfy linear recurrences:

• counts of 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 combinatorial counting sequences.

You do not need to understand the recurrence — just verify it predicts terms beyond those it was fitted to.

vector<long long> rec = berlekampMassey(first200Terms);
long long ans = kitamasa(first200Terms, rec, n);

From a matrix to a recurrence

If the DP is $v_n=A{v}_{n-1}$ with $A$ a $k\times k$ matrix, its characteristic polynomial (degree $k$) gives a linear recurrence for every entry of $v_n$ by Cayley-Hamilton. Compute it with Faddeev-LeVerrier ($O(k^4)$) or the Hessenberg method ($O(k^3)$), then use Kitamasa — reducing $O(k^3\log n)$ to $O(k^2\log n)$.

For sparse $A$, Wiedemann finds the minimal polynomial in $O(k\cdot \mathrm{nnz})$.

Non-homogeneous terms

$f_n=\sum c_i{f}_{n-i}+g(n)$:

$g(n)$	Fix
constant $c$	extend the state with a permanent 1
polynomial of degree $d$	extend with $n^d,\dots ,n,1$ and encode the binomial updates
$r^n$	extend with $r^n$ (multiplied by $r$ each step)
$n^d{r}^n$	extend with all products

The recurrence order grows, but the method is unchanged.

Prefix sums

$S_n={\sum }_{i\le n}{f}_i$ also satisfies a linear recurrence of order $k+1$ — add $S$ to the state with $S_n=S_{n-1}+f_n$. So “sum the first $10^{18}$ terms” is the same problem.

See also: Bostan-Mori · Berlekamp-Massey · Linear Recurrences