浅谈矩阵及其基本运算

· · 算法·理论

这篇文章主要用于讲解矩阵及其应用。

由于作者太菜，如有失误欢迎指出。

不一定更好的阅读体验

本文主要讲解以下几个内容：

0. 高斯消元

模板题见 P3389。

高斯消元法是一个用来求解线性方程组的算法。

可以解出形如

\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}

的 n 元一次方程组。

高斯消元主要是运用两种消元方法：加减消元法和代入消元法，和初中生解一次方程组一样。

高斯消元的方法是，使用主元法，先进行加减消元，然后我们可以先求得 x_1，然后可以逐层往回代，依次可以得到 x_2,x_3\cdots x_n。

举个例子，方程

\begin{cases} x + 3y + 4z = 5 \\ x + 4y + 7z = 3 \\ 9x + 3y + 2z = 2 \end{cases}

先写进矩阵：

1 & 3 & 4 & | & 5 \\ 1 & 4 & 7 & | & 3 \\ 9 & 3 & 2 & | & 2 \end{bmatrix}

然后通过主元法，系数化为 1，上下加减，得到：

0.25 & 0.75 & 1 & | & 1.25 \\ 0.14 & 0.57 & 1 & | & 0.43 \\ 4.5 & 1.5 & 1 & | & 1 \end{bmatrix}

加减消元：

\begin{cases} 0.11x + 0.18y = 0.82 \\ -4.36x -0.57y = -0.57\end{cases}

使用相同的方法：

可得 x=3.98\times -\dfrac{1}{4.09}=0.97

但是怎么回代呢？

我们每次操作的时候留下一个当前的方程，矩阵会变成这样：

-4.09 & 1 & 1 & | & 3.98 \\ 4.69 & 1 & 1 & | & 0.62 \\ 4.5 & 1.5 & 1 & | & 1 \end{bmatrix}

回代，解得

\begin{cases} x = -0.97 \\ y = 5.18 \\ z = -2.39 \end{cases}

如果一个方程没有了系数不为 0 的项，则方程无解或者有无数解，输出 No Solution。

复杂度 \mathcal{O}(n^3)，可以过。

参考代码：

#include<bits/stdc++.h>
using namespace std;
const int N = 110;
const double eps = 1e-6;
int n; double a[N][N];
bool gauss() {
    int r = 0;
    for(int c = 0; c < n; c++) {
        int t = r;
        for(int i = r; i < n; i++)
            if(fabs(a[i][c]) > fabs(a[t][c]))
                t = i;
        if(fabs(a[t][c]) < eps) continue;
        for(int i = c; i < n + 1; i++) swap(a[t][i], a[r][i]);
        for(int i = n; i >= c; i--) a[r][i] /= a[r][c];
        for(int i = r + 1; i < n; i++)
            if(fabs(a[i][c]) > eps)
                for(int j = n; j >= c; j--)
                    a[i][j] -= a[r][j] * a[i][c];
        r++;
    }
    if(r < n) return 1;
    for(int i = n - 1; i >= 0; i--) for(int j = i + 1; j < n; j++) a[i][n] -= a[j][n] * a[i][j];
    return 0;
}
int main() {
    cin >> n;
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n + 1; j++)
            cin >> a[i][j];
    bool t = gauss();
    if(t == 0) for(int i = 0; i < n; i++) printf("%.2lf\n", a[i][n]);
    else puts("No Solution");
    return 0;
}

类似的还有 P2455。

高斯消元在后面的逆矩阵和矩阵行列式中有用到，虽然不属于矩阵范畴，但还是提一下。

1. 矩阵的定义

进入正题：

在了解矩阵的运算之前，我们需要了解矩阵是啥。

一个 n \times m 的矩阵是一个由 n 行 m 列元素排列成的矩形阵列。即形如

A = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,m} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n,1} & a_{n,2} & \cdots & a_{n,m} \end{bmatrix} \text{.}

矩阵其实就是一个二维数组。

把数按照方形排列起来，就是一个是矩阵了。

与二维数组不同的是，矩阵具有许多操作。

