矩阵是干什么的呢?一句话来说就是,知道相邻两个函数的递推关系和第一个数,让你递推到第n个数。显然,如果n很大,那么一个一个递推过去是会超时的。所以矩阵就是用来解决这种快速递推的问题的。
比方说斐波那契数列就是一个递推的典型。
先丢题目链接:我是题目!
那么问题的关键就变成了如何找递推关系的中介矩阵temp了。如果题目告诉了你递推关系,题目就变得很简单了。但是告诉你递推关系大致可分为两类,一类是加法的递推,像斐波那契数列,temp矩阵的每个元素都是常数。另外一种,就是有乘法的递推关系了,这类关系一般需要凑出来,有时候隐晦的话也不是那么容易的。比如说上面题目中的I题,就是需要拼凑的。
总之先交代模板,先看H题。这题求矩阵sum(A)=A^1+A^2+A^3+...+A^n。方法的话是二分,左边这个式子如果n是偶数的话可以拆成(1+A^(n/2))*(A^1+A^2+...+A(n/2)),然后右边部分又可以递归,不断二分地递归即可。当然,如果n是奇数要再加上A^n;另外上面式子中的1在这里指的是同阶的单位矩阵。代码如下(也可以直接当做矩阵的模板了):
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <stdio.h> 2 #include <algorithm> 3 #include <math.h> 4 #include <vector> 5 #include <map> 6 #include <set> 7 #include <iostream> 8 #include <string.h> 9 #pragma comment(linker, "/STACK:1024000000,1024000000") 10 using namespace std; 11 typedef long long ll; 12 typedef pair<int,int> pii; 13 const int mod = 10; 14 15 //int mod; 16 17 void add(int &a,int b) 18 { 19 a += b; 20 if(a < 0) a += mod; 21 //if(a >= mod) a -= mod; 22 a %= mod; 23 } 24 25 struct matrix 26 { 27 int e[50+5][50+5],n,m; 28 matrix() {} 29 matrix(int _n,int _m): n(_n),m(_m) {memset(e,0,sizeof(e));} 30 matrix operator * (const matrix &temp)const 31 { 32 matrix ret = matrix(n,temp.m); 33 for(int i=1;i<=ret.n;i++) 34 { 35 for(int j=1;j<=ret.m;j++) 36 { 37 for(int k=1;k<=m;k++) 38 { 39 add(ret.e[i][j],1LL*e[i][k]*temp.e[k][j]%mod); 40 } 41 } 42 } 43 return ret; 44 } 45 matrix operator + (const matrix &temp)const 46 { 47 matrix ret = matrix(n,m); 48 for(int i=1;i<=n;i++) 49 { 50 for(int j=1;j<=m;j++) 51 { 52 add(ret.e[i][j],(e[i][j]+temp.e[i][j])%mod); 53 } 54 } 55 return ret; 56 } 57 void getE() 58 { 59 for(int i=1;i<=n;i++) 60 { 61 for(int j=1;j<=m;j++) 62 { 63 e[i][j] = i==j?1:0; 64 } 65 } 66 } 67 }; 68 69 matrix qpow(matrix temp,int x) 70 { 71 int sz = temp.n; 72 matrix base = matrix(sz,sz); 73 base.getE(); 74 while(x) 75 { 76 if(x & 1) base = base * temp; 77 x >>= 1; 78 temp = temp * temp; 79 } 80 return base; 81 } 82 83 void print(matrix p) 84 { 85 int n = p.n; 86 int m = p.m; 87 for(int i=1;i<=n;i++) 88 { 89 for(int j=1;j<=m;j++) 90 { 91 printf("%d ",p.e[i][j]); 92 } 93 cout << endl; 94 } 95 } 96 97 matrix solve(matrix a, int k) 98 { 99 if(k == 1) return a; 100 int n = a.n; 101 matrix temp = matrix(n,n); 102 temp.getE(); 103 if(k & 1) 104 { 105 matrix ex = qpow(a,k); 106 k--; 107 temp = temp + qpow(a,k/2); 108 return temp * solve(a,k/2) + ex; 109 } 110 else 111 { 112 temp = temp + qpow(a,k/2); 113 return temp * solve(a,k/2); 114 } 115 } 116 117 int main() 118 { 119 int n,k; 120 while(scanf("%d%d",&n,&k)==2 && n) 121 { 122 matrix a = matrix(n,n); 123 for(int i=1;i<=n;i++) 124 { 125 for(int j=1;j<=n;j++) {scanf("%d",&a.e[i][j]);a.e[i][j] %= mod;} 126 // 这里矩阵的每个元素在输入以后就要mod一下,不然可能会因为输入的元素太大导致在后面爆int(?) 127 } 128 matrix ans = solve(a,k); 129 for(int i=1;i<=n;i++) 130 { 131 for(int j=1;j<=n;j++) 132 { 133 printf("%d%c",ans.e[i][j],j==n?'\n':' '); 134 } 135 } 136 puts(""); 137 } 138 }
还要讲的几点是:1.E题中的类型,A是1000*6的矩阵,B是6*1000的矩阵,现在求(A*B)^N,显然1000*1000的矩阵是要超时的,那么不妨化成A*(B*A)^(n-1)*B,这样B*A是6*6的矩阵,就快多了;2.矩阵不要无脑开的很大,因为我发现有时候矩阵开大了时间会增加很多;3.关于循环矩阵的问题,比方说一个矩阵,它的每一行都是和第一行一样的,只是位置每行错开一个位置(当然列也是一样),这样的矩阵就叫做循环矩阵,在计算的过程中只要保存一行即可(即使是一行也可以知道整个矩阵的内容了),这样的题目出现在J题。这题的代码如下:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <stdio.h> 2 #include <algorithm> 3 #include <math.h> 4 #include <vector> 5 #include <map> 6 #include <set> 7 #include <iostream> 8 #include <string.h> 9 #pragma comment(linker, "/STACK:1024000000,1024000000") 10 using namespace std; 11 typedef pair<int,int> pii; 12 //const int mod = 1000; 13 14 int n,mod,d,k; 15 16 void add(int &a,int b) 17 { 18 a += b; 19 if(a < 0) a += mod; 20 //if(a >= mod) a -= mod; 21 a %= mod; 22 } 23 24 struct matrix 25 { 26 int e[5][500+5],n,m; 27 matrix() {} 28 matrix(int _n,int _m): n(_n),m(_m) {memset(e,0,sizeof(e));} 29 matrix operator * (const matrix &temp)const 30 { 31 matrix ret = matrix(1,m); 32 for(int i=0;i<m;i++) 33 { 34 for(int j=0;j<m;j++) 35 { 36 add(ret.e[0][i],1LL*e[0][j]*temp.e[0][((j-i)%m+m)%m]%mod); 37 } 38 } 39 return ret; 40 } 41 42 }; 43 44 matrix qpow(matrix temp,int x) 45 { 46 int sz = temp.m; 47 matrix base = matrix(1,sz); 48 base.e[0][0] = 1; 49 while(x) 50 { 51 if(x & 1) base = base * temp; 52 x >>= 1; 53 temp = temp * temp; 54 } 55 return base; 56 } 57 58 void print(matrix p) 59 { 60 int n = p.n; 61 int m = p.m; 62 for(int i=0;i<n;i++) 63 { 64 for(int j=0;j<m;j++) 65 { 66 printf("%d ",p.e[i][j]); 67 } 68 cout << endl; 69 } 70 } 71 72 int main() 73 { 74 while(scanf("%d%d%d%d",&n,&mod,&d,&k)==4) 75 { 76 matrix ans = matrix(1,n); 77 for(int i=0;i<n;i++) scanf("%d",&ans.e[0][i]); 78 matrix temp = matrix(1,n); 79 //for(int i=0;i<n;i++) 80 { 81 temp.e[0][0]=1; 82 for(int j=1;j<=d;j++) 83 { 84 int now = 0+j; 85 if(now >= n) now -= n; 86 temp.e[0][now] = 1; 87 } 88 for(int j=1;j<=d;j++) 89 { 90 int now = 0-j; 91 if(now < 0) now += n; 92 temp.e[0][now] = 1; 93 } 94 } 95 //print(temp); 96 temp = qpow(temp,k); 97 ans = ans * temp; 98 for(int i=0;i<n;i++) 99 { 100 printf("%d%c",ans.e[0][i],i==n-1?'\n':' '); 101 } 102 } 103 }
要注意之所以这份代码里面的qpow函数的单位矩阵可以直接相乘,是因为单位矩阵也是一个循环矩阵,且循环性质和题中所给的矩阵一样。
另外,这份题目中其他题目的好题如下:B,C,M,N,P,R。
以上几题想稍微补充一下的是,M题,ceil里面的数加上它的共轭函数刚好是一个整数,证明如下:展开以后,有根号b的部分,如果其次数是偶数,那么一定是整数,如果是奇数,那么共轭的两个这个部分刚好一正一负,相抵消。这点证明了以后,又因为(a-1) 2< b < a2,因此(a-根号b)在0到1之间,所以共轭的两部分相加刚好是前面部分取ceil的值,即所求的没有%m的部分。对R题,虽然题目不算难,但是要注意的是,A^B%C,不能直接取余,要结合数论中的知识来,如下:
1 /* 2 A^B %C 这题的C是质数,而且A,C是互质的。 3 所以直接A^(B%(C-1)) %C 4 5 比较一般的结论是 A^B %C=A^( B%phi(C)+phi(C) ) %C , B>=phi(C) 6 */
最后想说的是,D题和L题并没有做,而且P题最后部分也没法很好的证明,P题还是留个代码吧:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <stdio.h> 2 #include <algorithm> 3 #include <math.h> 4 #include <vector> 5 #include <map> 6 #include <set> 7 #include <iostream> 8 #include <string.h> 9 using namespace std; 10 typedef long long ll; 11 typedef pair<int,int> pii; 12 //const int mod = 1e9 + 7; 13 14 int mod; 15 16 void add(int &a,int b) 17 { 18 a += b; 19 if(a < 0) a += mod; 20 //if(a >= mod) a -= mod; 21 a %= mod; 22 } 23 24 struct matrix 25 { 26 int e[10][10]; 27 int n,m; 28 matrix() {} 29 matrix(int _n,int _m): n(_n),m(_m) {memset(e,0,sizeof(e));} 30 matrix operator * (const matrix &temp)const 31 { 32 matrix ret = matrix(n,temp.m); 33 for(int i=1;i<=ret.n;i++) 34 { 35 for(int j=1;j<=ret.m;j++) 36 { 37 for(int k=1;k<=m;k++) 38 { 39 add(ret.e[i][j],1LL*e[i][k]*temp.e[k][j]%mod); 40 } 41 } 42 } 43 return ret; 44 } 45 matrix operator + (const matrix &temp)const 46 { 47 matrix ret = matrix(n,m); 48 for(int i=1;i<=n;i++) 49 { 50 for(int j=1;j<=m;j++) 51 { 52 add(ret.e[i][j],(e[i][j]+temp.e[i][j])%mod); 53 } 54 } 55 return ret; 56 } 57 void getE() 58 { 59 for(int i=1;i<=n;i++) 60 { 61 for(int j=1;j<=m;j++) 62 { 63 e[i][j] = i==j?1:0; 64 } 65 } 66 } 67 }; 68 69 matrix qpow(matrix temp,int x) 70 { 71 int sz = temp.n; 72 matrix base = matrix(sz,sz); 73 base.getE(); 74 while(x) 75 { 76 if(x & 1) base = base * temp; 77 x >>= 1; 78 temp = temp * temp; 79 } 80 return base; 81 } 82 83 void print(matrix p) 84 { 85 int n = p.n; 86 int m = p.m; 87 for(int i=1;i<=n;i++) 88 { 89 for(int j=1;j<=m;j++) 90 { 91 printf("%d ",p.e[i][j]); 92 } 93 cout << endl; 94 } 95 } 96 97 void Set(matrix &temp,int x,int y,int op) 98 { 99 temp.e[x][y] = op; 100 } 101 102 matrix solve(matrix temp,int x) 103 { 104 if(x==1) return temp; 105 matrix ret = matrix(2,2); 106 ret.getE(); 107 if(x & 1) 108 { 109 return (ret+qpow(temp,x/2)) * solve(temp,x/2) + qpow(temp,x); 110 } 111 else 112 { 113 return (ret+qpow(temp,x/2)) * solve(temp,x/2); 114 } 115 } 116 117 void getAns(int r) { 118 //printf("Yes\n"); 119 int ans[205][205]; 120 memset(ans, -1, sizeof(ans)); 121 for (int i = r / 2 + 1; i <= r; i++) 122 ans[i][i - (r / 2)] = 0; 123 for (int i = r / 2 + 2; i <= r; i++) 124 for (int j = 1; j <= (i - r / 2 - 1); j++) 125 ans[i][j] = 1; 126 for (int i = r / 2 + 1; i <= r; i++) 127 for (int j = r - i + 1; j <= r; j++) 128 ans[j][i] = 1; 129 for (int i = 1; i <= r; i++) { 130 for (int j = 1; j <= r; j++) { 131 printf("%d", ans[i][j]); 132 printf(j == r ? "\n" : " "); 133 } 134 } 135 } 136 137 int main() 138 { 139 /* 140 http://www.cnblogs.com/Torrance/p/5412802.html 141 */ 142 143 int n; 144 int T;cin >>T; 145 for(int kase=1;kase<=T;kase++) 146 { 147 scanf("%d%d",&n,&mod); 148 printf("Case %d: ",kase); 149 /*matrix ans = matrix(1,2); 150 ans.e[1][1] = ans.e[1][2] = 1; 151 152 matrix temp = matrix(2,2); 153 temp.e[1][2] = temp.e[2][1] = temp.e[2][2] = 1; 154 temp = solve(temp,n-1); 155 156 matrix t = matrix(2,2);t.getE(); 157 temp = temp + t; 158 159 ans = ans * temp; 160 int sum = ans.e[1][1];*/ 161 162 // 上面的方法是这样的: 163 /* 164 (F1,F2)* temp(构造矩阵) = (F2,F3) 165 同理的:(F1,F2)* temp^2 = (F3,F4) 166 ... 167 那么只要(F1,F2)*(temp^0+temp^1+...+temp^(n-1)),这个矩阵的第一个元素即是sum 168 这是矩阵的分配律。而后面部分只要二分即可求得 169 */ 170 171 // 下面是常规的方法(?) 172 /* 173 (F1,F2,sum(初始化为F1))*temp^(n-1) 174 sum可以理解为加到F1(第一个元素)为止的元素的和,第二个元素可以理解为当前要加的元素 175 temp为: 176 0 1 0 177 1 1 1 178 0 0 1 179 */ 180 181 matrix ans = matrix(1,3); 182 for(int i=1;i<=3;i++) Set(ans,1,i,1); 183 matrix temp = matrix(3,3); 184 Set(temp,1,2,1);Set(temp,2,1,1);Set(temp,2,2,1);Set(temp,2,3,1);Set(temp,3,3,1); 185 186 temp = qpow(temp,n-1); 187 ans = ans * temp; 188 int sum = ans.e[1][3]; 189 190 //cout << sum <<endl; 191 192 cout << (sum%2 || sum ==0 ? "No":"Yes") <<endl; 193 if(sum % 2 == 0 && sum != 0) getAns(sum); 194 } 195 }
D题有个很好的博客在这儿:D题。关于Q题我找出的递推关系如下:
g(x)=g(x-1)+g(x-2)+1,(f(x)要调用一次,所以要加1)。
似乎还有一种更好的方法但是我不是很理解。。
那么,矩阵就暂时告一段落了~
————————————————隔日的分界线——————————————————————
这里写于第二天,因为发现昨天总结的过于匆忙,有个东西忘记讲了。关于E题,这种既需要1000*6的矩阵又需要6*1000的矩阵的,如果是结构体,那么就需要开1000*1000从而照顾他们都能实现,但是1000*1000的话结构体是开不下的(程序会爆栈而直接崩溃)。所以这里需要一种新的写法,vector的写法,这样的话,想开多少开多少即可。但是我不是很喜欢这种写法,所以这里只是顺带提供一下这样写法的模板,以备不时之需。E题代码如下:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <stdio.h> 2 #include <algorithm> 3 #include <math.h> 4 #include <vector> 5 #include <map> 6 #include <set> 7 #include <iostream> 8 #include <string.h> 9 using namespace std; 10 typedef pair<int,int> pii; 11 const int mod = 6; 12 13 int n,k; 14 15 typedef vector<int> vec ; 16 typedef vector<vec> mat ; 17 18 inline void add (int &a ,int b) { 19 a += b ; if(a<0) a+=mod; a%=mod ; 20 } 21 22 mat mul (const mat &a , const mat &b) { 23 mat ret(a.size(),vec(b[0].size())) ; 24 for (int i=0 ; i<a.size() ; i++) 25 for (int j=0 ; j<b[0].size() ; j++) 26 for (int z=0 ; z<b.size() ; z++) 27 add(ret[i][j],a[i][z]*b[z][j]%mod) ; 28 return ret ; 29 } 30 31 int main () { 32 while(scanf("%d%d",&n,&k)==2) 33 { 34 if(n==0 && k==0) break; 35 mat a(n,vec(k)); 36 mat b(k,vec(n)); 37 for(int i=0;i<n;i++) 38 { 39 for(int j=0;j<k;j++) scanf("%d",&a[i][j]); 40 } 41 for(int i=0;i<k;i++) 42 { 43 for(int j=0;j<n;j++) scanf("%d",&b[i][j]); 44 } 45 mat c =mul(b,a); 46 mat base(c.size(),vec(c.size())); 47 for(int i=0;i<c.size();i++) base[i][i] = 1; 48 /* 49 如果直接算A*B的话,得到的是1000*1000的矩阵,所以会一直超时。 50 因为C^(n*n)=(A*B)^(n*n)=A*(B*A)^(n*n-1)*B,由于B*A是6*6的矩阵,再用矩阵快速幂来求就行了。 51 */ 52 int temp = n*n-1; 53 while(temp) 54 { 55 if(temp & 1) base = mul(base,c); 56 temp >>= 1; 57 c = mul(c,c); 58 } 59 mat ans = mul(mul(a,base),b); 60 int sum = 0; 61 for(int i=0;i<n;i++) 62 { 63 for(int j=0;j<n;j++) sum += ans[i][j]; 64 } 65 cout << sum << endl; 66 } 67 return 0 ; 68 }
————————————————2017/9/4的分界线————————————————————
考虑到矩阵稀疏的情况,有些题目(如poj3735)需要进行稀疏矩阵的乘法优化:若矩阵A*B,枚举A的每个位置A.e[i][j],若其不为0,则能做出贡献,在乘法A.e[i][j]*B.e[j][k]做出贡献,那么接着枚举k即可。
另外上面的模板需要注意的一个地方是,由于每次乘法都会创造一个ret矩阵且伴随着一次memset,因此矩阵大小建议按照实际的大小设置(由于1base所以大小=size+1)否则可能出现TLE的情况。