矩阵是干什么的呢？一句话来说就是，知道相邻两个函数的递推关系和第一个数，让你递推到第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在这里指的是同阶的单位矩阵。代码如下（也可以直接当做矩阵的模板了）：
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?' ':' '); 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题。这题的代码如下：
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?' ':' '); 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题还是留个代码吧：
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 "); 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 ? " " : " "); 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题代码如下：
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的情况。