但是矩阵的存储方法就是用二维数组存储的。

一般为了方便，可能会定义成一个结构体的形式：

struct mtr {
    int a[N][N];
};

2. 矩阵加法

两个 n\times m 的矩阵 A,B，现在定义矩阵加法如下：

不妨设 A+B=C，则 C_{i,j}=A_{i,j}+B_{i,j}(1\leq i\leq n,1\leq j \leq m)。

如果 A 与 B 的尺寸大小不同，则无法进行加法运算。

矩阵加法具有交换律和结合律。

mtr add(mtr a, mtr b) {
    mtr ans;
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= m; j++) {
            ans.a[i][j] = a.a[i][j] + b.a[i][j];
        }
    }
    return ans;
}

例题：B2104

太过于水，不给代码。

复杂度 \mathcal{O}(nm)。

3. 矩阵减法

矩阵减法是矩阵加法的逆运算。

不妨设 A-B=C，则 C_{i,j}=A_{i,j}-B_{i,j}(1\leq i\leq n,1\leq j \leq m)。

这个也很好理解，复杂度 \mathcal{O}(nm)。

mtr add(mtr a, mtr b) {
    mtr ans;
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= m; j++) {
            ans.a[i][j] = a.a[i][j] - b.a[i][j];
        }
    }
    return ans;
}

4. 矩阵乘法

两个大小分别为 m \times n 和 n \times p 的矩阵 A, B 相乘的结果为一个大小为 m \times p 的矩阵。将结果矩阵记作 C，则

c_{i,j} = \sum\limits_{k = 1}^{n} a_{i,k} \times b_{k,j} \text{,\qquad($1 \le i \le m$, $1 \le j \le p$).}

当然，一般的矩阵乘法都用函数和结构体封装起来会更好。

mtr mul(mtr a, mtr b) {
    mtr ans; for(int i = 1; i <= m; i++) for(int j = 1; j <= p; j++) ans.a[i][j] = 0;
    for(int i = 1; i <= m; i++)
        for(int k = 1; k <= n; k++)
            for(int j = 1; j <= p; j++)
                ans.a[i][j] += a.a[i][k] * b.a[k][j];
    return ans;
}

而如果 A 的列数与 B 的行数不相等，则无法进行乘法。

矩阵乘法满足结合律，即 (A\times B)\times C = A\times(B\times C)。

这里介绍一个 Strassen 矩阵乘法，复杂度 \mathcal{O}(n^{2.81})，略好一筹。

早在清朝的时候，我国的古书《群经评议·周官二》就提到：“凡邦之有疾病者，疕疡者造焉，则使医分而治之，是亦不自医也。”

“分而治之”的思想，在 OI 中是非常重要的。

可以考虑用分块进行优化，使用分块矩阵，没错，将一个矩阵分成四块计算。

递归最简单的情况是 A 和 B 都是 2\times 2 的方阵，A_{1,1} 等分块实际就是一个数，即 2 个 2 阶方阵的乘积可以直接计算出来。

共需要 8 次乘法和 4 次加法。当子矩阵的阶大于 2 时，为求 2 个子矩阵的积，可以继续将子矩阵分块，直到子矩阵的阶降为 2。

由此产生分治降阶的递归算法。

依此算法，计算 2 个 n 阶方阵的乘积转化为计算 8 个 \dfrac{n}{2} 阶方阵的乘积和 4 个 \dfrac{n}{2} 阶方阵的加法。

