0x34.数学 - 矩阵乘法

声明：
本系列博客是《算法竞赛进阶指南》+《算法竞赛入门经典》+《挑战程序设计竞赛》的学习笔记，主要是因为我三本都买了按照《算法竞赛进阶指南》的目录顺序学习，包含书中的少部分重要知识点、例题解题报告及我个人的学习心得和对该算法的补充拓展，仅用于学习交流和复习，无任何商业用途。博客中部分内容来源于书本和网络（我尽量减少书中引用），由我个人整理总结（习题和代码可全都是我自己敲哒）部分内容由我个人编写而成，如果想要有更好的学习体验或者希望学习到更全面的知识，请于京东搜索购买正版图书：《算法竞赛进阶指南》——作者李煜东，强烈安利，好书不火系列，谢谢配合。

下方链接为学习笔记目录链接（中转站）

学习笔记目录链接

ACM-ICPC在线模板

矩阵

下方截图链接：https://www.cnblogs.com/Mark-X/p/11864285.html#%E7%9F%A9%E9%98%B5
大佬@Baikal Birch总结的很好%%%

对于矩阵 A ，主对角线是指 Ai,iA_{i, i}Ai,i的元素。
一般用 III 来表示单位矩阵，就是主对角线上为 1，其余位置为 0。

矩阵的逆

AAA 的逆矩阵 PPP 是使得A×P=IA \times P = IA×P=I的矩阵。
逆矩阵可以用高斯消元的方式来求。

由于线性递推式可以表示成矩阵乘法的形式，也通常用矩阵快速幂来求线性递推数列的某一项。

设 AAA 是 n×mn×mn×m 的矩阵 BBB 是 m×pm×pm×p 的矩阵，则C=A×BC=A×BC=A×B是n×pn×pn×p的矩，则∀i∈[1,n]\forall i \in[1, n]∀i∈[1,n]，∀j∈[1,p]\forall j\in \lbrack1, p \rbrack∀j∈[1,p]
Ci,j=∑k=1mAi,k×Bk,jC_{i,j} = \sum_{k = 1}^{m} A_{i, k} \times B_{k, j}Ci,j=k=1∑mAi,k×Bk,j

luogu P3390 【模板】矩阵快速幂

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>using namespace std;
typedef long long ll;
const int N = 107, mod = 1e9 + 7;int n, m;
int a[N][N];
int ans[N][N];
ll k;void mul(int c[][N], int a[][N], int b[][N])
{int tmp[N][N] = {0};for(int i = 1; i <= n; ++ i){for(int j = 1; j <= n; ++ j){for(int k = 1; k <= n; ++ k){tmp[i][j] = (tmp[i][j] + (ll)a[i][k] * b[k][j]) % mod;}}}memcpy(c, tmp, sizeof tmp);
}int main(){scanf("%d%lld", &n, &k);for(int i = 1; i <= n; ++ i)for(int j = 1; j <= n; ++ j)scanf("%d", &a[i][j]);for(int i = 1; i <= n; ++ i)ans[i][i] = 1;while(k){if(k & 1) mul(ans, ans, a);mul(a, a, a);k >>= 1;}for(int i = 1; i <= n;puts(""), ++ i)for(int j = 1; j <= n; ++ j)printf("%d ", ans[i][j]);return 0;
}

AcWing 205. 斐波那契

具体解析见《算法竞赛进阶指南》P157

时间复杂度：O(23logn)O(2^3logn)O(23logn)，其中232^323是矩阵乘法的复杂度，如果矩阵为n×nn\times nn×n，那么矩阵乘法的复杂度为O(n3)O(n^3)O(n3)

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>using namespace std;
const int N = 5007, mod = 10000;int n, m;void mul(int f[2], int a[2][2])
{int c[2];memset(c, 0, sizeof c);//因为f只有一维，相当于只有j,k没有ifor(int j = 0; j < 2; ++ j){for(int k = 0; k < 2; ++ k){c[j] = (c[j] + (long long)f[k] * a[k][j]) % mod;}}memcpy(f, c, sizeof c);
}void mulself(int a[2][2])
{int c[2][2];memset(c, 0, sizeof c);for(int i = 0; i < 2; ++ i){for(int j = 0; j < 2; ++ j){for(int k = 0; k < 2; ++ k){c[i][j] = (c[i][j] + (long long)a[i][k] * a[k][j]) % mod;}}}memcpy(a, c, sizeof c);
}int main()
{while(scanf("%d", &n) != EOF && n != -1){int f[2] = {0, 1};int a[2][2] = {{0, 1}, {1, 1}};while(n){if(n & 1)mul(f, a);n >>= 1;mulself(a);}printf("%d\n", f[0]);}return 0;
}

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

