「知识总结」特征多项式与线性递推【矩阵乘法】

特征多项式

定义

定义 n×nn \times n 的矩阵 AA 的特征多项式为:

p(λ)=det(λIA)=λn+g1λn1+g2λn2++gn\begin{aligned} p(\lambda) &= \det(\lambda I - A) \\ &= \lambda^n + g_1 \lambda^{n - 1} + g_2 \lambda^{n - 2} + \cdots + g_n \end{aligned}

其中 IIn×nn \times n 的单位矩阵,λ\lambda 可以为实数、向量或矩阵。

性质

根据 Cayley-Hamilton\text{Cayley-Hamilton} 定理,AA 满足:

p(A)=Op(A) = O

其中 OO 为零矩阵。则 k\forall k

Ak=Ak mod p(A)A^k = A^k \text{ mod } p(A)

上一步本质上将 AkA^k 表示为以 AA 为变量、次数小于 nn 的多项式。

示例

举个例子,考虑下列矩阵:

A=[1234]A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}

其特征多项式为:

p(λ)=λ123λ4=λ25λ2p(\lambda) = \begin{vmatrix} \lambda - 1 & -2 \\ -3 & \lambda - 4 \end{vmatrix} = \lambda^2 - 5 \lambda - 2

可以验证 Cayley-Hamilton\text{Cayley-Hamilton} 定理:

A25A2I=O Ak=Ak mod (A25A2I)\begin{aligned} & A^2 - 5A - 2I = O \\ \Rightarrow\ & A^k = A^k \text{ mod } (A^2 - 5A - 2I) \end{aligned}

求法

n+1n + 1 个值分别作为 λ\lambda 求行列式,最后插值得到特征多项式的系数。

由于高斯消元求行列式是 O(n3)O(n^3) 的,所以总时间复杂度为 O(n4)O(n^4)

另外,存在 O(n3)O(n^3) 的做法,可以参考 特征多项式与常系数线性齐次递推

例题 11

题目链接【BZOJ 4162】Shlw loves matrix II

题目大意:求 k×kk \times k 的矩阵 AAnn 次幂。保证 n50,n \leqslant 50, k210000k \leqslant 2^{10000}

思路分析:总时间复杂度为 O(k4+k2logn)O(k^4 + k^2 \log n)

  1. 首先用上述方法求出特征多项式 p(λ)p(\lambda)。复杂度为 O(n4)O(n^4)
  2. λn mod p(λ)\lambda^n \text{ mod } p(\lambda),暴力快速幂 ++ 取模。复杂度为 O(k2logn)O(k^2 \log n)
  3. AA 代入 λn mod p(λ)\lambda^n \text{ mod } p(\lambda),把矩阵代入多项式。复杂度为 O(k4)O(k^4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include <bits/stdc++.h>
using namespace std;

const int P = 1000000007;
int n, a[55][55], b[55][55], c[55][55], e[55][55];
int d[110], p[55], tmp[110], ans[55][55];
char s[10010];

int qp(int x, int y) {
int z = 1;
for (; y; y >>= 1, x = 1LL * x * x % P) {
if (y & 1) z = 1LL * z * x % P;
}
return z;
}

int det(int x) {
memcpy(b, a, sizeof(b));
for (int i = 0; i < n; i++) {
b[i][i] = (b[i][i] - x) % P;
}
int ans = 1;
for (int i = 0; i < n; i++) {
for (int j = i + 1; j < n; j++) {
if (b[j][i]) { ans = -ans, swap(b[i], b[j]); break; }
}
ans = 1LL * ans * b[i][i] % P;
for (int j = i + 1; j < n; j++) {
int r = 1LL * b[j][i] * qp(b[i][i], P - 2) % P;
for (int k = 0; k < n; k++) {
b[j][k] = (b[j][k] - 1LL * r * b[i][k]) % P;
}
}
}
return ans;
}

void gauss() {
for (int i = 0; i <= n; i++) {
for (int j = i + 1; j <= n; j++) {
if (e[j][i]) { swap(e[i], e[j]); break; }
}
for (int j = 0; j <= n; j++) if (i ^ j) {
int r = 1LL * e[j][i] * qp(e[i][i], P - 2) % P;
for (int k = 0; k <= n + 1; k++) {
e[j][k] = (e[j][k] - 1LL * r * e[i][k]) % P;
}
}
}
for (int i = 0; i <= n; i++) {
p[i] = 1LL * e[i][n + 1] * qp(e[i][i], P - 2) % P;
}
}

void mul() {
memset(c, 0, sizeof(c));
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
c[i][k] = (c[i][k] + 1LL * a[i][j] * b[j][k]) % P;
}
}
}
}