因此，上述分治法的复杂度 $T(n)$ 应满足： $$ T(n) = \left \{ \begin{aligned} &O(1)\ && n=2 \\ &8T(\dfrac{n}{2})+O(n^2)\ && n>2 \\ \end{aligned} \right. $$ 解得 $T(n)=O(n^3)$。 还是不行，要继续优化。 主要原因是：该方法并没有减少矩阵的乘法次数。而矩阵乘法耗费的时间要比矩阵加法的复杂度高。 所以优化算法的突破口，就是通过加减减少乘法运算。 Strassen 提出可以只进行 $7$ 次乘法运算，具体思想是： $$M_1=A_{1,1}(B_{1,2}-B_{2,2})$$ $$M_2=(A_{1,1}+A_{1,2})B_{2,2}$$ $$M_3=(A_{2,1}+A_{2,2})B_{1,1}$$ $$M_4=A_{2,2}(B_{2,1}-B_{1,1})$$ $$M_5=(A_{1,1}+A_{2,2})(B_{1,1}+B_{2,2})$$ $$M_6=(A_{1,2}-A_{2,2})(B_{2,1}+B_{2,2})$$ $$M_7=(A_{1,1}-A_{2,1})(B_{1,1}+B_{1,2})$$ 做完 $7$ 次乘法运算之后，加减可得 $$C_{1,1}=M_5+M_4+M_6-M_2$$ $$C_{1,2}=M_1+M_2$$ $$C_{2.1}=M_3+M_4$$ $$C_{2,2}=M_1+M_5-M_3-M_7$$ 乍一看好像复杂度更高，其实减少了一次乘法运算。 复杂度 $T(n)$ 应满足： $$ T(n) = \left \{ \begin{aligned} &O(1)\ && n=2 \\ &7T(\dfrac{n}{2})+O(n^2)\ && n>2 \\ \end{aligned} \right. $$ 解得 $T(n)=O(n^{\log_7})=O(n^{2.81})$，略好一筹，但是常数巨大。 其实斯特拉森算法在时效上与暴力差不多，当 $n>100$ 之后，斯特拉森算法略优于暴力。 不过，到目前为止还无法确切知道矩阵乘法的最优时间复杂度，但是复杂度是无论如何也突破不了 $\mathcal{O}(n^2)$。 Coppersmith-Winograd 乘法算法的渐近复杂度最佳，为 $\mathcal{O}(n^{2.37628639})$，但是正如 Sopel 所说，较小的指数不太可能值得在实际计算中使用，因为它们前面有很大的常数。 也正是如此，建议 $n=64$ 了直接暴力算。 例题：[B3615](https://www.luogu.com.cn/problem/B3615) # $5.$ 矩阵快速幂 模板题：[P3390](https://www.luogu.com.cn/problem/P3390) 矩阵想要进行幂的运算，就必须满足行 $=$ 宽，即以 $n\times n$ 的方式表现。 一个大小为 $n \times n$ 的矩阵 $A$ 可以与自身进行乘法，得到的仍是大小为 $n \times n$ 的矩阵，记作 $A^2 = A \times A$。进一步地，还可以递归地定义任意高次方 $A^k = A \times A^{k - 1}$，或称 $A^k = \underbrace{A \times A \times \cdots \times A}_{k \text{ 次}}$。 特殊地，定义 $A^0$ 为单位矩阵 $I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}$。 显然，由于矩阵乘法满足结合律，矩阵的幂也满足 $A^m\times A^n=A^{m+n}$。 不难发现 $IA=A$，所以 $I=\dfrac{A}{A}=A^{1-1}=A^0$。 求 $A^k$，$k$ 太大，可以考虑使用快速幂之类的方案。 矩阵乘法就可以用直接的暴力 $\mathcal{O}(n^3)$ 的求法，复杂度 $\mathcal{O}(n^3\log k)$，可以通过。 ```cpp #include<bits/stdc++.h> using namespace std; #define int long long int n, b; const int N = 110, p = 1e9 + 7; struct mtr {int a[N][N];}; mtr mul(mtr a, mtr b) { mtr ans; for(int i = 1; i <= n; i++) for(int j = 1; j <= n; j++) ans.a[i][j] = 0; for(int i = 1; i <= n; i++) for(int k = 1; k <= n; k++) for(int j = 1; j <= n; j++) ans.a[i][j] = (ans.a[i][j] + (a.a[i][k] % p * b.a[k][j] % p)) % p; return ans; } signed main() { scanf("%lld%lld", &n, &b); mtr A, ans; for(int i = 1; i <= n; i++) for(int j = 1; j <= n; j++) scanf("%lld", &A.a[i][j]); for(int i = 1; i <= n; i++) ans.a[i][i] = 1; while(b) {if(b & 1) ans = mul(ans, A); A = mul(A, A); b >>= 1;} for(int i = 1; i <= n; i++) {for(int j = 1; j <= n; j++) printf("%lld ", ans.a[i][j]); printf("\n");} return 0; } ``` 矩阵快速幂的应用范围很广泛，比如说求斐波那契数的时候。 例题：[P1962](https://www.luogu.com.cn/problem/P1962) 可以发现 $f_n$ 最早从 $f_{n-2}$ 来，所以开一个 $2\times 2$ 的矩阵即可。 第一行想要得到 $f_{n}$，由 $f_{n-1}+f_{n-2}$ 而来，所以第一行为 $[1,1]$。 第二行需要得到 $f_{n-1}$，由 $f_{n-2}+f_{n-3}$ 而来，所以第二行为 $[1,0]$。 答案为矩阵 $\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-1}$ 的头一个元素。 # $6.$ 矩阵的逆元，矩阵除法 想求矩阵除法 $\dfrac{A}{B}$，可以化为 $A\times B^{-1}

那 B^{-1} 怎么求呢？

模板题 P4783

定义：对于一个矩阵 A，需要找到一个矩阵 B，满足 AB=I（I 是单位矩阵），则称 A 是 B 的逆矩阵，B 是 A 的逆矩阵。

求逆矩阵的思路：

求 A 的逆矩阵，把 A 和单位矩阵 I 放在一个矩阵里

对 A 进行加减消元使 A 化成单位矩阵

此时原来单位矩阵转化成逆矩阵。

比如说，求

2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{bmatrix}^{-1}

首先，在右边添加 I，得出

2 & -1 & 0 & 1 & 0 & 0\\ -1 & 2 & -1 & 0 & 1 & 0\\ 0 & -1 & 2 & 0 & 0 & 1 \end{bmatrix}

对左边消元可得

2 & -1 & 0 & 1 & 0 & 0\\\\ 0 & \dfrac{3}{2} & -1 & \dfrac{1}{2} & 1 & 0\\\\ 0 & 0 & \dfrac{4}{3} & \dfrac{1}{3} & \dfrac{2}{3} & 1 \end{bmatrix}

回代，得

2 & 0 & 0 & \dfrac{3}{2} & 1 & \dfrac{1}{2}\\\\ 0 & \dfrac{3}{2} & 0 & \dfrac{3}{4} & \dfrac{1}{2} & \dfrac{3}{4}\\\\ 0 & 0 & \dfrac{4}{3} & \dfrac{1}{3} & \dfrac{2}{3} & 1 \end{bmatrix}

每项除以系数，右半边就是答案。

1 & 0 & 0 & \dfrac{3}{4} & \dfrac{1}{2} & \dfrac{1}{4}\\\\ 0 & 1 & 0 & \dfrac{1}{2} & 1 & \dfrac{1}{2}\\\\ 0 & 0 & 1 & \dfrac{1}{4} & \dfrac{1}{2} & \dfrac{3}{4} \end{bmatrix}

再将分数用费马小定理对 p 取模。

750000006 & 500000004 & 250000002\\\\ 500000004 & 1 & 500000004\\\\ 250000002 & 500000004 & 750000006 \end{bmatrix} 
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int P = 1e9 + 7, N = 410;
int Pow(int a,int b) {
    int ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % P;
        a = a * a % P; b >>= 1;
    }
    return ans;
}
int a[N][N * 2];
signed main() {
    int n, m; scanf("%lld", &n); m = 2 * n;
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) scanf("%lld", &a[i][j]);
        a[i][i + n] = 1;
    }
    for(int i = 1; i <= n; i++) {
        int pos = i;
        for(int j = i + 1; j <= n; j++) if(abs(a[j][i]) > abs(a[pos][i])) pos = j;
        if(i != pos) swap(a[i], a[pos]); if(!a[i][i]) puts("No Solution"), exit(0);
        int inv = Pow(a[i][i], P - 2);
        for(int j = 1; j <= n; j++)
            if(i != j) {
                int x = a[j][i] * inv % P;
                for(int k = i; k <= m; k++)
                    a[j][k] = ((a[j][k] - a[i][k] * x) % P + P) % P;
            }
        for(int j = 1; j <= m; j++) a[i][j] = a[i][j] * inv % P;
    }
    for(int i = 1; i <= n; i++) {
        for(int j = n + 1; j <= m; j++) printf("%lld ", a[i][j]); puts("");
    }
    return 0;
}