https://www.luogu.com.cn/problem/P1939

AcWing 206. 石头游戏

石头游戏在一个 nnn 行 mmm 列的网格上进行，每个格子对应一种操作序列，操作序列至多有 101010 种，分别用 0∼90\sim90∼9 这 101010 个数字指明。

操作序列是一个长度不超过 666 且循环执行、每秒执行一个字符的字符串。

每秒钟，所有格子同时执行各自操作序列里的下一个字符。

序列中的每个字符是以下格式之一：

数字 0∼9：表示拿 0∼90\sim90∼9 个石头到该格子。
NWSE：表示把这个格子内所有的石头推到相邻的格子，NNN 表示上方，WWW 表示左方，SSS 表示下方，EEE 表示右方。
DDD：表示拿走这个格子的所有石头。
给定每种操作序列对应的字符串，以及网格中每个格子对应的操作序列，求石头游戏进行了 ttt 秒之后，石头最多的格子里有多少个石头。

在游戏开始时，网格是空的。

Solution

由于每一个点的操作周期是不一样的，考虑全部扩大到他们的 lcm(1,2,3,4,5,6)=60\mathrm{lcm}(1,2,3,4,5,6)=60lcm(1,2,3,4,5,6)=60。

构造 606060 个转移矩阵，这样先把所有转移矩阵乘起来转移 ⌊t60⌋\left\lfloor\cfrac t {60}\right\rfloor⌊60t⌋ 次，剩下的直接暴力乘 606060 次就行了。

Code

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int maxn = 10 + 7, maxm = 60 + 7, maxs = 10 + 7;
#define num(i, j) ((i - 1) * m + j)
int n, m, t, act;
ll ans;
int LCM, op[maxm][maxm], len[maxs];
char s[maxn][maxs];struct Matrix
{int n, m;ll a[maxm][maxm];Matrix(){memset(a, 0, sizeof a);}Matrix operator * (const Matrix &t) const {Matrix C;C.n = n, C.m = t.m;//memset(C.a, 0, sizeof C.a);for (int i = 0; i <= n; ++ i) for (int j = 0; j <= C.m; ++ j) for (int k = 0; k <= m; ++ k) C.a[i][j] = C.a[i][j] + a[i][k] * t.a[k][j];return C;}
} A, F, a[maxm];int lcm(int a, int b)
{return a / __gcd(a, b) * b;
}int main()
{scanf("%d%d%d%d", &n, &m, &t, &act);for (int i = 1; i <= n; ++ i) for (int j = 1; j <= m; ++ j) {/*char c;cin >> c;op[i][j] = c - '0';*/scanf("%1d", &op[i][j]);}LCM = 1;for (int i = 0; i < act; ++ i) {scanf("%s", s[i] + 1);len[i] = strlen(s[i] + 1);LCM = lcm(LCM, len[i]);}  for (int k = 1; k <= LCM; ++ k) {a[k].a[0][0] = 1;a[k].n = n * m, a[k].m = n * m;for (int i = 1; i <= n; ++ i) {for (int j = 1; j <= m; ++ j) {int l = len[op[i][j]];if (l == 0) {cout << i << j << endl;cout << op[i][j] << endl;continue;}char c = s[op[i][j]][((k - 1) % l) + 1];if(c >= '0' && c <= '9') a[k].a[0][num(i, j)] = c - '0', a[k].a[num(i, j)][num(i, j)] = 1;else if(c == 'N' && i > 1) a[k].a[num(i, j)][num(i - 1, j)] = 1;else if(c == 'S' && i < n) a[k].a[num(i, j)][num(i + 1, j)] = 1;else if(c == 'W' && j > 1) a[k].a[num(i, j)][num(i, j - 1)] = 1;else if(c == 'E' && j < m) a[k].a[num(i, j)][num(i, j + 1)] = 1; }}} int r = t % LCM;F.n = 0, F.m = n * m, F.a[0][0] = 1;A = a[1];for (int i = 2; i <= LCM; ++ i)A = A * a[i];int q = t / LCM;while (q) {if (q & 1) F = F * A;A = A * A;q >>= 1;}ans = 0;for (int i = 1; i <= r; ++ i)F = F * a[i];for (int i = 1; i <= n * m; ++ i)ans = max(ans, F.a[0][i]);cout << ans << endl; return 0;
}