【学习笔记】 常系数齐次线性递推
问题描述
给定一个 \(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\)，整理一下就得到：
由 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;
}
