抱歉,您的浏览器无法访问本站
本页面需要浏览器支持(启用)JavaScript
了解详情 >

matrix.h

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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
#include <cassert>
#include <iostream>
#include <vector>
#include <cstring>
//template<typename Num>
#define Num double
struct Matrix{
std::string Message;
std::vector<std::vector<Num> >M;
int row,col;//行数、列数
Matrix(int rows,int cols){
//初始化为零矩阵
// std::cout<<rows<<" "<<cols;
row=rows,col=cols;
M.resize(row+1);
for (int i=1;i<=row;++i){
M[i].resize(col+1,0);
}
Message="OK";
// std::cout<<"OK"<<std::endl;
}
void init(int rows,int cols){
row=rows,col=cols;
M.resize(row+1);
for (int i=1;i<=row;++i){
M[i].resize(col+1,0);
}
Message="OK";
}
Matrix(){
row=col=0;
}
void Input(){
std::cout<<"Please enter the row and col of the matrix"<<std::endl;
std::cin>>row>>col;
assert(row>=1 && col >=1);
std::cout<<"Please enter the matrix"<<std::endl;
init(row,col);
for (int i=1;i<=row;++i){
for (int j=1;j<=col;++j){
std::cin>>M[i][j];
}
}
}
void Output(){
std::cout<<"--------------"<<std::endl;
for (int i=1;i<=row;++i){
for (int j=1;j<=col;++j){
std::cout<<M[i][j]<<" ";
}
std::cout<<std::endl;
}
}
std::vector<Num>&operator[](int i){
return M[i];
}
void identityMatrix(){
assert(row==col);
for (int i=1;i<=row;++i){
M[i][i]=1;
}
}
void diag(std::vector<Num>V){
Matrix((int)V.size(),(int)V.size());
for (int i=0;i<(int)V.size();++i){
M[i+1][i+1]=V[i];
}
}
//矩阵的初等行、列变换,type='R'行,type='C'列
void swap(char type,int x,int y){
if (type=='R'){
assert(1<=x && x<=row && 1<=y && y<=row);
for (int i=1;i<=col;++i){
std::swap(M[x][i],M[y][i]);
}
}
else if (type=='C'){
assert(1<=x && x<=col && 1<=y && y<=col);
for (int i=1;i<=row;++i){
std::swap(M[i][x],M[i][y]);
}
}
else {
assert(0);
}
}
void times(char type,int x,Num t){
if (type=='R'){
assert(1<=x && x<=row);
for (int i=1;i<=col;++i){
M[x][i]=t*M[x][i];
}
}
else if (type=='C'){
assert(1<=x && x<=col);
for (int i=1;i<=row;++i){
M[i][x]=t*M[i][x];
}
}
else {
assert(0);
}
}
void addtimes(char type,int x,Num t,int y){//x行/列的t倍加到y上
if (type=='R'){
assert(1<=x && x<=row && 1<=y && y<=row);
for (int i=1;i<=col;++i){
M[y][i]=M[y][i]+t*M[x][i];
}
}
else if (type=='C'){
assert(1<=x && x<=col && 1<=y && y<=col);
for (int i=1;i<=row;++i){
M[i][y]=M[i][y]+t*M[i][x];
}
}
else {
assert(0);
}
}
Num sum(char type,int x){
if (type=='R'){
Num s=0;
for (int i=1;i<=col;++i){
s=s+M[x][i];
}
return s;
}
else if (type=='C'){
Num s=0;
for (int i=1;i<=row;++i){
s=s+M[i][x];
}
return s;
}
else {
assert(0);
}
return 0;
}
Matrix transpose(){//矩阵的转置
Matrix B(col,row);
for (int i=1;i<=B.row;++i){
for (int j=1;j<=B.col;++j){
B[i][j]=M[j][i];
}
}
return B;
}
};
Matrix inverse(Matrix A){//求逆矩阵
assert(A.row==A.col);
int n=A.row;
Matrix B(A.row,A.col);
B.identityMatrix();
for (int i=1;i<=n;++i){
int r=i;
for (int j=i+1;j<=n;++j){
if (abs(A[j][i])>abs(A[r][i])) r=j;
}
A.swap('R',i,r),B.swap('R',i,r);
if (A[i][i]==0){
B.Message="Not full rank!";
return B;
}
Num inv=((Num)(1))/A[i][i];
for (int j=1;j<=n;++j){
if (j==i) continue;
Num t=A[j][i]*inv;
A.addtimes('R',i,t,j);
B.addtimes('R',i,t,j);
}
}
for (int i=1;i<=n;++i){
Num inv=((Num)(1))/A[i][i];
A.times('R',i,inv);
B.times('R',i,inv);
}
B.Message="OK";
return B;
}
Matrix operator * (Matrix A,Matrix B){
assert(A.col==B.row);
Matrix C(A.row,B.col);
for (int i=1;i<=A.row;++i){
for (int j=1;j<=B.col;++j){
for (int k=1;k<=A.col;++k){
C[i][j]+=A[i][k]*B[k][j];
}
}
}
return C;
}
Matrix operator * (Num lambda,Matrix A){
Matrix B(A.row,A.col);
for (int i=1;i<=A.row;++i){
for (int j=1;j<=B.col;++j){
B[i][j]=lambda*A[i][j];
}
}
return B;
}
Matrix operator + (Matrix A,Matrix B){
assert(A.row==B.row && A.col==B.col);
Matrix C(A.row,A.col);
for (int i=1;i<=A.row;++i){
for (int j=1;j<=A.col;++j){
C[i][j]=A[i][j]+B[i][j];
}
}
return C;
}
Matrix operator / (Num x,Matrix A){
return x*inverse(A);
}
Matrix operator / (Matrix A,Matrix B){
return A*inverse(B);
}
Matrix operator ^ (Matrix A,int k){
assert(A.row==A.col);
k--;
Matrix base=A;
while (k){
if (k&1) base=base*A;
A=A*A;
k>>=1;
}
return base;
}
Num tr(Matrix A){
assert(A.row==A.col);
Num ans=0;
for (int i=1;i<=A.row;++i){
ans=ans+A[i][i];
}
return ans;
}
Matrix addHorizontal(Matrix A,Matrix B){
assert(A.row==B.row);
Matrix C(A.row,A.col+B.col);
for (int i=1;i<=A.row;++i){
for (int j=1;j<=A.col;++j){
C[i][j]=A[i][j];
}
for (int j=1;j<=B.col;++j){
C[i][A.col+j]=B[i][j];
}
}
return C;
}
Matrix addVertical(Matrix A,Matrix B){
assert(A.col==B.col);
Matrix C(A.row+B.row,A.col);
for (int i=1;i<=A.col;++i){
for (int j=1;j<=A.row;++j){
C[j][i]=A[j][i];
}
for (int j=1;j<=B.row;++j){
C[A.row+j][i]=B[j][i];
}
}
return C;
}

