一些些数学的算法笔记
好好好，直接进入正题（
首先我们先要讲讲矩阵，矩阵你可以理解成 \(n\times m\) 的一个二维数组，我们如下表示它：
一、矩阵乘法
1 操作方法
矩阵最重要的一个运算就是矩阵乘法，矩阵乘法就是给定一个 \(n\times m\) 的矩阵 \(A\)，和一个 \(m\times l\) 的矩阵 \(B\)，得到一个 \(n\times l\) 的矩阵 \(C\)，而 \(C\) 的每一位 \(c_{i,j}=\sum_{k=1}^ma_{i,k}\times b_{k,j}\)。简单来说就是 \(c_{i,j}\) 等于 \(A\) 矩阵的第 \(i\) 行，与 \(B\) 矩阵的第 \(j\) 列的元素顺序分别相乘再相加得到的数。
矩阵我们一般会把他封装到一个结构体里，接下来看看一道题吧。
给定 \(n\times n\) 的矩阵 \(A\)，求 \(A^k\)。
首先我们要声明一件事，矩阵乘法满足结合律，但不满足交换律。前者可以直接爆算得出（或者用矩阵就是向量的变换得出），后者可以举例得出。我们知道任何满足结合律的运算都可以用快速幂求解，所以我们只要写出矩阵乘法的代码和快速幂的代码即可。
#include <bits/stdc++.h>
using namespace std;
const int kMaxN = 105, kM = 1e9 + 7;
int n;
long long k;
struct M { // 矩阵结构体
int n, m;
long long a[kMaxN][kMaxN];
M(int h, int w) { n = h, m = w, memset(a, 0, sizeof(a)); }
M() { memset(a, 0, sizeof(a)); } // 这两个都是初始化矩阵
M operator*(M b) { // 重载乘法运算符，实现矩阵乘法
M r(n, b.m);
for (int k = 1; k <= m; k++) {
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= b.m; j++) {
r.a[i][j] = (r.a[i][j] + a[i][k] * b.a[k][j]) % kM;
}
}
}
return r;
}
} A;
M fpow(M A, long long b) { // 矩阵快速幂
M r(n, n);
for (int i = 1; i <= r.n; i++) {
r.a[i][i] = 1;
}
for (; b; b >>= 1, A = A * A) {
if (b & 1) r = r * A;
}
return r;
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> n >> k, A = M(n, n);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cin >> A.a[i][j];
}
}
A = fpow(A, k);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cout << A.a[i][j] << ' ';
}
cout << '\n';
}
return 0;
}
2 例题
介绍完矩阵乘法，你们是不是就要问：“这有什么用呢？”其实这个东西的用处很大，下面再讲一道题。
LG-P1962 斐波那契数列
给定 \(n\)，求斐波那契的第 \(n\) 项模 \(10^9+7\) 得到的值（\(n<2^{63}\)）。
我们设斐波那契数列第 \(i\) 项的值为 \(F_i\)，那么 \(F_i=F_{i-1}+F_{i-2}=1\times F_{i-1}+1\times F_{i-2}\)，考虑构建矩阵 \([F_i\;\;F_{i+1}]\)，若要将其转移成 \([F_{i+1}\;\;F_{i+2}]\) 即 \([1\times F_{i+1}\;\;1\times F_i+1\times F_{i+1}]\)，我们可以使用矩阵乘法进行转移。
所以我们可以通过多次乘系数矩阵，得到 \(F_n\) 的值。由于矩阵乘法满足结合律，所以我们你可以运用矩阵快速幂快速计算 \(F_n\)。
#include <bits/stdc++.h>
using namespace std;
const int kM = 1e9 + 7;
long long n;
struct M {
int n, m;
long long a[5][5];
M(int h, int w) { n = h, m = w, memset(a, 0, sizeof(a)); }
M() { memset(a, 0, sizeof(a)); }
M operator*(M b) {
M r(n, b.m);
for (int k = 1; k <= m; k++) {
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= b.m; j++) {
r.a[i][j] = (r.a[i][j] + a[i][k] * b.a[k][j]) % kM;
}
}
}
return r;
}
} A(1, 2), B(2, 2);
M fpow(M A, long long b) {
M r(2, 2);
for (int i = 1; i <= r.n; i++) {
r.a[i][i] = 1;
}
for (; b; b >>= 1, A = A * A) {
if (b & 1) r = r * A;
}
return r;
}
int main() {
cin >> n;
A.a[1][1] = A.a[1][2] = B.a[1][2] = B.a[2][1] = B.a[2][2] = 1;
A = A * fpow(B, n - 1);
cout << A.a[1][1];
return 0;
}
看完这道题，我们也可以总结出一些规律，我们能够发现一些简单递推，如果范围较大，可以使用矩阵乘法求解，例如：
像上述递推式可以如下求解：
\(B\) 为系数矩阵，自己求。
T234597 路径问题（困难）
有一个 \(n \times m\) 的网格，初始时有个人在 \((1, 1)\) 处，他每次只能向右走一步或者向下走一步，而且不能走出网格。
请问到达 \((n, m)\) 有多少种方案，请输出方案数模 \(10^9 + 7\) 的结果。（\(1 ≤ n ≤ 10^{18}\)，\(1 ≤ m ≤ 100\)）
设 \(F_{i,j}\) 为到达 \((i, j)\) 的方案数，我们很容易推出转移方程：
由于这样的转移方程牵扯到前一行的状态，所以我们将每一行的 \(F\) 作为矩阵表示出来，做相应的矩阵乘法即可。
\(F[i,j]\) 满足下面的递推式:
求 \(F_{i,j}\) 除以 \(10^9+7\) 的余数 （\(1 \le n,m \le 10^{1\,000\,000}\)；\(1 \le a,b,c,d \le 10^9\)）
对于上述递推式，我们先考虑每一行，每一行的变化，就是将上一个数 \(\times a+b\)。考虑初始矩阵为 \([1\;\;F_{i,j}]\)，若要将其变为 \([1\;\;a\times F_{i,j}+b]\)，转移方程为：
那么考虑列的变化的话就是将系数矩阵的 \(b\) 改成 \(d\)，\(a\) 改成 \(c\)，即可。
设
\(A\) 为初始矩阵，\(B,C\) 为系数矩阵。
那么
由于 \(n,m\) 太大，\(n-1,m-1\) 需要用到高精度，快速幂也需要用到 \(10\) 进制快速幂即可。
#include <bits/stdc++.h>
using namespace std;
const int kM = 1e9 + 7;
string n, m;
string J(string s) {
for (int i = s.size() - 1; ~i; i--) {
if (s[i] == '0') {
s[i] = '9';
} else {
s[i]--;
break;
}
}
return s;
}
struct M {
int n, m;
long long a[5][5];
M(int h, int w) { n = h, m = w, memset(a, 0, sizeof(a)); }
M() { memset(a, 0, sizeof(a)); }
M operator*(M b) {
M r(n, b.m);
for (int k = 1; k <= m; k++) {
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= b.m; j++) {
r.a[i][j] = (r.a[i][j] + a[i][k] * b.a[k][j]) % kM;
}
}
}
return r;
}
};
M fpow(M A, string b) {
M r(A.n, A.m), T;
for (int i = 1; i <= r.n; i++) {
r.a[i][i] = 1;
}
for (int i = b.size() - 1; ~i; i--) {
T = M(A.n, A.m);
if (b[i] == '9') {
T = A * A, T = T * T, T = T * T, r = r * T * A;
} else {
for (int j = 1; j <= b[i] - '0'; j++) {
r = r * A;
}
}
A = A * A, T = A, A = A * A, A = A * A, A = A * T;
}
return r;
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
M A(1, 2), B(2, 2), C(2, 2);
cin >> n >> m >> B.a[2][2] >> B.a[1][2] >> C.a[2][2] >> C.a[1][2];
A.a[1][1] = A.a[1][2] = B.a[1][1] = C.a[1][1] = 1;
A = A * fpow(fpow(B, J(m)) * C, J(n)) * fpow(B, J(m));
cout << A.a[1][2];
return 0;
}
二、高斯消元
1 操作方法
学习完矩阵乘法，我们可以发现一个东西，任何一个 \(n\) 元一次方程都可以写成矩阵乘法的形式：
可以写成（上述方程为有 \(n\) 个方程的 \(n\) 原一次方程）：
是不是很神奇，虽然上面的内容和我们下面要讲的东西没有太大关系。（
我们在数学中时常遇到解这种线性方程组，但是元数不高，方程数不多，如果我们要解上述方程该怎么办呢？接下来我们用矩阵（分块矩阵）的思想解决这个问题。
我们设
考虑将 \(A\) 和 \(B\) 塞进一个矩阵得到 \([A\;\;B]\)，接下来引入矩阵的一种变换，它不会改变矩阵的本质，叫做初等行变换。初等行变化有三种变换方式：
- 将矩阵的任意一行中的每一个数乘常数 \(c(c\neq 0)\)。
- 将矩阵的任意两行交换位置。
- 将矩阵的任意一行中的每一个数乘常数 \(c(c\neq 0)\)，在将这一行的每一个数加到另一行对应的数上。
我们发现这三种操作就是方程的性质。
接下来我们将之前的矩阵 \([A\;\;B]\) 进行初等行变换，考虑我们之前学过的解多元一次方程组是如何解的，首先要消元！
我们考虑先消第一个元，如果 \(a_{1,1}=0\) 我们就把下面的一行换上来即可，如果下面全都是 \(0\)，那么就是无解或者无数多组解（后面会提到如何判定这两种情况），这样操作完后，我们就可以将第一行下面的所有行的第一列都消成 \(0\)，然后看第二行，找到第二列，跟第一行做初等行变换消元，然后继续如此操作，继续消，以此类推。不出意外最后我们可能应该也许会消成一个阶梯矩阵，即：
得到阶梯矩阵，我们就可以从第 \(n\) 层开始消，将上一行得到的值代入然后解出所有的未知数。
或者说可以用倒数第二行与倒数第一行进行消元，倒数第三行与倒数第一行和倒数第二行进行消元以此类推，最终可以变为：
这样就可以解得每一个 \(x\) 啦。
