浅谈矩阵加速
OneLastNOIP
·
·
个人记录
前置知识：矩阵乘法，快速幂
很喜欢我们科学老师的一句话:
有时候证明可以写一黑板，但可以想。
这篇文章某些深奥的内容可能要感性理解。
有些时候，我们经常遇到类似函数的求值：
\large f(n)=\sum_{i=1}^{3}A_in^{X_i}f(n-i)
这些类型的递推被称为线性递推。一个朴素的算法是直接递推。但世界上有个好东西，叫矩阵。
入门
我们以斐波那契数列为例。
构造 A_0=\begin{bmatrix}
1 & 1
\end{bmatrix}，即 f(2) 和 f(1)。
再是 B=\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}，设 A_i = A_0\times B^i，可以发现 A_i=\begin{bmatrix}
f(i + 2) & f(i + 1)
\end{bmatrix}。
来观察一下发生的过程。
设 A^\prime = A_i，C=A^\prime \times B，则 C_{1,1} = B_{1, 1}\times A^\prime_{1,1}+B_{2,1}\times A^\prime_{2,1}=f(i+1)+f(i+2)=f(i+3),C_{1,2}=B_{1,2}\times A^\prime_{1,2}=f(i+2)。
因此，f(i)=A_{i-2_{1,1}}。
那怎么计算 A_i 呢？使用快速幂计算 B^i 即可。
现在，省选与你(可能)有缘。
构造矩阵的技巧
假设我们用一行存储系数矩阵(A_i)，则若第 i 项的求和方法中包含 X\times A_{j}，则 B_{j,i}=X，在相乘时 A_j 与 X 相乘累加到答案的第 i 项上，这个过程需要对矩阵乘法有较直观的理解，建议深入思考一下。
Wikioi # 10375 矩阵加速入门 6
非常经典。如果只有 p\times f(n-1)+q\times f(n-2)（入门 3），应用技巧，A_0=\begin{bmatrix}
f(0),f(1)
\end{bmatrix}，A_{1,2} 要保存上一次的 A_{1,1}，即 A^\prime_{1,2}=1\times A_{1,1}，所以 B_{1,2}=1。再是 A^\prime_{1,1}=p\times A_{1,1}+q\times A_{1,2}，所以 B_{1,1}=p，B_{2,1}=q。
现在处理那三个系数。
给出两个矩阵。
```cpp
int a[1][5] = {1, 1, 1, 3, 9};// f(n-1) f(n) 1 n n^2
int b[5][5] =
{{0, b, 0, 0, 0},
{1, a, 0, 0, 0},
{0, 1, 1, 1, 1},
{0, 1, 0, 1, 2},
{0, 1, 0, 0, 1}
};
```
矩阵的模板下面代码里有。
## 练习
### LG P1962 斐波那契数列
思路见上。
```cpp
#include<bits/stdc++.h>
#define int long long
int mod = 1000000007 ;
using namespace std;
struct Matrix{
private:
int r,c;
vector<vector<int> >mat;
public:
void init(int x, int y, void *p) {
r = x, c = y;
mat.clear();
int (*q)[y] = (int (*)[y])(p);
for(int i = 0;i < x;i ++) {
vector<int>vec;
for(int j = 0;j < y;j ++)
vec.push_back(q[i][j]);
mat.push_back(vec);
}
}
void init(int x, int y) {
r = x, c = y;
mat.clear();
for(int i = 0;i < x;i ++) {
vector<int>vec;
for(int j = 0;j < y;j ++)
vec.push_back(0);
mat.push_back(vec);
}
}
Matrix operator*(Matrix &b) {
Matrix ans;
ans.init(r, b.c);
if(c != b.r) {
cout << "MulError: C != R" << endl;
exit(0);
}
for(int i = 0;i < r;i ++)
for(int j = 0;j < b.c;j ++) {
int sum = 0;
for(int k = 0;k < c;k ++)
sum += mat[i][k] * b.mat[k][j], sum %= mod;
ans.mat[i][j] = sum;
}
return ans;
}
void put(int x, int y) {
cout << mat[x][y];
}
};
int x[1][2] = {1, 1};
int b[2][2] =
{{0, 1},
{1, 1},
};
Matrix qpow(int a) {
Matrix ans, base;
ans.init(1, 2, x);
base.init(2, 2, b);
while(a) {
if(a % 2)
ans = ans * base;
base = base * base;
a >>= 1;
}
return ans;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
int n;
cin >> n;
if(n == 0)
cout << 1;
if(n <= 2)
cout << x[0][n - 1];
else
qpow(n - 2).put(0, 1);
cout << endl;
return 0;
}
```
### LG P4838 P哥破解密码
~~可惜啊，以前是紫的。~~
设 $f(i)$ 表示答案，推一下发现 $f(i)=f(i-1)+f(i-2)+f(i-3)$。算一下 $f(1),f(2),f(3)$，然后就是模板级别的应用。
### LG P4910 帕秋莉的手环
每次，我们可以在手环的末尾放入一个金属性、木属性的珠子。可以发现，如果末尾是木属性，下一个珠子就不会是木属性，因此我们把金属性和木属性绑定，每次只能添加一个金属性或“金属性+木属性”，可以证明不影响答案。现在就变成了斐波那契数列，只不过 $f(2)=3$ （金木，金金，木金）。
### AT ABC293E Geometric Progression
等比数列求和。一眼定真，$M$ 可能不为质数，$f(i)=A\times f(i-1)+1$，因为每次增加一个 $A^i$，相当于原等比数列每一项乘 $A$，再在头部加上 $1$。