特征多项式与常系数齐次线性递推
定义
对于 \(n\) 阶矩阵，若存在非零列向量 \(x\) 和数 \(\lambda\) 满足 \(Ax=\lambda x\)，则称 \(\lambda\) 和 \(x\) 为一组对应的特征值和特征向量。
在确定了特征值之后，可以得到对应的无穷多个解。
求解特征值和特征向量
求解特征值和特征向量：
容易发现，\(\lambda\) 是一个特征值，只需要满足 \(Ax=\lambda x\) 有解，以 \(x\) 为元容易列出方程，其常数项为均 \(0\)，系数矩阵为。
其中\(I\)是单位矩阵。
这个方程有非零解的充要条件是：\(det(\lambda I-A)=0\) （因为如果不为 \(0\)，则矩阵满秩，所有向量线性无关，无法得到0向量）
而 \(det(\lambda I-A)\) 是一个 \(n\) 次多项式 \(p(\lambda)\)，称为特征多项式，所有的特征值 \(\lambda\) 就是 \(p(\lambda)\) 的根。
应用
加速矩阵乘法：
由 \(Ax=\lambda x\)，迭代该式可以得到 \(A^nx=\lambda^nx\) 。
特殊矩阵的特征值
上三角矩阵
带入行列式即可知道 \(\displaystyle det(\lambda I-A)=\prod (\lambda -A_{i,i})\) 。
也就是说，主对角线上所有的 \(A_{i,i}\) 都是 \(det(\lambda I-A)=0\) 的根。
零化多项式
对于一个矩阵 \(A\)，它的一个零化多项式 \(\varphi(\lambda)\) 是满足 \(\varphi(A)=0\) 的多项式，定义域包含矩阵。
最小多项式：次数最低的零化多项式。
\(\text{Cayley-Hamilton}\) 定理
特征多项式：\(p(\lambda)=|\lambda I-A|\)，\(\lambda\) 定义域不止是 \(\mathbb{ R }\)，还可以是矩阵。
\(\text{Cayley-Hamilton}\) 定理指出：矩阵的特征多项式也是它的零化多项式。
即令：
则有：
证明：
将 \(\varphi(\lambda)\) 改写为：
由定理：任意的 \(n\) 阶矩阵都能相似为上三角矩阵 可知，存在可逆矩阵 \(P\)，使得：
将 \(PAP^{-1}\) 代入 \(\varphi(\lambda)\) 得到：
计算：
即 \(\varphi(PAP^{-1})=P\varphi(A)P^{-1}=O\)，故有 \(\varphi(A)=O\)
求解特征多项式
带入 \(n\) 个数，求出得 \(det(x I_n-A)\)，得到 \(n\) 个矩阵，通过高斯消元可以 \(O(n^3)\) 地求出行列式。
然后可 \(O(n^2)\) 拉格朗日插值求出原来的多项式，总复杂度受限于高斯消元，为 \(O(n^4)\) 。
求解最小多项式
构造矩阵序列 \(a_i=A^i\)。
求出它的一个线性递推 \(r_i\)，即：
所以可以由 \(r_i\) 翻转得到 \(f(\lambda)\) 。
求解 \(a_i\) 前 \(n\) 项的复杂度受限于矩阵乘法为 \(O(n^4)\)，求解递推式的复杂度为 \(O(n^3)\) 。
考虑到实际求解递推式时，随机生成了两个向量 \(u,v\) 。
实际是计算标量序列 \(\{uA^iv\}\) 的递推式，所以实际每次求出 \(uA^i\) 复杂度应为 \(O(n^2)\) 。
求这个递推式需要用到 \(a_i\) 前 \(2n\) 项，求解复杂度为 \(O(n^3)\) 。
因此总复杂度为 \(O(n^3)\) 。
(但是如果只是求出来并没有什么用，因为求解方法是随机的，甚至连检查一次保证正确都需要 \(O(n^2(n+e))\) 的时间(\(e\) 为矩阵非 \(0\) 位置个数))
求解稀疏方程组
设方程系数用矩阵 \(A\) 表示，右侧每个方程的常数用向量 \(b\) 表示，答案用向量 \(x\) 表示，则满足关系式
\(Ax=b\)，即 \(x=A^{-1}b\) 。
求出 \(\{A^ib\}\) 线性递推式，反推出 \(A^{-1}b\) 即可。
反推方法：
带入线性递推的 \(m\) 项，则 \(\sum_{i=0}^{m}\limits A^{m-i}b\cdot r_i=0\)
两边同乘 \(A^{-1}\)，得到 \(A^{-1}b\cdot r_m +\sum_{i=0}^{m-1}\limits A^{m-i}br_i=0\)
求解矩阵 k 次幂
我们要求解 \(A^k\)，常规做法是直接用快速幂
设矩阵 \(A\) 的一个零化多项式是 \(f(\lambda)\) ，显然，\(A^k\) 可以用一个多项式表示 \(A^k=\sum_0^k w_i A^i\)。
\(\{w_i\}\) 构成了一个 \(k+1\) 次多项式 \(F_k(x)\)。
存在一种合法的表示是 \(F_k(x)=x^k\)。
\( \because f(A)=0\\ \therefore \forall i, f(A)A^i=0\\ \)
所以对于任意实数 \(T\)，\(G_k(x)=x^k-Tf(x)\) 也合法，也就是相当于我们要求出 \(x^k\) 对于 \(f(x)\) 这个 \(n+1\) 多项式取模。
显然可以通过类似快速幂的方式倍增求解这个多项式，每次对 \(f(x)\) 取模复杂度是 \(O(n\log n)\) ，总时间复杂度 \(O(n\log m\log n)\) 。
最后得到的 \(F(x)\) 是一个 \(n\) 次多项式，带入就可以快速求出 \(A_k\)，可以认为这个复杂度是受限于求解 \(A^0,A^1,\cdots,A^{n-1}\) 的 \(O(n^4)\) 。
对于元矩阵 \(A\) 为稀疏矩阵的情况，设其包含 \(e\) 个非零位置，那么求解 \(B\cdot A\) 的过程是 \(O(n\cdot e)\) 的，求解 \(A_0,A^1,\cdots,A^{n-1}\) 的过程，是 \(O(n^2e)\) 的。
求解零化多项式的复杂度也是 \(O(n^2(n+e))\) 的，因此总复杂度为 \(O(n^2(n+e))\)。
而一般的矩阵快速幂是 \(O(n^3\log k)\) 的，这种方法适用情况非常特殊。
另外，对于并不需要知道整个矩阵的答案，并且 \(A^0,A^1,\cdots,A^{n-1}\) 特殊的具体问题，这个方法也十分有效。
求解常系数线性齐次递推
问题是要求数列 \(f_i=\sum _{j=1}^{n}a_j\cdot f_{i-j}\) 。
给出 \(f_0,f_1,\cdots,f_{n-1}\)，求第 \(k\) 项的值。
线性递推显然可以用初始向量列与转移矩阵的幂次的乘积表示，即\(f_i=(S \cdot A^i)_n\)，其中\(A\)为转移矩阵，\(S\)为初始向量列，我们求的是第\(n\)项。
对于\(n=4\)的情况，我们的转移矩阵。
鉴于它的特殊性，我们可以直接求出它的特征多项式表达式。
由
带入行列式最暴力的求法
枚举一个排列\(p_i\)，设排列\(p\)的逆序对为\(f(p)\)，\(|A|=\sum (-1)^{f(p)} \Pi A_{i,p_i}\)
实际上合法的排列只有\(n\)个，枚举 \(p_i=n\)
那么：
当 \(i=n\) 时，\((-1)^{f(p)} \prod A_{i,p_i}=\lambda ^n-a_1\lambda ^{n-1}\)
当 \(i\neq n\) 时：
综上，转移矩阵 \(A\) 的特征多项式有简单的表达
假设有 \(f_0\) 这一项(不需要知道是多少)，那么认为初始向量列为 \(S=(f_{-(n-1)},f_{-(n-2)},\cdots ,f_{0})\) 。
这个问题，我们要求的是 \(S\cdot A^k\) 的第 \(n\) 项，不需要知道整个矩阵。
类似求出 \(A^k\) 的过程，求出 \(F_k(x)\mod p(\lambda)\)。
我们要求解 \((S\cdot A^k)_n=\sum_1^{n}[x^i]{F(x)}(S\cdot A^i)_n\)，而 \((S\cdot A^i)_n=f_i\) 已知，求出 \(F(x)\) 后直接带入即可
需要用到多项式取模，求解这个表达式是 \(O(n\log n\log k)\) 的，求完直接带入即可。
模板代码
#4161. Shlw loves matrixI
非 \(\text{NTT}\) 模数，\(k\le2000\) ，可以 \(O(k^2)\) 卷积，总复杂度 \(O(k^2\log n)\)。
点击查看代码
#include <bits/stdc++.h>
using namespace std;
int n, m;
const int md = 1e9 + 7;
inline int pwr(int x, int y) {
int res = 1;
while (y) {
if (y & 1)
res = 1ll * res * x % md;
x = 1ll * x * x % md;
y >>= 1;
}
return res;
}
struct Poly {
vector<int> dp;
Poly(int x = 0) { dp.resize(x); }
inline auto size() { return dp.size(); }
inline int& operator[](int t) { return dp[t]; }
inline Poly Reverse() {
Poly res = *this;
reverse(res.dp.begin(), res.dp.end());
return res;
}
inline Poly operator-(Poly b) {
Poly res(max(dp.size(), b.size()));
for (int i = 0; i < dp.size(); i++) res[i] = (res[i] + dp[i] + md) % md;
for (int i = 0; i < b.size(); i++) res[i] = (res[i] - b[i] + md) % md;
return res;
}
inline Poly operator*(int x) {
Poly res = *this;
for (auto& it : res.dp) it = 1ll * it * x % md;
return res;
}
inline Poly operator*(Poly y) {
Poly res(dp.size() + y.size() - 1);
for (int i = 0; i < dp.size(); i++) {
for (int j = 0; j < y.size(); j++) res[i + j] = (res[i + j] + 1ll * dp[i] * y[j]) % md;
}
return res;
}
inline Poly Inv(int len) {
if (len == 1) {
Poly res(1);
res[0] = pwr(dp[0], md - 2);
return res;
}
Poly res = this->Inv((len + 1) / 2), tmp = *this;
tmp.dp.resize(len);
res = res * 2 - tmp * res * res;
res.dp.resize(len);
return res;
}
inline Poly operator/(Poly b) {
int N = dp.size(), M = b.size();
if (N < M)
return Poly();
Poly res = this->Reverse() * (b.Reverse().Inv(N - M + 1));
res.dp.resize(N - M + 1);
return res.Reverse();
}
inline Poly operator%(Poly b) {
Poly G = *this / b, res = *this - (G * b);
res.dp.resize(min(dp.size(), b.size()));
return res;
}
};
inline Poly Polypower(Poly x, int y, Poly MD) {
Poly res(1);
res[0] = 1;
while (y) {
if (y & 1)
res = res * x % MD;
x = x * x % MD;
y >>= 1;
}
return res;
}
int a[40005], k;
int main() {
scanf("%d%d", &k, &n);
Poly F(n + 1), X(2);
F[n] = 1;
X[1] = 1;
for (int i = n - 1; ~i; i--) scanf("%d", &F[i]), F[i] = (md - F[i]) % md;
for (int i = 0; i < n; i++) scanf("%d", &a[i]);
Poly tmp = Polypower(X, k, F);
int res = 0;
for (int i = 0; i < n; i++) res = (res + 1ll * a[i] * tmp[i]) % md;
printf("%d\n", (res + md) % md);
return 0;
}
P4723 【模板】常系数齐次线性递推
\(\text{NTT}\) 模数，\(k\le32000\) ，\(O(k\log k)\) 卷积，总复杂度 \(O(k\log k\log n)\)。
点击查看代码
#include <bits/stdc++.h>
using namespace std;
int n, m;
const int md = 998244353, G = 3, Gi = (md + 1) / 3;
int r[1 << 20];
inline int pwr(int x, int y) {
int res = 1;
while (y) {
if (y & 1)
res = 1ll * res * x % md;
x = 1ll * x * x % md;
y >>= 1;
}
return res;
}
inline void NTT(vector<int> &dp, int lim, int W) {
for (int i = 0; i < (1 << lim); i++)
if (i < r[i])
swap(dp[i], dp[r[i]]);
for (int i = 0; i < lim; i++) {
int w = pwr(W, (md - 1) / (1 << (i + 1)));
for (int j = 0; j < (1 << lim); j += (1 << (i + 1))) {
int Pw = 1;
for (int t = 0; t < (1 << i); t++) {
int x = dp[j + t], y = 1ll * Pw * dp[j + (1 << i) + t] % md;
dp[j + t] = (x + y) % md;
dp[j + (1 << i) + t] = (x - y + md) % md;
Pw = 1ll * Pw * w % md;
}
}
}
}
struct Poly {
vector<int> dp;
Poly(int x = 0) { dp.resize(x); }
inline auto size() { return dp.size(); }
inline int &operator[](int t) { return dp[t]; }
inline Poly Reverse() {
Poly res = *this;
reverse(res.dp.begin(), res.dp.end());
return res;
}
inline Poly operator-(Poly b) {
Poly res(max(dp.size(), b.size()));
for (int i = 0; i < dp.size(); i++) res[i] = (res[i] + dp[i] + md) % md;
for (int i = 0; i < b.size(); i++) res[i] = (res[i] - b[i] + md) % md;
return res;
}
inline Poly operator*(int x) {
Poly res = *this;
for (auto &it : res.dp) it = 1ll * it * x % md;
return res;
}
inline Poly operator*(Poly y) {
Poly x = *this, res(dp.size() + y.size() - 1);
int lim = 0;
while ((1 << lim) <= res.size()) lim++;
for (int i = 0; i < (1 << lim); i++) r[i] = (r[i >> 1] >> 1) + ((i & 1) << (lim - 1));
vector<int> a(1 << lim), b(1 << lim), c(1 << lim);
for (int i = 0; i < dp.size(); i++) a[i] = dp[i];
for (int i = 0; i < y.size(); i++) b[i] = y[i];
NTT(a, lim, G);
NTT(b, lim, G);
for (int i = 0; i < (1 << lim); i++) c[i] = 1ll * a[i] * b[i] % md;
NTT(c, lim, Gi);
int inv = pwr(1 << lim, md - 2);
for (int i = 0; i < res.size(); i++) res[i] = 1ll * c[i] * inv % md;
return res;
}
inline Poly Inv(int len) {
if (len == 1) {
Poly res(1);
res[0] = pwr(dp[0], md - 2);
return res;
}
Poly res = this->Inv((len + 1) / 2), tmp = *this;
tmp.dp.resize(len);
res = res * 2 - tmp * res * res;
res.dp.resize(len);
return res;
}
inline Poly operator/(Poly b) {
int N = dp.size(), M = b.size();
if (N < M)
return Poly();
Poly res = this->Reverse() * (b.Reverse().Inv(N - M + 1));
res.dp.resize(N - M + 1);
return res.Reverse();
}
inline Poly operator%(Poly b) {
Poly G = *this / b, res = *this - (G * b);
res.dp.resize(min(dp.size(), b.size()));
return res;
}
};
inline Poly Polypower(Poly x, int y, Poly MD) {
Poly res(1);
res[0] = 1;
while (y) {
if (y & 1)
res = res * x % MD;
x = x * x % MD;
y >>= 1;
}
return res;
}
int a[40005], k;
int main() {
scanf("%d%d", &k, &n);
Poly F(n + 1), X(2);
F[n] = 1;
X[1] = 1;
for (int i = n - 1; ~i; i--) scanf("%d", &F[i]), F[i] = (md - F[i]) % md;
for (int i = 0; i < n; i++) scanf("%d", &a[i]);
Poly tmp = Polypower(X, k, F);
int res = 0;
for (int i = 0; i < n; i++) res = (res + 1ll * a[i] * tmp[i]) % md;
printf("%d\n", (res + md) % md);
return 0;
}
