CP Notebook

← all snippets

LinearRecurrence

Generates the k'th term of an n-order linear recurrence S[i] = Σ_j S[i-j-1]tr[j], given S[0 … ≥ n-1] and tr[0 … n-1]. Faster than matrix multiplication. Useful together with Berlekamp--Massey.

Time: O(n² log k) 26 lines bruteforce-tested mod 5 for n <= 5

Usage: linearRec(0, 1, 1, 1, k) // k'th Fibonacci number

content/numerical/LinearRecurrence.h — Lucian Bicsi, source: Chinese material

typedef vector<ll> Poly;
ll linearRec(Poly S, Poly tr, ll k) {
	int n = sz(tr);

	auto combine = [&](Poly a, Poly b) {
		Poly res(n * 2 + 1);
		rep(i,0,n+1) rep(j,0,n+1)
			res[i + j] = (res[i + j] + a[i] * b[j]) % mod;
		for (int i = 2 * n; i > n; --i) rep(j,0,n)
			res[i - 1 - j] = (res[i - 1 - j] + res[i] * tr[j]) % mod;
		res.resize(n + 1);
		return res;
	};

	Poly pol(n + 1), e(pol);
	pol[0] = e[1] = 1;

	for (++k; k; k /= 2) {
		if (k % 2) pol = combine(pol, e);
		e = combine(e, e);
	}

	ll res = 0;
	rep(i,0,n) res = (res + pol[i + 1] * S[i]) % mod;
	return res;
}