T1

间接胜就是矩阵乘法。如果矩阵是一次操作,代表将原先的队伍(单位矩阵)转换为其胜的队伍,那么两次操作就代表把原先的队伍转换为间接胜的队伍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include "matrix.h"
#include <algorithm>
int main(){
Matrix A(5,5);
A[1][2]=1,A[1][3]=1;
A[2][3]=1,A[2][4]=1,A[2][5]=1;
A[4][1]=1,A[4][3]=1,A[4][5]=1;
A[5][1]=1,A[5][3]=1;
Matrix B=A*A;
A.Output(),(A*A).Output();
std::vector<std::pair<int,std::pair<int,int> > >rk;
for (int i=1;i<=5;++i){
rk.push_back(std::make_pair(-A.sum('R',i),std::make_pair(-B.sum('R',i),i)));
}
std::sort(rk.begin(),rk.end());
for (int i=0;i<(int)rk.size();++i){
std::cout<<rk[i].second.second<<" ";
}
}

T2

经典的模型,简单的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include "matrix.h"
#include <algorithm>
int main(){
Matrix A(7,6);
A.Input();//实在不想手敲了。
/*
1 1 0 1 0 0
1 0 1 1 0 0
1 1 1 1 0 0
1 1 0 1 0 0
1 1 0 1 0 0
1 0 0 0 1 0
0 0 0 0 1 1
*/
Matrix alpha(6,1);
alpha[1][1]=1,alpha[4][1]=1,alpha[2][1]=1;
for (int i=1;i<=7;++i){
if (alpha.sum('C',1)==(A*alpha).sum('R',i)){
std::cout<<i<<" matches"<<std::endl;
}
}
}

T3,4,5

简单解方程。

T6

见:https://gaisaiyuno.github.io/archives/f2c32ff8.html

T8

用矩阵操作,将 α=(a1,b1,c1),β=(a2,b2,c2)\alpha=(a_1,b_1,c_1),\beta=(a_2,b_2,c_2) 转换为子代的基因型,考虑实二次型 f=xTAxf=x^TAx,这里,将 xTx^Txx 分别改成 αT\alpha^Tβ\beta 即可。

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
#include "matrix.h"
#include <algorithm>
Matrix A1(3,3),A2(3,3),A3(3,3);
Matrix solve(Matrix alpha,Matrix beta){
Matrix ret(3,1);
ret[1][1]=(alpha.transpose()*A1*beta)[1][1];
ret[2][1]=(alpha.transpose()*A2*beta)[1][1];
ret[3][1]=(alpha.transpose()*A3*beta)[1][1];
return ret;
}
int main(){
A1[1][1]=1,A1[2][2]=0.25,A1[1][2]=A1[2][1]=0.5;
A2[1][2]=A2[2][3]=A2[2][1]=A2[3][2]=0.5,A2[1][3]=A2[3][1]=1,A2[2][2]=0.5;
A3[3][3]=1,A3[2][2]=0.25,A3[2][3]=A3[3][2]=0.5;
//Plan 1
Matrix start(3,1);
start[1][1]=start[2][1]=start[3][1]=1.0/3.0;
for (int i=1;i<=10;++i){
Matrix alpha(3,1);
alpha[1][1]=1;
Matrix beta(3,1);
beta=start;
start=solve(alpha,beta);
start.Output();
}
//plan 2
start[1][1]=start[2][1]=start[3][1]=1.0/3.0;
for (int i=1;i<=10;++i){
Matrix alpha(3,1);
alpha[2][1]=1;
Matrix beta(3,1);
beta=start;
start=solve(alpha,beta);
start.Output();
}
//plan 3
start[1][1]=start[2][1]=start[3][1]=1.0/3.0;
for (int i=1;i<=10;++i){
Matrix a(3,1),b(3,1),c(3,1);
a[1][1]=start[1][1],b[2][1]=start[2][1],c[3][1]=start[3][1];
start=solve(a,a)+solve(b,b)+solve(c,c);
start=(1.00/start.sum('C',1))*start;
start.Output();
}
}

T10

1
2
3
4
5
6
7
8
9
10
11
12
#include "matrix.h"
#include <algorithm>
int main(){
Matrix A(2,2);
A[1][1]=1,A[1][2]=1,A[2][1]=1;
Matrix Base(2,2);
Base.identityMatrix();
for (int i=1;i<=10;++i){
Base=Base*A;
Base.Output();
}
}

评论