0.前言
感谢远古神猴tz1带来的30分钟极速版矩阵乘法讲解,tql!
1.矩阵作用
矩阵本质是对一个向量的进行变换,也就是描述转移的东西,因此我们常常用其来加速转移过程。
2.技巧
让我们来结合几道题目来谈谈吧。
2.1 优化矩阵运算
[CF1970E3] Trails
首先,暴力矩阵转移是显然的,记\(p_i = l_i + s_i\) , 我们可以很轻松得到以下式子
\[dp^{'}_i = \sum\limits_j dp_j \times (p_i p_j - l_i l_j) \]然而转移矩阵太大,我们考虑简化。
对于转移矩阵
是等价于
\[ \left[\begin{array}{cccc} p_1 & l_1\\ p_2 & l_2\\ \vdots & \vdots \\ p_n & l_n \end{array}\right] \left[\begin{array}{ccc} p_1 & p_2 & \cdots & p_n \\ -l_1 & -l_2 & \cdots & -l_n \\ \end{array}\right] \]不妨分别记作\(A , B\)
但是\((AB)^n\)仍然没法算,此时可以利用奇淫技巧,把\((AB)^n\)改写为\(A\times(BA)^{n - 1}\times B\)
\(BA\)是一个\(2 \times 2\) 的矩阵,可以矩阵快速幂计算。
2.2 矩阵加速转移
这个用的比较多,具有代表性的有交换值。
[THUSCH2017] 大魔法师
算是矩阵运算的一个大杂烩,有点专门恶心人的意思。
当你看见三种基本等地位的变量与变量之间互相转移时,矩阵便呼之欲出了。
可以注意到每一种操作都可以比较容易的用矩阵描述,在实际中多开一个位置存常量控制区间加,乘,赋值操作,然后在线段树上每个位置开一个\(lazy\)矩阵标记即可。
巨丑代码
#include <bits/stdc++.h>
#define SG Segment_Tree
using namespace std;
const int Maxn = 2.5e5 + 3 , Mod = 998244353;
int n , A[Maxn][3];
struct Matrix {
long long nums[4][4] ;
int n , m ;
void init() {
n = m = 0;
memset(nums, 0 , sizeof nums);
}
Matrix operator* (const Matrix& cur) {
Matrix res;res.init();
res.n = n , res.m = cur.m;
for(int i = 0 ; i < res.n ; i++) {
for(int j = 0 ; j < res.m ; j++) {
for(int k = 0 ; k < m ; k++) {
res.nums[i][j] = (res.nums[i][j] + 1ll * nums[i][k] * cur.nums[k][j]) % Mod;
}
}
}
return res;
}
Matrix operator+ (const Matrix& cur) {
Matrix res;res.init();
res.n = n , res.m = m;
for(int i = 0 ; i < res.n ; i++) {
for(int j = 0 ; j < res.m ; j++) {
res.nums[i][j] = (nums[i][j] + cur.nums[i][j]) % Mod;
}
}
return res;
}
}Unit , Spd;
namespace Segment_Tree {
Matrix trans[Maxn << 2] , val[Maxn << 2];
void Push_Up(int p) {
for(int i = 0 ; i < 3 ; i++) {
val[p].nums[0][i] = (val[p << 1].nums[0][i] + val[p << 1 | 1].nums[0][i]) % Mod;
}
}
void build(int l = 1 , int r = n ,int p = 1) {
trans[p] = Unit , val[p].init();
val[p].n = 1 , val[p].m = 4;
val[p].nums[0][3] = r - l + 1;
if(l == r) {
for(int j = 0 ; j < 3 ; j++) {
val[p].nums[0][j] = A[l][j];
}
return ;
}int mid =(l + r) >> 1;
build(l , mid , p << 1) , build(mid + 1 , r, p << 1 | 1);
Push_Up(p);
}
void change(int p , Matrix& cur) {
val[p] = val[p] * cur;
trans[p] = trans[p] * cur;
}
void Push_Down(int p) {
bool f =0;
for(int i = 0 ; i < 4 ; i++) {
for(int j = 0 ; j< 4; j++) {
if(trans[p].nums[i][j] != Unit.nums[i][j]) {
f = 1;
}
}
}
if(f) {
change(p << 1 , trans[p]);
change(p << 1 | 1 , trans[p]);
trans[p] = Unit;
}
}
inline void sec(int tarl ,int tarr , Matrix& cur , int l = 1 , int r = n , int p = 1) {
if(tarl == l && tarr == r) {
change(p , cur);
// printf("-> %d %d %d\n" , val[p].nums[0][0] , val[p].nums[0][1] , val[p].nums[0][2]);
return ;
}int mid = (l + r) >> 1;
Push_Down(p);
if(tarr <= mid) sec(tarl , tarr , cur , l , mid , p << 1);
else if(tarl > mid) sec(tarl , tarr , cur , mid + 1 , r , p << 1 | 1);
else sec(tarl , mid , cur , l , mid , p << 1) , sec(mid + 1 , tarr , cur, mid + 1 , r , p << 1 | 1);
Push_Up(p);
}
inline void Mul(int tarl , int tarr , int op , int v , int l = 1 , int r = n , int p = 1) {
if(tarl == l && tarr == r) {
if(op == 1) {
// puts("SA");
trans[p].nums[3][0] = (trans[p].nums[3][0] + v) % Mod;
val[p].nums[0][0] = (val[p].nums[0][0] + 1ll * v * (r - l + 1)) % Mod;
}else if(op == 2) {
for(int i = 0 ; i < 4; i++) trans[p].nums[i][1] = (1ll * trans[p].nums[i][1] * v) % Mod;
val[p].nums[0][1] = (1ll * val[p].nums[0][1] * v) % Mod;
}else {
for(int i = 0 ; i < 4; i++) trans[p].nums[i][2] = 0;
trans[p].nums[3][2] = v;
val[p].nums[0][2] = (1ll * v * (r - l + 1)) % Mod;
}
return ;
}int mid = (l + r) >> 1;
Push_Down(p);
if(tarr <= mid) Mul(tarl , tarr , op , v , l , mid , p << 1);
else if(tarl > mid) Mul(tarl , tarr , op , v , mid + 1 , r , p << 1 | 1);
else Mul(tarl , mid , op , v , l , mid , p << 1) , Mul(mid + 1 , tarr , op , v, mid + 1 , r , p << 1 | 1);
Push_Up(p);
}
inline Matrix ask(int tarl , int tarr , int l = 1 , int r = n , int p = 1) {
if(tarl == l && tarr == r) {
return val[p];
}int mid =(l + r ) >> 1;
Push_Down(p);
if(tarr <= mid) return ask(tarl , tarr , l , mid , p << 1);
else if(tarl > mid) return ask(tarl , tarr , mid + 1 , r , p << 1 | 1);
else return (ask(tarl , mid , l , mid, p << 1) + ask(mid + 1 , tarr , mid + 1 , r , p << 1 | 1));
}
}
inline int Read() {
int x = 0;
char sign = getchar();
bool f =0;
for(; sign < '0' || sign > '9' ; sign = getchar()) f |= (sign == '-');
for(; sign <= '9' && sign >= '0' ; sign = getchar()) x = (x << 1) + (x << 3) + (sign ^ 48);
return f ? -x : x;
}
int main() {
Unit.init();
Unit.n = 4 , Unit.m = 4;
for(int i= 0 ; i< 4 ; i++) {
Unit.nums[i][i] = 1;
}
n = Read();
for(int i= 1 ; i <= n ; i++) {
for(int j = 0 ; j < 3 ; j++) {
A[i][j] = Read();
}
}SG::build();
int q , op , l , r , v;
q = Read();
while(q--) {
op = Read() , l = Read() , r = Read(), Spd.init();Spd.n = Spd.m = 4;
if(op == 1) {
Spd.nums[0][0] = 1 , Spd.nums[1][0] = 1 , Spd.nums[1][1] = 1 , Spd.nums[2][2] = 1 , Spd.nums[3][3] = 1;
SG::sec(l , r , Spd);
}else if(op == 2){
Spd.nums[0][0] = 1 , Spd.nums[1][1] = 1 , Spd.nums[2][1] = 1 , Spd.nums[2][2] = 1 , Spd.nums[3][3] = 1;
SG::sec(l , r , Spd);
}else if(op == 3) {
Spd.nums[0][0] = 1 , Spd.nums[1][1] = 1 , Spd.nums[0][2] = 1 , Spd.nums[2][2] = 1 , Spd.nums[3][3] = 1;
SG::sec(l , r , Spd);
}else if(op == 4) {
v =Read();
SG::Mul(l , r , 1 , v);
}else if(op == 5) {
v = Read();
SG::Mul(l , r , 2 , v);
}else if(op == 6) {
v = Read();
SG::Mul(l , r , 3 , v);
}else {
Matrix res = SG::ask(l , r);
printf("%d %d %d\n" , res.nums[0][0],res.nums[0][1] , res.nums[0][2]);
}
}
}
[CF1458C] Latin Square
最主要的就是这个逆排列的操作。
不妨记操作前后排列为\(p ,q\),记$S_1 = \{ (i , p_i) \} $ , $S_2 = \{ (q_i , i) \} $ , 不难发现,\(S_1 = S_2\), 这说明对于你排列操作,只是把值与下标互换,这可以用矩阵实现,而且其他操作用矩阵同样较为方便,对于一个三元组\((x, y ,z)\)表示位置\((x , y)\)的值为\(z\),我们对其进行矩阵变换,可以先预处理出所有操作矩阵的叠加矩阵,再把所有的\((x , y , z)\)一个一个套进去就可以,实现时可以多开一位常量处理平移,边界用模\(n\)实现头和尾的交换。
巨丑代码
#include <bits/stdc++.h>
#define SG Segment_Tree
using namespace std;
const int Maxn = 1000 + 3 , inf = 0x3f3f3f3f;
struct Matrix {
int nums[5][5] ;
int n , m;
void init() {
memset(nums , 0 ,sizeof nums);
}
}Speed;
int T;
int n , m , a[Maxn][Maxn] , p[Maxn][Maxn];
Matrix operator*(const Matrix& A , const Matrix& B) {
Matrix res;
res.init();
res.n = A.n , res.m = B.m;
for(int i = 0 ; i < A.n ; i++)
{for(int k = 0 ; k < A.m ; k++) {
int c = A.nums[i][k];
for(int j = 0 ; j < B.m ; j++) {
res.nums[i][j] = (res.nums[i][j] + 1ll * c * B.nums[k][j] - 1) % n + 1;
}
}
}
for(int i = 0 ;i < A.n ; i++) {
for(int j = 0 ; j < B.m ; j++) {
res.nums[i][j] = (res.nums[i][j] + n - 1) % n + 1;
}
}
return res;
}
char s[100005];
int main() {
scanf("%d" , &T);
while(T--) {
scanf("%d%d" , &n, &m);
for(int i = 1 ; i <= n ; i++) {
for(int j = 1 ; j<= n ; j++) {
scanf("%d" ,&a[i][j]);
}
}Speed.init();
Speed.n =Speed.m = 4;
for(int i = 0 ; i < 4 ;i++) Speed.nums[i][i] = 1;
scanf("%s" ,s + 1);
for(int i = 1 ; i<= m ; i++) {
Matrix tmp;
tmp.init();
tmp.n = tmp.m = 4;
for(int i =0 ;i < 4 ; i++) tmp.nums[i][i] = 1;
if(s[i] == 'U') {
tmp.nums[3][0] = -1;
}else if(s[i] == 'D') {
tmp.nums[3][0] = 1;
}else if(s[i] == 'L') {
tmp.nums[3][1] = -1;
}else if(s[i] == 'R') {
tmp.nums[3][1] = 1;
}else if(s[i] == 'I') {
swap(tmp.nums[1][1] , tmp.nums[1][2]);
swap(tmp.nums[2][2] , tmp.nums[2][1]);
}else {
swap(tmp.nums[0][0] , tmp.nums[2][0]);
swap(tmp.nums[0][2] , tmp.nums[2][2]);
}
Speed = Speed * tmp;
}
for(int i = 1 ;i <= n ; i++) {
for(int j = 1 ; j <= n ; j++) {
Matrix cur;
cur.n = 1 , cur.m = 4;
cur.nums[0][0] = i , cur.nums[0][1] = j , cur.nums[0][2] = a[i][j] , cur.nums[0][3] = 1;
cur = cur * Speed;
p[cur.nums[0][0]][cur.nums[0][1]] = cur.nums[0][2];
}
}
for(int i = 1 ;i <= n ; i++) {
for(int j = 1 ; j <= n ; j++) {
printf("%d " , p[i][j]);
}puts("");
}puts("");
}
}
2.3 一类特殊运用——邻接矩阵
设在图上的邻接矩阵的\(n\)次自乘后得到的矩阵为\(A\) ,则\(A_{i , j}\) 为\(i\) 到\(j\) 的所有\(n\)步不同路径(不一定是简单路径)
HH去散步
考虑到边很少,我们不妨对一条边开出两个新点,表示从这条边的另一点到当前点的状态,可以走到除了该边的其他边连接的点,再开\(n\)个孤点表示起点(没有上一个点)。总点数不超过\(200\) ,可以解决问题
code
#include <bits/stdc++.h>
using namespace std;
const int Maxn = 200+ 3 , Mod = 45989;
const long long inf = 1e18;
struct Matrix {
long long nums[Maxn][Maxn] ;
int n , m;
void init() {
memset(nums , 0 ,sizeof nums);
}
}Speed ;
Matrix operator* (const Matrix& A , const Matrix& B) {
Matrix res;res.init();
res.n =A.n , res.m = B.m;
for(int i = 1 ; i <= A.n ; i++) {
for(int k = 1 ; k <= A.m ; k++) {
int c = A.nums[i][k];
for(int j = 1 ; j <= B.m ; j++) {
res.nums[i][j] = (res.nums[i][j] + 1ll* c * B.nums[k][j]) % Mod;
}
}
}
return res;
}
int T , n, m, A , B , cnt;
struct Node {
int u , v;
}a[Maxn];
vector < pair < int , int > > vt[Maxn];
map < pair < int , int > , int > mp;
pair < int , int > rfs[Maxn];int ans = 0;
void qpow(int x) {
Matrix res;res.init();
res.n = res.m = cnt;
for(int i = 1; i <= cnt ; i++) res.nums[i][i] = 1;
while(x) {
if(x & 1) {
res = res * Speed;
}
Speed = Speed * Speed;
x >>= 1;
}
for(int i = 1; i <= cnt ; i++) {
if(rfs[i].second == B) {
ans = (ans + res.nums[mp[make_pair(0 , A)]][i]) % Mod;
}
}
}
int main() {
scanf("%d%d%d%d%d" ,&n , &m , &T , &A , &B);
for(int i= 1 ;i <= m ; i++) {
scanf("%d%d" , &a[i].u , &a[i].v);
vt[a[i].u].push_back(make_pair(a[i].v , i));
vt[a[i].v].push_back(make_pair(a[i].u , i));
mp[make_pair(i , a[i].v)] = ++cnt;
rfs[cnt] = make_pair(i , a[i].v);
mp[make_pair(i , a[i].u)] = ++cnt;
rfs[cnt] = make_pair(i , a[i].u);
}
for(int i = 0; i < n ; i++) {
mp[make_pair(0 , i)] = ++cnt;
rfs[cnt] = make_pair(0, i);
}Speed.n = Speed.m = cnt;
Speed.init();
for(int i = 1 ; i <= cnt ; i++) {
int u = rfs[i].second;
for(auto x : vt[u]) {
if(x.second == rfs[i].first) {
continue;
}
Speed.nums[i][mp[make_pair(x.second , x.first)]] = 1;
}
}qpow(T);
printf("%d\n" , ans);
}
[CF593E] Strange Calculation and Cats
其实与上一题差不了多少,这里需要实时维护一下邻接矩阵即可
code
#include <bits/stdc++.h>
#define SG Segment_Tree
using namespace std;
const int Maxn = 20 + 3 , inf = 0x3f3f3f3f , Mod = 1e9 + 7;
struct Matrix {
int nums[Maxn][Maxn] ;
int n , m;
void init() {
memset(nums , 0 ,sizeof nums);
}
}Speed ,res;
int T;
int n , m;
Matrix operator*(const Matrix& A , const Matrix& B) {
Matrix res;
res.init();
res.n = A.n , res.m = B.m;
for(int i = 1 ; i <= A.n ; i++)
{for(int k = 1 ; k <= A.m ; k++) {
int c = A.nums[i][k];
for(int j = 1 ; j <= B.m ; j++) {
res.nums[i][j] = (res.nums[i][j] + 1ll * c * B.nums[k][j]) % Mod;
}
}
}
return res;
}
inline void qpow(int times) {
while(times) {
if(times & 1) {
res = res * Speed;
}
Speed = Speed * Speed;
times >>= 1;
}
}
int Name[Maxn][Maxn] ,cnt , op , x , y , p , path[Maxn][Maxn] ,a[Maxn][Maxn];
const int dx[4] = {1 , 0 , -1 , 0} , dy[4] = {0 , 1 , 0 , -1};
int main() {
scanf("%d%d%d", &n , &m,&T );
for(int i = 1 ; i <= n ; i++) {
for(int j = 1 ; j <= m ; j++) {
Name[i][j] = ++cnt;
}
}
for(int i = 1 ; i <= n ; i++) {
for(int j = 1 ; j <= m ; j++) {
path[Name[i][j]][Name[i][j]] = 1;
for(int k = 0 ; k < 4 ; k++) {
int tx = i + dx[k] , ty = j + dy[k];
if(tx < 1 || tx > n || ty < 1 || ty > m) continue;
path[Name[i][j]][Name[tx][ty]] = 1;
}
}
}
int lai = 1;
res.n = res.m = cnt;
for(int i = 1 ; i <= cnt; i++) {
res.nums[i][i] = 1;
}
while(T--) {
scanf("%d%d%d%d" , &op ,&x , &y , &p);
Speed.n = Speed.m = cnt;
for(int i = 1 ; i <= cnt ; i++) {
for(int j = 1 ; j <= cnt ;j++) {
Speed.nums[i][j] = path[i][j];
}
}
qpow(p - lai);lai = p;
if(op == 1) {
printf("%d\n" , res.nums[Name[1][1]][Name[x][y]]);
}else if(op == 2){
for(int i = 0 ; i < 4 ; i++) {
int tx = x + dx[i] , ty = y + dy[i];
if(tx < 1 || tx > n || ty < 1 || ty > m) {
continue;
}
path[Name[x][y]][Name[tx][ty]] = path[Name[tx][ty]][Name[x][y]] = 0;
}
path[Name[x][y]][Name[x][y]] = 0;
a[x][y] = 1;
}else {
a[x][y] = 0;
for(int i = 0 ; i < 4 ; i++) {
int tx = x + dx[i] , ty = y + dy[i];
if(tx < 1 || tx > n || ty < 1 || ty > m) {
continue;
}
if(!a[x][y] && !a[tx][ty]) path[Name[x][y]][Name[tx][ty]] = path[Name[tx][ty]][Name[x][y]] = 1;
} path[Name[x][y]][Name[x][y]] = 1;
}
}
}
未完待续......
标签:Matrix,nums,int,res,矩阵,Maxn,浅谈 From: https://www.cnblogs.com/Cmr-/p/18346087