求出了矩阵的逆元，矩阵除法就非常简单了，只需要再做一次乘法就可以了。

复杂度 \mathcal{O}(n^3)。

矩阵的逆元有以下几个性质：

所以我们就可以求 A^k 的逆元和 AB 的逆元了。

7. 广义矩阵乘法

假设有矩阵 A:n\times p 和矩阵 B:p \times m，设 C=A\times B 则：

C_{i,j}=\bigotimes_{k=1}^p(A_{i,k}\bigoplus B_{k,j})

若 \otimes 满足对 \oplus 有分配律，则这个运算满足对矩阵的结合律。例如说，乘法满足对加法有分配律，因为 a(b+c)=ab+ac。

在图论中，一种常见的模式是令 \oplus=\min，\otimes=+，进行矩阵运算。事实上这样与 Floyd 算法的代码几乎无异。

例题：P2886

考虑用计数类的 dp。令原图邻接矩阵为 G，令 dp_{u,v,k} 为 u 到 v 长度为 k 的方案数，则可得转移方程

dp_{u,v,k}=\min_{i=1}^{n}{dp_{u,w,k-1}+G_{w,v}}

但是复杂度 \mathcal{O}(n^2k) 太大，不太行。

此时可以用广义矩阵乘法去做。

得 dp_k=G^k，矩阵乘法的 \otimes=\min,\oplus=+，复杂度 \mathcal{O}(n^3\log k)。

