一些些数学的算法笔记
好好好，直接进入正题（
首先我们先要讲讲矩阵，矩阵你可以理解成 \(n\times m\) 的一个二维数组，我们如下表示它：
矩阵最重要的一个运算就是矩阵乘法，矩阵乘法就是给定一个 \(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\) 列的元素顺序分别相乘再相加得到的数。
矩阵我们一般会把他封装到一个结构体里，接下来看看一道题吧。
LG-P3390 【模板】矩阵快速幂
给定 \(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;
}
介绍完矩阵乘法，你们是不是就要问：“这有什么用呢？”其实这个东西的用处很大，下面再讲一道题。
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\) 作为矩阵表示出来，做相应的矩阵乘法即可。
P1397 [NOI2013] 矩阵游戏
\(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;
}