标签:理解 公交 数组 实现 amp 矩阵加速 mst 使用 tps
矩阵(Matrix)是一个由若干个数有序排列组成的集合。在这里我们只从 OI 的角度研究它。
由 \(m \times n\) 个数组成的 \(m\) 行 \(n\) 列的数表称为 \(m\) 行 \(n\) 列的矩阵,简称 \(m \times n\) 矩阵。
上面的描述指出矩阵本质上是个二维数组。可以把它想象为一个 \(m \times n\) 的长方形棋盘,棋盘的每一个方格中都放有一个数。在数学中,我们称这 \(m \times n\) 个数为矩阵的元素,而位于第 \(i\) 行第 \(j\) 列的元素可记做 \(a_{i,j}\)。
常用两个括号将矩阵的所有元素括起,表示一个矩阵。比如,一个 \(4 \times 5\) 的矩阵 \(A\) 可以这样表示:
\[A = \left[ \begin{array}{ccc} a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & a_{1,5}\a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & a_{2,5}\a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5}\a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5}
\end{array}\right]\]
对于两个矩阵 \(A,B\),如果它们的行数和列数均相等,就称 \(A\) 和 \(B\) 为同型矩阵;如果一个矩阵 \(A\) 的行数和列数都为 \(n\),就称 \(A\) 为 \(n\) 阶方阵。
这些定义都非常简单易懂,有的小学生都能理解。这么简单的东西,究竟有什么用呢?
如果没有运算,矩阵只是一个“路标”。因为有了运算,矩阵的力量无限。
这段话出自人教版必修 4 数学课本,原本被用于描述向量。我认为这段话用来描述矩阵也十分恰当。矩阵的运算赋予矩阵灵魂,让它在各个学科都发挥无限的力量。从信息学的角度来讲,这些运算都是非常便于实现的。
template <class Tp>
class Matrix{
private:
Tp x[MAXN][MAXN];
public:
Matrix(){
std::memset(x, 0, sizeof x);
}
Tp *operator[](const int& b){
return x[b];
}
};
矩阵的基本运算包括加法,只有同型矩阵间才能进行加法运算。运算的得到的新矩阵 \(C\) 中的每一个元素 \(c_{i,j}\) 恰是原矩阵 \(A,B\) 对应位置的两个元素的和,即 \(c_{i,j}=a_{i,j}+b_{i,j}\)。
\[A + B = \left[ \begin{array}{ccc} a_{1, 1} + b_{1, 1} & ... & ...\... & ... & ...\... & ... & a_{m, n} + b_{m, n} \end{array}\right]\]
矩阵的加法运算满足结合律,即 \[A + B + C = (A + B) + C = A + (B + C)\]
矩阵的加法运算也满足交换律,即 \[A + B = B + A\]
Matrix::Matrix operator+(const Matrix& b)const{
Matrix res;
for(int i = 0; i < MAXN; ++i)
for(int j = 0; j < MAXN; ++j)
res[i][j] = x[i][j] + b[i][j];
return res;
}
矩阵的基本运算包括减法。它的性质和加法的基本一致,不过不满足交换律,即 \(A - B \neq B - A\)。
\[A - B = \left[ \begin{array}{ccc} a_{1, 1} - b_{1, 1} & ... & ...\... & ... & ...\... & ... & a_{m, n} - b_{m, n} \end{array}\right]\]
Matrix::Matrix operator-(const Matrix& b)const{
Matrix res;
for(int i = 0; i < MAXN; ++i)
for(int j = 0; j < MAXN; ++j)
res[i][j] = x[i][j] - b[i][j];
return res;
}
和向量一样,矩阵可以乘上一个数,得到一个新的矩阵。设 \(C = x \cdot A\),那么对于 \(C\) 中的每一个元素,有 \(c_{i,j} = x \cdot a_{i,j}\)。
\[ x \cdot A = \left[ \begin{array}{ccc} x \cdot a_{1, 1} & ... & ...\... & ... & ...\... & ... & x \cdot a_{m, n} \end{array} \right] \]
Matrix scm(Matrix mat, const int& x){
for(int i = 0; i < MAXN; ++i)
for(int j = 0; j < MAXN; ++j)
mat[i][j] *= x;
return mat;
}
这是在 OI 中,矩阵最常用、最重要的运算。和加法、减法和数乘不一样,它的定义与我们的惯常的认知不符,却拥有更为强大的作用。
矩阵 \(A\) 和 \(B\) 能够相乘,当且仅当 \(A\) 的列数和 \(B\) 的行数相等。它们相乘的结果 \(C\) 的行数与 \(A\) 的行数相等,列数与 \(B\) 的列数相等。设 \(A, B, C\) 的行数分别为 \(m_a, m_b, m_c\),列数分别为 \(n_a, n_b, n_c\),这个关系可以抽象成:
\[ \left\{ \begin{aligned} n_a = m_b \m_c = m_a \n_c = n_b \end{aligned} \right. \]
\(C\) 是 \(A\) 和 \(B\) 的积,记做 \(C = AB\) 或 \(C = A \times B\)。对于 \(C\) 中的每一个元素,有:\(c_{i,j} = \sum_{k=1}^na_{i,k}\cdot b_{k,j}\)。
这个奇怪的定义常常令初学者疑惑。不过随着对矩阵了解的加深,你一定会感叹它的巧妙。
非常特别的,矩阵的乘法不满足交换律,即 \[A \times B \neq B \times A\]
不过,矩阵乘法满足结合律,即 \[A \times B \times C = (A \times B) \times C = A \times (B \times C)\]
Matrix::Matrix operator*(const Matrix& b)const{
Matrix res;
for(int i = 0; i < MAXN; ++i)
for(int j = 0; j < MAXN; ++j)
for(int k = 0; k < MAXN; ++k)
res[i][j] += x[i][k] * b[k][j];//是不是有点儿像 Floyd?
return res;
}
从代码看,矩阵乘法与那些线性运算的最大区别就是它有三层循环,时间复杂度为 \(O(n ^ 3)\)。
矩阵乘法非常值得研究,往下的内容几乎都是围绕它展开的。
符合我们的常识,矩阵的幂运算 \(A^k\) 的意义为 \(k\) 个 \(A\) 相乘,即
\[
A^k= \underbrace{ A \times A \times ... \times A \times A } \\qquad \, k 个 A
\]
快速幂是一种非常经典的算法,它运用倍增的思维加速了运算。
template <typename Tp>
Tp pow(Tp a, int b){
Tp res = E;//E 是什么?
for(; b; b >>= 1, a = a * a)
if(b & 1)
res = res * a;
return res;
}
它成立的根本是乘法的结合律,正如 Dijkstra 算法成立的根本是边权非负。矩阵乘法也满足结合律,所以矩阵的幂也可以用快速幂求出。
不知道你是否思考过一个问题:\(A^0\) 是什么?我们都知道,任何非零数的 \(0\) 次方都等于 \(1\),任何数乘上 \(1\) 都等于它本身,矩阵中是否存在这样的“1”呢?
这就是单位矩阵,常用 \(E_n\) 或 \(I_n\) 表示。任何矩阵的 \(0\) 次方都是相应的单位矩阵,任何矩阵乘上相应的单位矩阵都等于它本身。
\[E_3 = \left[ \begin{array}{ccc} 1 & 0 & 0\0 & 1 & 0\0 & 0 & 1 \end{array} \right]\]
\[ E_4 = \left[ \begin{array}{ccc} 1 & 0 & 0 & 0\0 & 1 & 0 & 0\0 & 0 & 1 & 0\0 & 0 & 0 & 1 \end{array} \right] \]
\(E_n\) 是一个左上角至右下角对角线(称为主对角线)上的 \(n\) 个数为 \(1\),其他数为 \(0\) 的 \(n\) 阶方阵。
为了避免单位矩阵的构造,常常写这样的快速幂:
Matrix pow(Matrix a, int b){
Matrix res = a; --b;//a ^ b = a * a ^ (b - 1)
while(b){
if(b & 1)
res = res * a;
b >>= 1, a = a * a;
}
return res;
}
一些动态规划的转移方程可以结合矩阵优化。
乘法与方案数有着不解之缘。如果你没有学过小学奥数,建议先自行百度了解一下乘法原理(基本的乘法原理对于高中生来说应该很好理解)。
给一张 \(n\) 个点 \(m\) 条边的无向图,初始时你在 \(1\) 号点。每秒你都可以停在原地、去往一个相邻的点或结束游戏。求在 \(t\) 秒的时间内你有多少种行为方案?
\(1 \leq n \leq 30, 0 < m < 100, 1 \leq t \leq 10^6\)
不难想出朴素 DP:
\(\quad\) \(\bullet\) 设 \(dp_{i, j}\) 表示 \(i\) 秒时到达 \(j\) 的方案数。
\(\quad\) \(\bullet\) \(dp_{i, j} = \sum_{k = i ∨ k \text{ 能够到达 } j} dp_{i - 1, k}\)
答案就是 \(\sum_{i = 1} ^ t \sum_{j = 1} ^ n dp_{i, j}\),时间复杂度为 \(O(m \cdot t)\),于本题是可接受的。
有没有更好的做法呢?我们发现,这样的转移方程也是完全正确的:
\(\quad\) \(\bullet\) 设 \(dp_{t, i, j}\) 表示 \(t\) 秒时从 \(i\) 走到 \(j\) 的方案数。
\(\quad\) \(\bullet\) 根据计数原理中的乘法原理,\(dp_{t, i, j} = \sum_{k = 1} ^ n dp_{p, i, k} \cdot dp_{t - p, k, j}\)
其中 \(p\) 为任意满足 \(1 \leq p \leq t\) 的数。用文字语言描述,这个转移方程的含义是:花费 \(t\) 秒从 \(i\) 走到 \(j\) 的方案数,等于花费 \(p\) 秒从 \(i\) 走到 \(k\) 的方案数与花费 \(t - p\) 秒从 \(k\) 走到 \(j\) 的方案数之积。
那么有
\[dp_{t, i, j} = \sum_{k = 1} ^ n dp_{t - 1, i, k} \cdot dp_{1, k, j}\]
单次转移的时间复杂度是 \(O(n ^ 3)\)。
设矩阵 \(A_t\) 代表 \(dp_t\),根据矩阵乘法的定义 \(c_{i,j} = \sum_{k=1}^na_{i,k}\cdot b_{k,j}\) 可得 \(A_t = A_{t - p} \times A_p\)。
\[\begin{aligned}
A_t & = A_{t - 1} \times A_1\& = A_{t - 2} \times A_1 \times A_1\& = A_{t - 3} \times A_1 \times A_1 \times A_1\& = ......\& = {A_1} ^ t
\end{aligned}\]
\(A_1\) 是题目给定的,只要一个乘法快速幂就可以快速求出 \(A_t\) 了。由于矩阵乘法本身复杂度为 \(O(n ^ 3)\),总时间复杂度为 \(O(n ^ 3 \log t)\)。
#include <cstdio>
#include <cstring>
const int MAXN = 3e1 + 19, MOD = 2017;
struct Matrix{
int e[MAXN][MAXN];
Matrix(){
std::memset(e, 0, sizeof e);
}
Matrix operator*(const Matrix& b)const{
Matrix res;
for(int i = 0; i < MAXN; ++i)
for(int j = 0; j < MAXN; ++j)
for(int k = 0; k < MAXN; ++k)
res.e[i][j] = (res.e[i][j] + e[i][k] * b.e[k][j]) % MOD;
return res;
}
}ans, trans;
Matrix pow(Matrix a, int b){
Matrix res = a; --b;
while(b){
if(b & 1)
res = res * a;
b >>= 1, a = a * a;
}
return res;
}
int n, m, t, s;
int main(){
std::scanf("%d%d", &n, &m);
for(int i = 0; i <= n; ++i)
trans.e[i][0] = 1,//把自爆当作 0 号点,自爆就相当于去往 0 号点
trans.e[i][i] = 1;
while(m--){
int u, v;
std::scanf("%d%d", &u, &v);
trans.e[u][v] = 1, trans.e[v][u] = 1;/*
如果 u 和 v 之间有路,dp(u, v) = dp(v, u) = 1
*/
}
std::scanf("%d", &t);
ans = pow(trans, t);//ans = A1 ^ t
for(int i = 0; i <= n; ++i)
s = (s + ans.e[1][i]) % MOD;
std::printf("%d\n", s);
return 0;
}
\(8\) 个公交站台 \(A, B, C, D, E, F, G, H\) 按顺序成环状分布,每次乘车可以前往相邻站台。求从 \(A\) 出发,乘 \(n\) 次车到达 \(E\) 有多少种方案?(到达 \(E\) 后就不再继续乘车)
\(4 \leq n \leq 10 ^ 7\)
如果你学懂了上面一题,这道题应该能秒切。(由于数据范围小,朴素 dp 也能 AC)
int n;
Matrix trans;
int main(){
trans[1][2] = 1, trans[2][3] = 1, trans[3][4] = 1;
trans[2][1] = 1, trans[3][2] = 1, trans[4][3] = 1;
trans[1][5] = 1, trans[5][6] = 1, trans[6][7] = 1;
trans[5][1] = 1, trans[6][5] = 1, trans[7][6] = 1;
trans[7][0] = 1, trans[4][0] = 1;/*
1 - A, 0 - E
*/
std::scanf("%d", &n);
std::printf("%d\n", pow(trans, n)[1][0]);
return 0;
}
给一张包含 \(n\) 个点的有向图,每条边的边权均为 \(1\) 至 \(9\) 的正整数。求从 \(0\) 号点到达 \(n - 1\) 号点,且正好花费 \(T\) 时间的方案数。
\(2 \leq n \leq 10, 1 \leq T \leq 10 ^ 9\)
如果使每条边的边权都为 \(1\),这道题和可乐就一样了。怎么处理不同的边权呢?
观察到边权较小,可以暴力拆边。把一条长度为 \(9\) 的边强行拆成 \(9\) 条长度为 \(1\) 的边即可。
#include <cstdio>
#include <cstring>
const int MOD = 2009;
class Matrix{
private:
int x[100][100];
public:
Matrix(){
std::memset(x, 0, sizeof x);
}
int *operator[](const int& b){
return x[b];
}
Matrix operator*(const Matrix& b)const{
Matrix res;
for(int i = 0; i < 100; ++i)
for(int j = 0; j < 100; ++j)
for(int k = 0; k < 100; ++k)
res.x[i][j] = (res.x[i][j] + x[i][k] * b.x[k][j]) % MOD;
return res;
}
}mymatrix;
Matrix pow(Matrix a, int b){
Matrix res = a; --b;
while(b){
if(b & 1)
res = res * a;
b >>= 1;
a = a * a;
}
return res;
}
int n, t;
char s[19];
int main(){
std::scanf("%d%d", &n, &t);
for(int i = 1; i <= n; ++i){
std::scanf("%s", s + 1);
for(int j = 1; j <= n; ++j){
int len = s[j] - '0';
if(len == 0)
continue;
for(int k = 0; k < len - 1; ++k)
mymatrix[(i - 1) * 10 + k][(i - 1) * 10 + k + 1] = 1;
if(len > 0)
mymatrix[(i - 1) * 10 + len - 1][(j - 1) * 10] = 1;
}
}
std::printf("%d\n", pow(mymatrix, t)[0][(n - 1) * 10]);
return 0;
}
给一张有 \(T\) 条边的无向连通图 \(G\),和 \(G\) 上的两点 \(s, e\)。求从 \(s\) 到 \(e\) 经过 \(k\) 条边的最短路。
\(1 \leq n \leq 10 ^ 6, 2 \leq T \leq 100\),所有点的编号不大于 \(1000\),所有边的边权不大于 \(1000\)。
设 \(dist_{t, i, j}\) 代表从 \(i\) 到 \(j\),经过 \(t\) 条边的最短路。显然,
\[dist_{t, i, j} = \min_{k = 1} ^ n (dist_{t - 1, i, k} + dist_{1, k, j})\]
定义矩阵运算 \(\bigodot\):
若
\[C = A \bigodot B\]
那么对于 \(C\) 中的每一个元素 \(c_{i, j}\) 都有
\[c_{i, j} = \min_{k = 1} ^ n (a_{i, k} + b_{k, j})\]
可以证明,\(\bigodot\) 运算满足结合律,所以它也可以倍增快速幂计算。
#include <algorithm>
#include <cstdio>
#include <cstring>
const int MAXN = 1e2 + 19;
inline int min(const int& a, const int& b){
return a < b ? a : b;
}
class Matrix{
private:
int x[MAXN][MAXN];
public:
Matrix(){
std::memset(x, 0, sizeof x);
}
int *operator[](const int& b){
return x[b];
}
Matrix operator*(const Matrix& b)const{//这个乘法不是真正的矩阵乘法,是我们自定义的 min 运算
Matrix res;
std::memset(res.x, 0x3f, sizeof res.x);
for(int i = 0; i < MAXN; ++i)
for(int j = 0; j < MAXN; ++j)
for(int k = 0; k < MAXN; ++k)
res[i][j] = min(res[i][j], x[i][k] + b.x[k][j]);
return res;
}
void set(void){
std::memset(x, 0x3f, sizeof x);
}
}trans;
inline Matrix pow(Matrix a, int b){//由于 min 运算也满足结合律,它也可以用快速幂计算
Matrix res = a; --b;
while(b){
if(b & 1)
res = res * a;
b >>= 1, a = a * a;
}
return res;
}
int n, t, s, e;
int node[MAXN << 1];
int w[MAXN], u[MAXN], v[MAXN];
int main(){
std::scanf("%d%d%d%d", &n, &t, &s, &e);
for(int cnt = 0, i = 1; i <= t; ++i){
std::scanf("%d%d%d", w + i, u + i, v + i);
node[++cnt] = u[i], node[++cnt] = v[i];
}
std::sort(node + 1, node + t + t + 1);
int *end = std::unique(node + 1, node + t + t + 1);//离散化。
trans.set();
for(int i = 1; i <= t; ++i){
int a = std::lower_bound(node + 1, end, u[i]) - node - 1,
b = std::lower_bound(node + 1, end, v[i]) - node - 1;
trans[a][b] = w[i],
trans[b][a] = w[i];
}
std::printf("%d\n", pow(trans, n)[std::lower_bound(node + 1, end, s) - node - 1][std::lower_bound(node + 1, end, e) - node - 1]);
return 0;
}
举一反三,\(\max\) 在矩阵中也具有满足结合律的性质。这种拓展运算的题不多见,感兴趣的话可以再做一下 P3502 [POI2010]CHO-Hamsters。
最著名的线性递推就是斐波那契数列 \(f(x) = f(x - 1) + f(x - 2)\)。矩阵可以把递推加速到 \(O(\log n)\) 的复杂度内。
对于一个数列 \(\{a_n\}\),如果它的每一项都可由前面的某 \(n\) 项推得,那这个递推就是 \(n\) 阶的。学过高中数学的都知道,一阶线性递推可以直接求出通项公式,我在这里简要提一下。
对一阶线性递推式 \(a_n=p \cdot a_{n-1} + q\) 求 \(a_n\) 的通项公式:
\[
\frac {a_n} {p ^ n} = \frac {a_{n - 1}} {p^{n - 1}} + \frac q {p ^ n}
\]
令 \(b_n = \frac {a_n} {p^n}\),那么
\[
b_n = b_{n - 1} + \frac q{p ^ n}
\]
\[
\begin{aligned}\therefore
b_n - b_1 & = (b_n - b_{n - 1}) + (b_{n - 1} - b_{n - 2}) + ... + (b_2 - b_1)\& = \sum_{i = 2} ^ n (b_i - b_{i - 1})\& = \sum_{i = 2} ^ n \frac q{p ^ i}\& = q \cdot \sum_{i = 2} ^ n \frac 1 {p ^ i}\& = \frac {q (p ^ {n - 1} - 1)}{(p - 1) p ^ n}\\therefore b_n & = b_1 + \frac {q (p ^ {n - 1} - 1)} {(p - 1) p ^ n}\& = \frac {a_1} p + \frac {q (p ^ {n - 1} - 1)} {(p - 1) p ^ n}\\therefore a_n & = (a_1 + \frac q {p - 1}) p ^ {n - 1} - \frac q {p - 1}
\end{aligned}
\]
这只是众多递推公式求法之一,求单项的时间复杂度为 \(O(1)\)。由于它过于简单,并且没有优化空间,这里不做研究。
这种做法涉及分数运算,会损失精度。在常数项为 \(0\) 的特殊情况下,可以参照二阶常系数齐次线性递推的做法。
可以用待定系数法构造等比数列来求出通项公式。
对于数列 \[a_n = p \cdot a_{n - 1} + q \cdot a_{n - 2}\]
设有 \(r, s\) 满足
\[a_n - r \cdot a_{n - 1} = s (a_{n - 1} - r \cdot a_{n - 2})\]
那么
\[ \begin{aligned} a_n - r \cdot a_{n - 1} & = s ^ {n - 2} \cdot (a_2 - r \cdot a_1)\a_n & = r \cdot a_{n - 1} + s ^ {n - 2} \cdot (a_2 - r \cdot a_1) \end{aligned} \]
由于 \(r, s, a_1, a_2\) 都是已知的,这样我们就把二阶线性递推转化为一阶的了。
斐波那契数列的通项公式
\[ f_n = \frac 1 {\sqrt 5} \left[ \left( \frac {1 + \sqrt 5} 2 \right) ^ n - \left( \frac {1 - \sqrt 5} 2 \right) ^ n \right] \]
二阶线性递推的公式出现了幂,所以时间复杂度是 \(O(\log n)\),除此之外,由于涉及分数、无理数运算,会损失相当可观的精度。
为了避免浮点数运算,我们使用矩阵,将递推关系转换为矩阵的乘法。构造矩阵 \(A, B\):
\[A = \left[ \begin{array}{ccc} a_1 & a_2\0 & 0 \end{array} \right], B = \left[ \begin{array}{ccc} 0 & q\1 & p \end{array} \right]\]
根据矩阵乘法的意义
\[A \times B = \left[ \begin{array}{ccc} p \cdot a_2 + q \cdot a_1 & a_2\0 & 0 \end{array} \right] = \left[ \begin{array}{ccc} a_3 & a_2\0 & 0 \end{array} \right]\]
同样地
\[A \times B ^ 2 = \left[ \begin{array}{ccc} a_4 & a_3\0 & 0 \end{array} \right], A \times B ^ 3 = \left[ \begin{array}{ccc} a_5 & a_4\0 & 0 \end{array} \right], A \times B ^ n = \left[ \begin{array}{ccc} a_{n + 2} & a_{n + 1}\0 & 0 \end{array} \right]\]
要求 \(a_n\),只需求 \(A \times B ^ {n - 2}\) 即可。借助矩阵快速幂,可以保证 \(O(\log n)\) 的时间复杂度。
Matrix A, B;
int n, a1, a2;
int main(){
std::scanf("%d%d%d", &n, &a1, &a2);
A[0][0] = a1, A[0][1] = a2;
B[0][1] = q, B[1][0] = 1, B[1][1] = p;
std::printf("%d\n", (A * pow(B, n - 2))[0][0]);
return 0;
}
相关的模板题:
P1962 斐波那契数列
P1939 【模板】矩阵加速(数列)
事实上,更高阶的常系数齐次线性递推也可以求出通项公式。然而随着这些公式愈来愈复杂,矩阵的优越性也愈加体现出来了。
\[g_n = \sum_{i = 1} ^ k f_i \cdot g_{n - i}\]
高阶常系数齐次线性递推的矩阵该如何构造呢?我们发现,对于 \(C = A \times B\),\(c_{i, j} = \sum_{k = 1}^n a_{i, k} \cdot b_{k, j}\)。我们使 \(c_{i, n}\) 表示数列的第 \(n\) 项。而 \(g_n = \sum_{i = 1}^k g_i \cdot b_{i, n}\) 正好与递推公式相吻合。那么原矩阵 \(A\) 就可以构造成一个第一行依次为 \(a_1, a_2, a_3, ..., a_k\) 的矩阵,转移矩阵 \(B\) 构造成每一个元素都是相应的 \(f_i\) 的矩阵。
为了解决 \(k\) 阶的递推,矩阵也必须是 \(k\) 阶的,时间复杂度为 \(O(k ^ 3 \cdot \log n)\)。这样的算法不足以通过洛谷上的模板题,不过也向我们展示了矩阵的强大威力。
标签:理解 公交 数组 实现 amp 矩阵加速 mst 使用 tps
原文地址:https://www.cnblogs.com/natsuka/p/12449834.html