【学习笔记】 常系数齐次线性递推

问题描述

给定一个 \(k\) 阶常系数齐次线性递推数列的前 \(k\) 项 \(h_1, h_2, h_3..h_k\) 和线性递推式 \(h_n = \sum_{i = 1} ^ k a_i h_{n - i}\)， 求这个数列的第 \(n\) 项。

复杂度要求： \(O(k^2logn)\)

加强： \(O(klogklogn)\)

前置知识

矩阵的特征值

设 \(A\) 是 \(n\) 阶方阵，如果存在数 \(\lambda\) 和非零 \(n\) 维列向量 \(x\)，使得 \(Ax=\lambda x\) 成立，则称 \(\lambda\) 是矩阵 \(A\) 的一个特征值。

矩阵的特征多项式

设 \(E\) 为单位矩阵，\(n\) 阶方阵 \(A\) 的特征多项式为 \(|\lambda E - A|\)

解释一下：\(\lambda\) 为一个变量。 \(\lambda E - A\) 的行列式是一个关于 \(\lambda\) 的 \(n\) 次多项式，即 \(A\) 的特征多项式。

特征值是 \(|\lambda E - A| = 0\) 的根。

Cayley-Hamilton定理

设 \(A\) 的特征多项式为 \(p_A\)，\(O\) 为零矩阵。用矩阵 \(A\) 代替 \(\lambda\) 带入 特征多项式，有 \(p_A(A) = O\)

简要证明一下：把 \(A\) 直接带进特征多项式， \(P_A(A) = |AE - A| = 0\)

其实我并不会证这个东西。自学的线代等于没学

算法流程

现在我们已经熟背了Cayley-Hamilton定理。

考虑计算一下矩阵快速幂的转移矩阵 \(M\) 的特征多项式。

设 \(M_{x,y}\) 为去掉位置 \((x,y)\) 的代数余子式，\(m_{x,y}\) 为 \((x,y)\) 位置上的元素，直接对第 \(1\) 行拉普拉斯展开，\(M = \sum_{i=1}^k m_{1,i} M_{1,i}\) ，可以发现去掉第一行第 \(i\) 列之后留下的是一个下三角矩阵，\(M_{1,i} = (-1) ^ {i + 1} (-1)^{i-1}\lambda^{k - i}\),\(m_{1,1} =\lambda-a_1,m_{1,i} = -a_i\)，整理一下就得到：

\[|E\lambda - M| = \lambda ^ k - a_1 \lambda ^ {k - 1} - a_2 \lambda ^ {k - 2} - ... - a_n \]

由 Cayley-Hamilton定理，我们得到 \(p_M(M) = O\)

现在我们用 \(n\) 代替 \(n - k\)，我们要求出 \(M^n\)

\(M^n= p_M(M) A(M) + r(M)\)，其中\(r(M)\) 的次数不高于\(k - 1\)

因为 \(p_M(M) = O\)， 得到 \(M^n = r(M)\)

我们现在只需要求 \(M^n \mod p_M(M)\)

因为 \(AB \mod C = (A \mod C)(B \mod C) \mod C\)，所以可以快速幂计算 \(M^n \mod p_M(M)\)

直接暴力复杂度 \(O(k^2logn)\), 用多项式乘法和取模可以做到 \(O(klogklogn)\)

现在得到了 \(M^n = r(M) = \sum_{i = 0} ^ {k - 1} c_i M ^ i\)

再分析一下就可以得到 \(ans = \sum_{i = 0} ^ {k - 1} c_i h_{i + k}\)

再暴力或者 NTT 处理一下前 \(h\) 的前 \(2k\) 项就好了。

需要注意，当 \(k = 1\) 时，\(M^1\) 需要对 \(p_M(M)\) 取模。

模板

BZOJ4161: Shlw loves matrixI

这道题里面 \(h\) 的下标是从 \(0\) 开始的。

#pragma GCC optimize("2,Ofast,inline")
#include<bits/stdc++.h>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define LL long long
#define pii pair<int, int>
using namespace std;
const int mod = 1e9 + 7;

template <typename T> T read(T &x) {
	int f = 0;
	register char c = getchar();
	while (c > '9' || c < '0') f |= (c == '-'), c = getchar();
	for (x = 0; c >= '0' && c <= '9'; c = getchar())
		x = (x << 3) + (x << 1) + (c ^ 48);
	if (f) x = -x;
	return x;
}

inline void upd(int &x, int y) {
	(x += y) >= mod ? x -= mod : 0;
}

inline int add(int x, int y) {
	return (x += y) >= mod ? x - mod : x;
}

inline int dec(int x, int y) {
	return (x -= y) < 0 ? x + mod : x;
}

namespace Linear {
	static const int Maxn = 5005;
    
	int n, k;
	int a[Maxn], h[Maxn];
	int b[Maxn], c[Maxn], p[Maxn];

	void module(int *x) {
		for (int i = k * 2; i >= k; --i) {
			int tmp = x[i];
			for (int j = 0; j <= k; ++j) {
				x[i - j] = dec(x[i - j], 1LL * p[k - j] * tmp % mod);
			}
		}
	}
	
	void mul(int *x, int *y, int *z) {
		static int res[Maxn];
		for (int i = 0; i <= k * 2; ++i) res[i] = 0;
		for (int i = 0; i < k; ++i) {
			for (int j = 0; j < k; ++j) {
				upd(res[i + j], 1LL * x[i] * y[j] % mod);       
			}
		}
		module(res);
		for (int i = 0; i < k; ++i) z[i] = res[i];
	}

	void poly_pow(int p) {
		while (p) {
			if (p & 1) mul(b, c, c);
			p >>= 1;
			if (!p) break;
			mul(b, b, b);
		}
	}
    
	int solve() {
		if (n <= k) return h[n];
		p[k] = 1;
		for (int i = 0; i < k; ++i)
			p[i] = dec(0, a[k - i]);
		b[1] = 1; c[0] = 1;
		if (k == 1) module(b);
		poly_pow(n - k);
		for (int i = k + 1; i <= k * 2; ++i) {
			for (int j = 1; j <= k; ++j) {
				upd(h[i], 1LL * h[i - j] * a[j] % mod);
			}
		}
		int ans = 0;
		for (int i = 0; i < k; ++i)
			upd(ans, 1LL * c[i] * h[i + k] % mod);
		return ans;
	}
}
using namespace Linear;

int main() {
	read(n); read(k); ++n;
	for (int i = 1; i <= k; ++i) {
		read(a[i]);
		if (a[i] < 0) a[i] += mod;
	}
	for (int i = 1; i <= k; ++i) {
		read(h[i]);
		if (h[i] < 0) h[i] += mod;
	}
	cout << solve() << endl;
	return 0;
}