但是这里要注意的是需要将点离散化，否则 n\leq 1000，还是会超。

离散化之后 n\le 2T\leq 200,k\leq 10^6，\mathcal{O}(n^3\log k) 还可以接受。

参考代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1010; int id[N];
struct mtr {int a[N][N];} G; int n, m, k, s, t;
mtr mul(mtr a, mtr b) {
    mtr ans; memset(ans.a, 0x3f, sizeof ans.a);
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            for(int k = 1; k <= n; k++)
                ans.a[i][j] = min(ans.a[i][j], a.a[i][k] + b.a[k][j]);
    return ans;
}
mtr Pow(mtr a, int b) {
    mtr ans = a; b--;
    while(b) {
        if(b & 1) ans = mul(ans, a);
        a = mul(a, a); b >>= 1;
    }
    return ans;
}
signed main() {
    cin >> k >> m >> s >> t; memset(G.a, 0x3f, sizeof G.a);
    while(m--) {
        int a, b, c; cin >> c >> a >> b; int u, v;
        if(!id[a]) id[a] = ++n; u = id[a];
        if(!id[b]) id[b] = ++n; v = id[b];
        G.a[u][v] = G.a[v][u] = c;
    }
    s = id[s], t = id[t];
    mtr ans = Pow(G, k);
    cout << ans.a[s][t] << endl;
    return 0;
}

8. 矩阵行列式

模板题见 P7112

矩阵行列式的表示方法：一个 n\times n 的矩阵 A 的行列式记做 \det(A) 或者 |A|。

一个 2\times 2 的矩阵 \begin{bmatrix} a & b\\ c & d\\ \end{bmatrix} 行列式为 ad-bc。

矩阵行列式的定义：\det(A)=|A|=\sum\limits_{p}(-1)^{\tau(p)}\prod a_{i,p_i}，\tau(p) 指 p 的逆序对个数。

把一个 n 阶行列式中的元素 a_{i,j} 所在的第 i 行和第 j 列划去后，留下来的 n-1 阶行列式叫做元素 a_{i,j} 的余子式,记作 M_{i,j}。记 A_{i,j}=(-1)^{i+j}M_{i,j}，叫做元素 a_{i,j} 的代数余子式。

