矩阵
概念
定义 矩形的数组
特殊矩阵
零矩阵：所有元素均为0的矩阵
单位矩阵（E或I）：对角线元素为1，其他元素为0的矩阵（A*E=A）
方阵：列数与行数相等的矩阵
对称矩阵：转置矩阵与原矩阵相同的矩阵
在编程求解有关矩阵的题目时，可以使用下面的模板：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
| template<typename T>//使用泛型以适用不同的数据类型
struct MATRIX
{
T ary[maxn][maxn];
T* operator[](int pos)//这样定义后可以用matrix[i][j]直接访问矩阵元素
{
return ary[pos];
}
void clear()
{
for(int i=0;i<=n;++i)
{
for(int j=0;j<=n;++j)
{
ary[i][j]=0;
}
}
}
MATRIX()//结构体的构造函数
{
memset(ary,0,sizeof(ary));
}
void sete()//初始化为单位矩阵
{
clear();
for(int i=1;i<=n;++i)
ary[i][i]=1;
}
};
|
基本操作
初等变换
行初等变换
- 交换两行
- 用k乘i行（k!=0）
- i+=j*k(k!=0)
列初等变换类似。
转置
交换一个矩阵的行和列。
转置后得到的矩阵称为原矩阵的转置矩阵。
基本运算
矩阵加（减）法
概念
行列数相同的矩阵A、B，
C=A+B定义为C[i][j]=A[i][j]+B[i][j]
代码
1
2
3
4
5
6
7
8
| MATRIX operator+(MATRIX a,MATRIX b)
{
MATRIX ans;
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
ans[i][j]=(a[i][j]+b[i][j])%mod;
return ans;
}
|
矩阵乘法
概念
矩阵A(n*r)、B(r*m)
定义C(n*m)=A*B，
则有C[i][j]=∑A[i][k]*B[k][j],k=1~r
性质
l 若A(n*r)、B(r*m)，则相乘得到C(n*m)。
l 满足结合律但不满足交换律。
手算
如图，将相乘的两个矩阵A，B如图放置，对于结果矩阵中一个元素（紫色格子），它的值为A矩阵中同一行从左到右（红色部分）与B矩阵中同一列从上到下（蓝色部分）依次相乘后相加。
代码
1
2
3
4
5
6
7
8
9
| MATRIX operator*(MATRIX a,MATRIX b)
{
MATRIX ans;
for(int i=1;i<=n;++i)
for(int j=1;j<=m;++j)
for(int k=1;k<=r;++k)
ans[i][j]+=a[i][k]*b[k][j]%mod,ans[i][j]%=mod;
return ans;
}
|
逆矩阵
矩阵A的逆矩阵B，满足A*B=B*A=E
求法：
将A与E 写到一起，形如[A]|[E]
$$
\left ( \begin{bmatrix}
a & 0\\
1 & b
\end{bmatrix} |\begin{bmatrix}
1 & 0\\
0 & 1
\end{bmatrix}\right )
$$
而后利用矩阵初等变换将左边变换为单位矩阵，记总的矩阵变换为f
结果形如[E]|[B]，则B为A的逆矩阵。
应用
求解线性方程组
将线性方程组化简为：
a1x1+a2x2+a3x3+…+amxm=y1
b1x1+b2x2+b3x3+…+bmxm=y2
…
k1x1+k2x2+k3x3+…+kmxm=yn
之后，该方程组的解可以通过矩阵的方法来求。
几个概念：
系数矩阵：将所有方程左边的未知数系数按顺序排列得到的矩阵，如下：
$$
\begin{bmatrix}
a1 & a2 & ... & am\\
b1 & b2 & ... & bm\\
& ... & ... & \\
k1 & k2 & ... & km
\end{bmatrix}
$$
增广矩阵：将方程右边的常数项也加入系数矩阵中得到的矩阵。
$$
\begin{bmatrix}
a1 & a2 & ... & am &y1\\
b1 & b2 & ... & bm &y2\\
& ... & ... & &...\\
k1 & k2 & ... & km &ym
\end{bmatrix}
$$
首元：每个行向量中第一个非零元素。
解方程组：利用矩阵的行初等变换，对矩阵进行化简，最终到达这样一种状况：
如果该方程组有解：每一行的首元位置依次递增，且除了常数项和首元外，其他的元素均为0。
矩阵快速幂
基本同快速幂，一般使用非递归式。
注意：由矩阵乘法的条件，能进行矩阵快速幂的矩阵一般为方阵。
1
2
3
4
5
6
7
8
9
10
11
| MATRIX fpow(MATRIX w,lint p)
{
MATRIX ans;
for(ans.sete();p;p/=2)//sete():将矩阵初始化为单位矩阵
{
if(p&1)
ans=ans*w;
w=w*w;
}
return ans;
}
|
加速数列线性递推
利用矩阵快速幂，可以将O(n)的线性数列递推优化到O(logn)。
基本步骤
根据数列递推式构造矩阵
不断递推得到矩阵的幂指数
矩阵快速幂
计算答案
分步解析
由于能够进行矩阵快速幂的矩阵一般为方阵，所以我们只需要确定方阵的边长。
一般而言与方阵的大小有关的是递推式右边不同的项数，不过题目给出的初始值项数有时也会有一定提示。
一般构造如下所示的等式：
$$
\begin{bmatrix}
f[ n]\\
f[ n-1]\\
...\\
f[ n-k]
\end{bmatrix} =[ A] \times \begin{bmatrix}
f[ n-k-1]\\
f[ n-k-2]\\
\\\\
f[ n-k-( k+1)]
\end{bmatrix}
$$
其中A为根据递推式构造的k*k的矩阵。
至于k的选取，一般要求能够使得f[n]~f[n-k]能表示为
f[n-i]=a1*f[n-k-1]+a2*f[n-k-2]+...+ak*f[n-k-(k+1)]
的形式，这样才能将a1~ak填入A中第i行。
例1
对于斐波那契数列，f[1]=f[2]=1,f[n]=f[n-1]+f[n+1]，我们可以这样构造矩阵：
$$
\begin{bmatrix}
f[ n]\\
f[ n-1]
\end{bmatrix} =\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix} \times \begin{bmatrix}
f[ n-1]\\
f[ n-2]
\end{bmatrix}
$$
而矩阵 $\begin{bmatrix}
f[ n-1]\\
f[ n-2]
\end{bmatrix}$ 也可以用同样的方法展开，有：
$$
\begin{gathered}
\begin{bmatrix}
f[ n]\\
f[ n-1]
\end{bmatrix} =\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix} \times \begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix} \times \begin{bmatrix}
f[ n-2]\\
f[ n-3]
\end{bmatrix}\\
=\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}^{2} \times \begin{bmatrix}
f[ n-2]\\
f[ n-3]
\end{bmatrix}\\
=\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}^{k} \times \begin{bmatrix}
f[ n-k]\\
f[ n-k-1]
\end{bmatrix}\\
=\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix}^{n-2} \times \begin{bmatrix}
f[ 2]\\
f[ 1]
\end{bmatrix}\\
\end{gathered}
$$
代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
| /*
斐波那契数列矩阵加速
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
using namespace std;
typedef long long lint;
lint n;
const lint mod=lint(1e9+7);
struct MATRIX
{
lint ary[3][3];
void clear()
{
for(int i=0;i<=2;++i)
for(int j=0;j<=2;++j)
ary[i][j]=0;
}
MATRIX()
{
clear();
}
void init()
{
ary[1][1]=1,ary[1][2]=1,ary[2][1]=1,ary[2][2]=0;
}
lint* operator[](int pos)
{
return ary[pos];
}
void sete()
{
for(int i=1;i<=2;++i)
{
for(int j=1;j<=2;++j)
{
ary[i][j]=lint(i==j);
}
}
}
}matrix;
MATRIX operator*(MATRIX a,MATRIX b)
{
MATRIX ans;
for(int i=1;i<=2;++i)
for(int j=1;j<=2;++j)
for(int k=1;k<=2;++k)
ans[i][j]+=a[i][k]*b[k][j],ans[i][j]%=mod;
return ans;
}
MATRIX operator%(MATRIX a,lint mo)
{
for(int i=1;i<=2;++i)
for(int j=1;j<=2;++j)
a[i][j]%=mo;
return a;
}
MATRIX pow(MATRIX w,lint p)
{
MATRIX ans;
for(ans.sete();p;p/=2)
{
if(p&1)ans=(ans*w)%mod;
w=(w*w)%mod;
}
return ans;
}
int main()
{
scanf("%lld",&n);
if(n==1||n==2)printf("%d\n",1);
else
{
matrix.init();
matrix=pow(matrix,n-2);
printf("%lld\n",(matrix[1][1]+matrix[1][2])%mod);
}
return 0;
}
|
例2
（洛谷P1939）a[n]=a[n-3]+a[n-1]，a[1]=a[2]=a[3]=1。
初看时没有什么头绪，但是对递推公式整理可以得到：
a[x]=a[x-3]+a[x-1]=2a[x-3]+a[x-4]+a[x-5]
a[x-1]=a[x-4]+a[x-2]=a[x-3]+a[x-4]+a[x-5]
a[x-2]=a[x-3]+a[x-5]
这样a[x]~a[x-2]都可以用a[x-3]~a[x-5]表示，所以构造的矩阵为 $\begin{bmatrix}
2 & 1 & 1\\
1 & 1 & 1\\
1 & 0 & 1
\end{bmatrix}$。
但是还有一个问题，因为每次递推是以3为步长的，
所以推到最后等式右边不一定正好是 $\begin{bmatrix}
f[3] \\
f[2] \\
f[1]
\end{bmatrix}$。
那么根据x mod 3还要分成以下几类：
x mod 3=0： $\begin{bmatrix}
f[ 3]\\
f[ 2]\\
f[ 1]
\end{bmatrix} =\begin{bmatrix}
1\\
1\\
1
\end{bmatrix}$
x mod 3=1： $\begin{bmatrix}
f[ 4]\\
f[ 3]\\
f[ 2]
\end{bmatrix} =\begin{bmatrix}
2\\
1\\
1
\end{bmatrix}$
x mod 3=2： $\begin{bmatrix}
f[ 5]\\
f[ 4]\\
f[ 3]
\end{bmatrix} =\begin{bmatrix}
3\\
2\\
1
\end{bmatrix}$
手工计算得到f[4]=2,f[5]=3，于是我们可以将这三类递推到最后等式右边的值记录下来，存储为：
int a[3][4]={ {0,1,1,1}, {0,2,1,1}, {0,3,2,1}};
还有一个问题：这样一来矩阵的幂指数是不确定的。
我们观察到，对于所要求的项a[x]，
矩阵的指数k应当满足如下条件：
3<=x-3k<=5，且k为整数。
解得：(x-5)/3 <=k<= (x-3)/3
所以在程序中，可以用
k=ceil(((double)x-5.0)/3.0) 或 k=floor(((double)x-3.0)/3.0)
来得到k的值。
代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
| /*
luogu1939
a[x]=a[x-3]+a[x-1]=2a[x-3]+a[x-4]+a[x-5]
a[x-1]=a[x-4]+a[x-2]=a[x-3]+a[x-4]+a[x-5]
a[x-2]=a[x-3]+a[x-5]
the matrix is:
|2,1,1|
|1,1,1|
|1,0,1|
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
using namespace std;
typedef long long lint;
lint n,T;
const lint mod=lint(1e9+7);
struct MATRIX
{
lint ary[4][4];
void clear()
{
for(int i=0;i<=3;++i)
for(int j=0;j<=3;++j)
ary[i][j]=0;
}
MATRIX()
{
clear();
}
void init()
{
for(int i=1;i<=3;++i)
{
for(int j=1;j<=3;++j)
{
ary[i][j]=1;
if(i==1&&j==1)ary[i][j]=2;
if(i==3&&j==2)ary[i][j]=0;
}
}
}
lint* operator[](int pos)
{
return ary[pos];
}
void sete()
{
for(int i=1;i<=3;++i)
{
for(int j=1;j<=3;++j)
{
ary[i][j]=lint(i==j);
}
}
}
}matrix;
MATRIX operator*(MATRIX a,MATRIX b)
{
MATRIX ans;
for(int i=1;i<=3;++i)
for(int j=1;j<=3;++j)
for(int k=1;k<=3;++k)
ans[i][j]+=a[i][k]*b[k][j],ans[i][j]%=mod;
return ans;
}
MATRIX operator%(MATRIX a,lint mo)
{
for(int i=1;i<=3;++i)
for(int j=1;j<=3;++j)
a[i][j]%=mo;
return a;
}
MATRIX pow(MATRIX w,lint p)
{
MATRIX ans;
for(ans.sete();p;p/=2)
{
if(p&1)ans=(ans*w)%mod;
w=(w*w)%mod;
}
return ans;
}
int a[3][4]={{0,1,1,1},{0,2,1,1},{0,3,2,1}};
int main()
{
//freopen("luogu1939_1.in","r",stdin);
//freopen("luogu1939.out","w",stdout);
for(scanf("%lld",&T);T;--T)
{
scanf("%lld",&n);
if(n==1||n==2||n==3)
printf("%d\n",1);
else if(n==4)
printf("%d\n",2);
else if(n==5)
printf("%d\n",3);
else
{
lint k=lint(ceil((double(n-5.0)/3.0)));
matrix.init();
matrix=pow(matrix,k);
lint tans=
(matrix[1][1]*a[n%3][1]
+matrix[1][2]*a[n%3][2]
+matrix[1][3]*a[n%3][3])%mod;
printf("%lld\n",tans);
}
}
return 0;
}
|
合并点的坐标变换
对于一些平面直角坐标系中的坐标的变换，可以构造出一个特定的矩阵成与矩阵$\begin{bmatrix}
x\\
y
\end{bmatrix}$相乘。而利用矩阵乘法，可以将一系列坐标变换变成一个。
常见点的变换
平移
给定点P(x,y)，平面向量$\vec{d} =( p,q) ,P'=P+\vec{d} ,P'( x+p,y+q)$。设构造的矩阵[A]，
则有：$[ A] \times \begin{bmatrix}x\\y\end{bmatrix} =\begin{bmatrix}x+p\\y+q\end{bmatrix}$ 。
设[A]大小为n*m，根据矩阵乘法的条件，A(n*m)×B(2*1)=C(2*1)，可以得到n=m=2。
设$[ A]=\begin{bmatrix}a&b\\c&d\end{bmatrix} $，则有：
ax+by=x+p,cx+dy=y+q。
发现不可能找到一组常数a,b,c,d使得上式成立，即不存在这样的[A]。
So？
我们将坐标矩阵加一行，变成 $\begin{bmatrix}x\\y\\1\end{bmatrix}$ ，那么相应的 $[A]=\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix}$ 。有：
$$
[A]\times\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}x+p\\y+p\\1\end{bmatrix}
$$
注意：[A]应该写在前面，否则不符合矩阵乘法的条件。
同理得到：
$$
\begin{cases}ax+by+c=x+p\\dx+ey+f=y+q\\gx+hy+i=1\end{cases}
$$
解得：$[A]=\begin{bmatrix}1&0&p\\0&1&q\\0&0&1\end{bmatrix}$。
类似地，可以得到其他几个坐标变换对应的构造矩阵：
关于x轴对称
P’(x,-y)，则有： $\begin{bmatrix}1&0&0\\0&-1&0\\0&0&1\end{bmatrix}\times\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}x\\-y\\1\end{bmatrix}$
关于y轴对称
P’(-x,y)，则有： $\begin{bmatrix}-1&0&0\\0&1&0\\0&0&1\end{bmatrix}\times\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}-x\\y\\1\end{bmatrix}$
旋转
https://jingyan.baidu.com/article/2c8c281dfbf3dd0009252a7b.html
根据数学知识，得到P(x,y)绕O(a,b)旋转 （弧度制）得到点P’(x’,y’)，则有：
$$
\begin{cases}x'=(x-a)cosθ-(y-b)sinθ+a\\y'=(x-a)sinθ+(y-b)cosθ+b\end{cases}
$$
则有：
$$
\begin{bmatrix}cosθ&-sinθ&-acosθ+bsinθ+a\\sinθ&cosθ&-asinθ-bcosθ+b\\0&0&1\end{bmatrix}\times\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}(x-a)cosθ-(y-b)sinθ+a\\(x-a)sinθ+(y-b)cosθ+b\\1\end{bmatrix}
$$
特殊地，当O为坐标轴圆点时，有：
$$
\begin{bmatrix}cosθ&-sinθ&0\\sinθ&cosθ&0\\0&0&1\end{bmatrix}\times\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}xcosθ-ysinθ\\xsinθ+ycosθ\\1\end{bmatrix}
$$
坐标放大/缩小
即将坐标乘上一个数k。 $\begin{bmatrix}k&0&0\\0&k&0\\0&0&1\end{bmatrix}×\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}kx\\ky\\1\end{bmatrix}$
要点
- 在将坐标矩阵[P](3*1)与构造的矩阵[A](3*3)相乘时，应该写成[P1]=[A1]*[P0]，否则无法相乘
2.[Pn]=[An]*[P(n-1)]=[An]*( [a(n-1)]*[p(n-2)] )
发现构造的矩阵A始终是写在前面，而且是按An*A(n-1)*A(n-2)…的逆序形式相乘的，
又由于矩阵乘法不满足交换律，所以不能顺序相乘。
有两种解决办法：一是计算出所有Ai后逆序累乘，二则可以写成proA=An*proA，即将An写在乘号的左边。
模板，解析略。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
| /*
nyoj298
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<cctype>
using namespace std;
const double eps=1e-9;
const double pi=acos(-1.0);
const int maxn=10003;
const int maxm=1000003;
int n,m;
template<typename T,int nn,int mm>
struct MATRIX
{
double ary[nn+1][mm+1];
T* operator[](int pos)
{
return ary[pos];
}
void clear()
{
for(int i=0;i<=nn;++i)
{
for(int j=0;j<=mm;++j)
{
ary[i][j]=0.0;
}
}
}
MATRIX(){clear();}
void sete()
{
clear();
for(int i=1;i<=3;++i)
ary[i][i]=1;
}
void initMove(double p,double q)
{
clear();
for(int i=1;i<=3;++i)
ary[i][i]=1.0;
ary[1][3]=p,ary[2][3]=q;
}
void initX()
{
clear();
for(int i=1;i<=3;++i)
{
if(i&1)ary[i][i]=1;
else ary[i][i]=-1;
}
}
void initY()
{
clear();
for(int i=1;i<=3;++i)
{
if(i!=1)ary[i][i]=1;
else ary[i][i]=-1;
}
}
void initZoom(double k)
{
clear();
for(int i=1;i<=3;++i)
{
if(i!=3)ary[i][i]=k;
else ary[i][i]=1;
}
}
void initRotate(double theta)
{
clear();
theta=theta/180.0*pi;
double cost=cos(theta);
double sint=sin(theta);
if(cost<eps&&cost>-eps)cost=0.0;
if(sint<eps&&sint>-eps)sint=0.0;
for(int i=1;i<=3;++i)
{
if(i!=3)ary[i][i]=cost;
else ary[i][i]=1;
}
ary[1][2]=-sint,ary[2][1]=sint;
}
};
template<typename T,int nn,int rr,int mm>
MATRIX<T,nn,mm> operator*(MATRIX<T,nn,rr> a,MATRIX<T,rr,mm> b)
{
MATRIX<T,nn,mm> ans;
for(int i=1;i<=nn;++i)
for(int j=1;j<=mm;++j)
for(int k=1;k<=rr;++k)
ans[i][j]+=a[i][k]*b[k][j];
return ans;
}
MATRIX<double,3,1> p[maxn];
void input()
{
scanf("%d%d",&n,&m);
double x,y;
for(int i=1;i<=n;++i)
{
scanf("%lf%lf",&x,&y);
p[i][1][1]=x;
p[i][2][1]=y;
p[i][3][1]=1;
}
}
MATRIX<double,3,3> trans,tmp;
void solve()
{
/*
首字符如果是M,则表示平移操作，该行后面将跟两个数x,y，表示把所有点按向量(x,y)平移;
首字符如果是X，则表示把所有点相对于X轴进行上下翻转;
首字符如果是Y，则表示把所有点相对于Y轴进行左右翻转;
首字符如果是S，则随后将跟一个数P,表示坐标放大P倍;
首字符如果是R，则随后将跟一个数A,表示所有点相对坐标原点逆时针旋转一定的角度A(单位是度)
*/
trans.sete();
for(int i=1;i<=m;++i)
{
char r;
do{r=getchar();}
while(!isupper(r));
if(r=='M')
{
double x,y;
scanf("%lf%lf",&x,&y);
tmp.initMove(x,y);
}
else if(r=='X')
{
tmp.initX();
}
else if(r=='Y')
{
tmp.initY();
}
else if(r=='S')
{
double k;
scanf("%lf",&k);
tmp.initZoom(k);
}
else if(r=='R')
{
double t;
scanf("%lf",&t);
tmp.initRotate(t);
}
trans=tmp*trans;
}
}
void output()
{
for(int i=1;i<=n;++i)
{
p[i]=trans*p[i];
printf("%.1lf %.1lf\n",p[i][1][1],p[i][2][1]);
}
}
int main()
{
freopen("nyoj298_1.in","r",stdin);
// freopen("nyoj298.out","w",stdout);
input();
solve();
output();
return 0;
}
|
合并线性数组的部分操作
假设有线性数组a[]，将a构造为形为$\begin{bmatrix}a1\\a2\\...\\an\end{bmatrix}$的矩阵，然后对于a[]的部分操作就可以构造成一个矩阵与上述矩阵相乘。
对于重复的大量操作，可以合并成一个矩阵来加快速度。
因为每次操作一般只改变a中很少的几个元素，所以各种操作的构造矩阵基本都是以单位矩阵为模板。
构造矩阵
说明：以下内容中，E表示单位矩阵，且构造的矩阵应放在乘号的左边。
1
2
3
4
| A[i]*=k：
E[i][i]=k
A[i]+=A[j]：
E[i][j]=1
|
推广：A[i]+=k，可以通过两种方法来实现：（未经验证）
一是扩大矩阵的大小,将A[n+1]设为k，然后按上述方法加法后再截取结果矩阵大小n的子矩阵；二是令t=A[0]，A[0]=k，按上述方法相加后再令A[0]=t。
调整元素位置
按希望得到的排列调整单位矩阵的各行。
特殊地：Swap(A[i],A[j])：E[i][i]=E[j][j]=0,E[i][j]=E[j][i]=1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
| /*
caioj1483
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<cctype>
#include<iomanip>
using namespace std;
typedef long long lint;
template<typename T>
void read(T& w)
{
char r;int f=1;
for(r=getchar();r!='-'&&(r<48||r>57);r=getchar());
if(r=='-')r=getchar(),f=-1;
for(w=0;r>=48&&r<=57;r=getchar())w=w*10+r-48;
w*=f;
}
const int maxn=103;
const int maxm=23;
lint n,m,t,s;
template<typename T>
struct MATRIX
{
T ary[maxn][maxn];
T* operator[](int pos)
{
return ary[pos];
}
void clear()
{
for(int i=0;i<=n;++i)
{
for(int j=0;j<=n;++j)
{
ary[i][j]=0;
}
}
}
void sete()
{
clear();
for(int i=1;i<=n;++i)
ary[i][i]=1;
}
void initA()
{
clear();
for(int i=1;i<=n;++i)
{
ary[i][1]=1;
}
}
void initD(lint i)
{
sete();
ary[i][i]=0;
}
void initR(lint i,lint k)
{
sete();
ary[i][i]=k;
}
void initC(lint i,lint j)
{
sete();
ary[i][j]=1;
}
void initT(lint i,lint j)
{
sete();
ary[i][j]=1;
ary[j][j]=0;
}
void initS(lint i,lint j)
{
sete();
ary[i][i]=ary[j][j]=0;
ary[j][i]=ary[i][j]=1;
}
void initM()
{
clear();
for(int i=2;i<=n;++i)
{
ary[i][i-1]=1;
}
ary[1][n]=1;
}
void out()
{
for(int i=1;i<=n;++i)
{
for(int j=1;j<=n;++j)
cout<<setw(4)<<ary[i][j];
cout<<endl;
}
cout<<endl;
}
};
template<typename T>
MATRIX<T> operator*(MATRIX<T> a,MATRIX<T> b)
{
MATRIX<T> ans;
ans.clear();
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
for(int k=1;k<=n;++k)
ans[i][j]+=a[i][k]*b[k][j]%s,ans[i][j]%=s;
return ans;
}
MATRIX<lint> box,op[maxm],pdt[maxm];
void input()
{
read(n),read(m),read(t),read(s);
box.initA(),pdt[0].sete();
lint x,y;
for(int i=1;i<=m;++i)
{
char r;
do{r=getchar();}
while(!islower(r));
read(x),read(y);
if(r=='d')
op[i].initD(x);
else if(r=='r')
op[i].initR(x,y);
else if(r=='c')
op[i].initC(x,y);
else if(r=='t')
op[i].initT(x,y);
else if(r=='s')
op[i].initS(x,y);
else if(r=='m')
op[i].initM();
pdt[i]=op[i]*pdt[i-1];
}
}
template<typename T>
MATRIX<T> fpow(MATRIX<T> w,lint p)
{
MATRIX<T> ans;
for(ans.sete();p;p>>=1)
{
if(p&1)
ans=ans*w;
w=w*w;
}
return ans;
}
MATRIX<lint> ans;
void solve()
{
ans=pdt[t%m]*fpow(pdt[m],t/m)*box;
}
void output()
{
for(int i=1;i<=n;++i)
{
printf("%lld\n",ans[i][1]);
}
}
int main()
{
freopen("caioj1483_1.in","r",stdin);
//freopen("caioj1483.out","w",stdout);
input();
solve();
output();
return 0;
}
|
练习题