int main() {
scanf("%s %d", s, &n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
scanf("%d", &a[i][j]);
}
}
for (int i = 0; i <= n; i++) {
for (int j = 0; j <= n; j++) {
e[i][j] = qp(i, j);
}
e[i][n + 1] = det(i);
}
gauss(), d[0] = 1;
for (int i = 0; s[i]; i++) {
memcpy(tmp, d, sizeof(d)), memset(d, 0, sizeof(d));
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
d[j + k] = (d[j + k] + 1LL * tmp[j] * tmp[k]) % P;
}
}
if (s[i] == '1') {
for (int j = 2 * n - 1; j; j--) {
d[j] = d[j - 1];
}
d[0] = 0;
}
for (int j = 2 * n - 1; j >= n; j--) {
for (int k = n; ~k; k--) {
d[j - k] = (d[j - k] - 1LL * d[j] * p[n - k]) % P;
}
}
}
memset(b, 0, sizeof(b));
for (int i = 0; i < n; i++) {
b[i][i] = 1;
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
ans[j][k] = (ans[j][k] + 1LL * d[i] * b[j][k]) % P;
}
}
mul(), memcpy(b, c, sizeof(c));
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
printf("%d%c", (ans[i][j] + P) % P, " \n"[j == n - 1]);
}
}
return 0;
}

常系数线性齐次递推

定义

给定数列 ff 的前 kkf0,f1,,fk1f_0, f_1, \cdots, f_{k - 1} 和系数 c0,c1,,ck1c_0, c_1, \cdots, c_{k - 1}

ik, fi=j=1kfijcj\forall i \geqslant k,\ f_i = \sum_{j = 1}^k f_{i - j} c_j。需要快速求 fnf_n 的值。

方法一

构造出转移矩阵,然后矩阵快速幂加速:

(fk1fk2fk3f0)[c1c2ck1ck100001000010]=(fkfk1fk2f1)\begin{pmatrix} f_{k - 1} \\ f_{k - 2} \\ f_{k - 3} \\ \vdots \\ f_0 \end{pmatrix} \begin{bmatrix} c_1 & c_2 & \cdots & c_{k - 1} & c_k \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix} = \begin{pmatrix} f_k \\ f_{k - 1} \\ f_{k - 2} \\ \vdots \\ f_1 \end{pmatrix}

时间复杂度为 O(k3logn)O(k^3 \log n)

方法二

将【方法一】中的转移矩阵设为 AA,其特征多项式为:

p(λ)=det(λIA)=det([λc1c2ck1ck1λ000100001λ])=λki=1kciλki\begin{aligned} p(\lambda) &= \det(\lambda I - A) \\ &= \det \left ( \begin{bmatrix} \lambda - c_1 & -c_2 & \cdots & -c_{k - 1} & -c_k \\ -1 & \lambda & \cdots & 0 & 0 \\ 0 & -1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & -1 & \lambda \end{bmatrix} \right ) \\ &= \lambda^k - \sum_{i = 1}^k c_i \lambda^{k - i} \end{aligned}

f=(fk1fk20),\vec{f} = \begin{pmatrix} f_{k - 1} \\ f_{k - 2} \\ \cdots \\ 0 \end{pmatrix}, g=λn mod p(λ)g = \lambda^n \text{ mod } p(\lambda)[a]k[\vec{a}]_k 表示列向量 a\vec{a} 的第 kk 项,则:

fn=[Anf]k=[i=0k1giAif]k=i=0k1gi[Aif]k=i=0k1gifi\begin{aligned} f_n &= \big[ A^n \vec{f} \big]_k \\ &= \bigg[ \sum_{i = 0}^{k - 1} g_i A^i \vec{f} \bigg]_k \\ &= \sum_{i = 0}^{k - 1} g_i \big[ A^i \vec{f} \big]_k \\ &= \sum_{i = 0}^{k - 1} g_i f_i \end{aligned}

所以只用求出 λn mod p(λ)\lambda^n \text{ mod } p(\lambda),暴力取模即可。

时间复杂度为 O(k2logn)O(k^2 \log n)

方法三

使用 NTT\text{NTT} 加速【方法二】中的多项式乘法和取模。

时间复杂度为 O(klogklogn)O(k \log k \log n)

例题 22

题目链接【BZOJ 4161】Shlw loves matrix I

题目大意:与【定义】中描述的问题相同。保证 n109,n \leqslant 10^9, k2000k \leqslant 2000

思路分析:使用【方法二】即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include <bits/stdc++.h>
using namespace std;

const int maxk = 4010, P = 1000000007;
int n, ans, k, p[maxk], a[maxk], b[maxk], c[maxk], d[maxk], t[maxk];

void fft(int *a, int *b, int *c) {
fill(t, t + k + k, 0);
for (int i = 0; i < k; i++) {
for (int j = 0; j < k; j++) {
t[i + j] = (t[i + j] + 1LL * a[i] * b[j]) % P;
}
}
for (int i = 2 * k - 2; i >= k; i--) {
for (int j = 1; j <= k; j++) {
t[i - j] = (t[i - j] - 1LL * t[i] * p[k - j]) % P;
}
}
copy(t, t + k, c);
}

int main() {
scanf("%d %d", &n, &k);
for (int i = 1; i <= k; i++) {
scanf("%d", &c[i]), p[k - i] = -c[i];
}
p[k] = b[1] = d[0] = 1;
for (; n; n >>= 1, fft(b, b, b)) {
if (n & 1) fft(d, b, d);
}
for (int i = 0, a; i < k; i++) {
scanf("%d", &a), ans = (ans + 1LL * a * d[i]) % P;
}
printf("%d\n", (ans + P) % P);
return 0;
}