一个 n\times n 矩阵的行列式等于其任意行(或列)的元素与对应的代数余子式乘积之和，即

\det(A)=a_{i,1}A_{i,1}+\cdots +a_{i,n}A_{i,n}=\sum\limits_{j=1}^{n} a_{i.j}(-1)^{i+j}M_{i,j}

矩阵行列式的性质：

接下来，问题是如何求矩阵行列式的值。

还是消元。

我们考虑一个情况，当一个矩阵任意一个位置出现 0，其对行列式的影响非常大，一旦选到 0 就对答案没有任何贡献了。

考虑可以通过消元和前面的性质使得矩阵中出现更多的 0。

显然瞎转化肯定是不行的，我们要让运算次数尽可能少，否则会 TLE。

我们知道求解线性方程组的算法高斯消元。其实高斯消元在做增广矩阵行（初等）变换为行最简形时的步骤和转化相当类似。

我们现在考虑将矩阵一行（列）消成只有最后一个元素非 0 该怎么做。

这里引入代数余子式的概念：

在一个 n 阶行列式 D 中选定 k 行 k 列可以组成一个 k 阶子行列式 A；

删除在 k 行 k 列后剩下的 n-k 阶行列式称为 A 对应的 n-k 阶余子式 M。

对于一个元素 a_{i,j} 的代数余子式为 A_{i,j}，余子式为 M_{i,j}，则

A_{i,j}=(-1)^{i+j}\times M_{i,j}

枚举第 i 行，找到第 i 行某个非 0 元素并将该列与第 i 列交换，然后把 a_{i,i} 变为 1，再用这一行消剩下所有行的第 i 列即可。

然而 a_{i,i} 在模 p 意义下不一定有逆元（p 不一定是质数）。考虑到可以任意相减，这个性质和辗转相除法很相似，可以考虑对两行进行辗转相除，这样一定可以消掉某行第 i 列，很像高斯消元。复杂度 \mathcal{O}(n^3)，卡一卡常，可以冲过去。如果还卡不过去，就把矩阵 A 用指针存储。

参考代码（1.92s 极限冲过）：

#include<bits/stdc++.h>
using namespace std;
#define I inline
#define W while
#define gc getchar
#define pc putchar
namespace SlowIO {
    I int read() {
        int x = 0; char ch = gc();
        W(ch < '0' || ch > '9') ch = gc();
        W(ch >= '0' && ch <= '9') x = x * 10 + (ch ^ 48), ch = gc();
        return x;
    }
    I void Read(int &x) {x = read();}
} using namespace SlowIO;
int p, n; const int N = 610; int *A[N], a[N][N];
int solve() {
    int flag = 1;
    for(int i = 1; i <= n; i++) {
        for(int j = i + 1; j <= n; j++) {
            while(A[j][i]) {
                int t = A[i][i] / A[j][i];
                for(int k = i; k <= n && t; k++)
                    A[i][k] = ((A[i][k] - 1ll * t * A[j][k] % p) % p + p) % p;
                swap(A[i], A[j]);
                flag = -flag;
            }
        }
    }
    int ret = flag;
    for(int i = 1; i <= n; i++) ret = 1ll * ret * A[i][i] % p;
    return (ret + p) % p;
}
int main() {
    cin >> n >> p;
    for(int i = 1; i <= n; i++) A[i] = a[i];
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++) {
            Read(A[i][j]); A[i][j] %= p;
        }
    cout << solve();
    return 0;
}

同时，矩阵树定理也需要用到行列式。

Exercise：

B2104 矩阵加法

B3615 测测你的矩阵乘法

P1962 斐波那契数列

P1939 【模板】矩阵加速（数列）

P2106 Sam数

P1707 刷题比赛

P2044 [NOI2012] 随机数生成器

P1306 斐波那契公约数

P1349 广义斐波那契数列

P5343 【XR-1】分块

P8944 The Da Vinci Code

P6178 【模板】Matrix-Tree 定理