Linear Recurrences

$f_n=c_1{f}_{n-1}+c_2{f}_{n-2}+\cdots +c_k{f}_{n-k}.$

Given the first $k$ terms and the coefficients, compute $f_n$ for enormous $n$.

The methods

Method	Time	When
Iterate	$O(nk)$	small $n$
Matrix exponentiation	$O(k^3\log n)$	$k\le 100$
Kitamasa	$O(k^2\log n)$	$k\le 5000$
Bostan-Mori	$O(k\log k\log n)$	$k$ large
Closed form (characteristic roots)	$O(1)$	distinct rational roots only

For $k=100$ and $n=10^{18}$: matrix exponentiation is $6\cdot 10^7$ operations, Kitamasa is $6\cdot 10^5$. Both are fine; for $k=5000$ only the latter two are.

The unifying idea

All three fast methods compute the same thing:

$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(\mathrm{mod}C(x))$, then

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

They differ only in how the polynomial multiplication is done: matrices ($k^3$), schoolbook polynomials ($k^2$), or FFT/NTT ($k\log k$).

Kitamasa

// f[0..k-1] initial terms, c[0..k-1] coefficients: f[n] = sum c[i] * f[n-1-i]
long long kitamasa(vector<long long> f, vector<long long> c, long long n, long long MOD) {
    int k = f.size();
    if (n < k) return f[n];
 
    // multiply two polynomials mod C(x)
    auto mulmod = [&](vector<long long> a, vector<long long> b) {
        vector<long long> r(2*k - 1, 0);
        for (int i = 0; i < k; i++)
            for (int j = 0; j < k; j++)
                r[i+j] = (r[i+j] + a[i] * b[j]) % MOD;
        for (int i = 2*k - 2; i >= k; i--)                  // reduce mod C(x)
            for (int j = 0; j < k; j++)
                r[i - 1 - j] = (r[i - 1 - j] + r[i] * c[j]) % MOD;
        r.resize(k);
        return r;
    };
 
    vector<long long> res(k, 0), base(k, 0);
    res[0] = 1;                                             // x^0
    if (k == 1) base[0] = c[0]; else base[1] = 1;           // x
    while (n) { if (n & 1) res = mulmod(res, base); base = mulmod(base, base); n >>= 1; }
 
    long long ans = 0;
    for (int i = 0; i < k; i++) ans = (ans + res[i] * f[i]) % MOD;
    return ans;
}

Finding the recurrence: Berlekamp-Massey

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

The combination is the point: brute-force the first ~200 terms, run Berlekamp-Massey to find the recurrence, then Kitamasa to jump to term $10^{18}$.

This works astonishingly often. Sequences that satisfy a linear recurrence include:

• 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 and permanents of banded matrices,
• many combinatorial counting sequences.

If your sequence has such a structure, you may not even need to understand it — just compute enough terms and let Berlekamp-Massey find the pattern.

Non-homogeneous recurrences

$f_n=\sum c_i{f}_{n-i}+g(n)$ with $g$ a polynomial of degree $d$: extend the state vector with $n^d,n^{d-1},\dots ,1$ and encode the binomial updates. The matrix grows to $(k+d+1)\times (k+d+1)$.

For $g(n)=r^n$, add $r^n$ to the state (one extra dimension, multiplied by $r$ each step).

Closed forms

If the characteristic polynomial $C(x)$ has distinct roots $r_1..r_k$, then

$f_n={\sum }_i{\alpha }_i{r}_i^n$

with the ${\alpha }_i$ determined by the initial conditions. A repeated root of multiplicity $m$ contributes $({\alpha }_0+{\alpha }_{1}n+\cdots +{\alpha }_{m-1}{n}^{m-1})r^n$.

Useful when the roots are nice (Fibonacci’s $\varphi$), but numerically dangerous for large $n$ and usually not worth it when a $\log n$ algorithm exists.

Sums of a linear recurrence

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

See also: Kitamasa · Berlekamp-Massey · Matrix Exponentiation