矩阵快速幂
poj 3070 Fibonacci http://poj.org/problem?id=3070
模板题,矩阵都给你写好了。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 class Matrix { ///矩阵 5 typedef int typev;///权值类型 6 static const int MV=2;///矩阵长度 7 static const int mod=10000;///%mod 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 9 Matrix ret; 10 ret.n=a.n; 11 ret.m=b.m; 12 ret.zero(); 13 if(a.m==b.n) { 14 for(int k=0; k<a.m; k++) { 15 for(int i=0; i<a.n; i++) { 16 if(a.val[i][k]) { 17 for(int j=0; j<b.m; j++) { 18 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 19 ret.val[i][j]%=mod;///看具体需要 20 } 21 } 22 } 23 } 24 } 25 return ret; 26 } 27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 28 Matrix ret; 29 ret.n=ret.m=a.n; 30 ret.unit(); 31 for(; b; a=a*a,b>>=1) { 32 if(b&1) { 33 ret=ret*a; 34 } 35 } 36 return ret; 37 } 38 public: 39 int n,m;///n行m列 40 typev val[MV][MV]; 41 void zero() { 42 mt(val,0); 43 } 44 void unit() { 45 zero(); 46 for(int i=0; i<MV; i++) 47 val[i][i]=1; 48 } 49 } A; 50 int main() { 51 int n; 52 while(~scanf("%d",&n),n>-1) { 53 A.n=A.m=2; 54 A.val[0][0]=1; 55 A.val[0][1]=1; 56 A.val[1][0]=1; 57 A.val[1][1]=0; 58 printf("%d\n",(A^n).val[0][1]); 59 } 60 return 0; 61 }
hdu 1575 Tr A http://acm.hdu.edu.cn/showproblem.php?pid=1575
模板题
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 const int MOD=9973; 5 class Matrix { ///矩阵 下标0开始 6 typedef int typev;///权值类型 7 static const int MV=16;///矩阵长度 8 static const int mod=MOD;///%mod 9 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 10 Matrix ret; 11 ret.n=a.n; 12 ret.m=b.m; 13 ret.zero(); 14 if(a.m==b.n) { 15 for(int k=0; k<a.m; k++) { 16 for(int i=0; i<a.n; i++) { 17 if(a.val[i][k]) { 18 for(int j=0; j<b.m; j++) { 19 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 20 ret.val[i][j]%=mod;///看具体需要 21 } 22 } 23 } 24 } 25 } 26 return ret; 27 } 28 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 29 Matrix ret; 30 ret.n=ret.m=a.n; 31 ret.unit(); 32 for(; b; a=a*a,b>>=1) { 33 if(b&1) { 34 ret=ret*a; 35 } 36 } 37 return ret; 38 } 39 public: 40 int n,m;///n行m列 41 typev val[MV][MV]; 42 void zero() { 43 mt(val,0); 44 } 45 void unit() { 46 zero(); 47 for(int i=0; i<MV; i++) 48 val[i][i]=1; 49 } 50 } A; 51 int main() { 52 int t,n,k; 53 while(~scanf("%d",&t)){ 54 while(t--){ 55 scanf("%d%d",&n,&k); 56 A.n=A.m=n; 57 for(int i=0;i<n;i++){ 58 for(int j=0;j<n;j++){ 59 scanf("%d",&A.val[i][j]); 60 } 61 } 62 A=A^k; 63 int ans=0; 64 for(int i=0;i<n;i++){ 65 ans+=A.val[i][i]; 66 ans%=MOD; 67 } 68 printf("%d\n",ans); 69 } 70 } 71 return 0; 72 }
Kiki & Little Kiki 2 http://acm.hdu.edu.cn/showproblem.php?pid=2276
矩阵先不说了。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 class Matrix { ///矩阵 下标0开始 5 typedef int typev;///权值类型 6 static const int MV=128;///矩阵长度 7 static const int mod=2;///%mod 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 9 Matrix ret; 10 ret.n=a.n; 11 ret.m=b.m; 12 ret.zero(); 13 if(a.m==b.n) { 14 for(int k=0; k<a.m; k++) { 15 for(int i=0; i<a.n; i++) { 16 if(a.val[i][k]) { 17 for(int j=0; j<b.m; j++) { 18 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 19 ret.val[i][j]%=mod;///看具体需要 20 } 21 } 22 } 23 } 24 } 25 return ret; 26 } 27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 28 Matrix ret; 29 ret.n=ret.m=a.n; 30 ret.unit(); 31 for(; b; a=a*a,b>>=1) { 32 if(b&1) { 33 ret=ret*a; 34 } 35 } 36 return ret; 37 } 38 public: 39 int n,m;///n行m列 40 typev val[MV][MV]; 41 void zero() { 42 mt(val,0); 43 } 44 void unit() { 45 zero(); 46 for(int i=0; i<MV; i++) 47 val[i][i]=1; 48 } 49 } A,T; 50 int main(){ 51 int n; 52 char s[128]; 53 while(~scanf("%d%s",&n,s)){ 54 int ls=strlen(s); 55 A.n=A.m=T.n=T.m=ls; 56 A.zero(); 57 T.unit(); 58 for(int i=0;i<ls;i++){ 59 if(s[i]=='1'){ 60 A.val[0][i]=1; 61 } 62 T.val[i][i+1]=1; 63 } 64 T.val[ls-1][0]=1; 65 A=A*(T^n); 66 for(int i=0;i<ls;i++){ 67 printf("%d",A.val[0][i]); 68 } 69 puts(""); 70 } 71 return 0; 72 }
Matrix Power Series http://poj.org/problem?id=3233
测模板,具体构造再说吧。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 int mod; 5 class Matrix { ///矩阵 下标0开始 6 typedef int typev;///权值类型 7 static const int MV=64;///矩阵长度 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 9 Matrix ret; 10 ret.n=a.n; 11 ret.m=b.m; 12 ret.zero(); 13 if(a.m==b.n) { 14 for(int k=0; k<a.m; k++) { 15 for(int i=0; i<a.n; i++) { 16 if(a.val[i][k]) { 17 for(int j=0; j<b.m; j++) { 18 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 19 ret.val[i][j]%=mod;///看具体需要 20 } 21 } 22 } 23 } 24 } 25 return ret; 26 } 27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 28 Matrix ret; 29 ret.n=ret.m=a.n; 30 ret.unit(); 31 for(; b; a=a*a,b>>=1) { 32 if(b&1) { 33 ret=ret*a; 34 } 35 } 36 return ret; 37 } 38 public: 39 int n,m;///n行m列 40 typev val[MV][MV]; 41 void zero() { 42 mt(val,0); 43 } 44 void unit() { 45 zero(); 46 for(int i=0; i<MV; i++) 47 val[i][i]=1; 48 } 49 } A; 50 int main(){ 51 int n,k; 52 scanf("%d%d%d",&n,&k,&mod); 53 A.n=A.m=n<<1; 54 for(int i=0;i<n;i++){ 55 for(int j=0;j<n;j++){ 56 scanf("%d",&A.val[i][j]); 57 A.val[i][j+n]=A.val[i][j]; 58 A.val[i+n][j]=0; 59 A.val[i+n][j+n]=0; 60 if(i+n==j+n) 61 A.val[i+n][j+n]=1; 62 } 63 } 64 A=A^k; 65 for(int i=0;i<n;i++){ 66 for(int j=n;j<A.n;j++) 67 printf("%d ",A.val[i][j]); 68 puts(""); 69 } 70 return 0; 71 }
Training little cats http://poj.org/problem?id=3735
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #define mt(a,b) memset(a,b,sizeof(a)) 5 using namespace std; 6 typedef __int64 LL; 7 class Matrix { ///矩阵 下标0开始 8 typedef LL typev;///权值类型 9 static const int MV=128;///矩阵长度 10 static const int mod=10000;///%mod 11 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 12 Matrix ret; 13 ret.n=a.n; 14 ret.m=b.m; 15 ret.zero(); 16 if(a.m==b.n) { 17 for(int k=0; k<a.m; k++) { 18 for(int i=0; i<a.n; i++) { 19 if(a.val[i][k]) { 20 for(int j=0; j<b.m; j++) { 21 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 22 // ret.val[i][j]%=mod;///看具体需要 23 } 24 } 25 } 26 } 27 } 28 return ret; 29 } 30 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 31 Matrix ret; 32 ret.n=ret.m=a.n; 33 ret.unit(); 34 for(; b; a=a*a,b>>=1) { 35 if(b&1) { 36 ret=ret*a; 37 } 38 } 39 return ret; 40 } 41 public: 42 int n,m;///n行m列 43 typev val[MV][MV]; 44 void zero() { 45 mt(val,0); 46 } 47 void unit() { 48 zero(); 49 for(int i=0; i<MV; i++) 50 val[i][i]=1; 51 } 52 } A; 53 int main(){ 54 int x,y,n,m,k; 55 char op[4]; 56 while(~scanf("%d%d%d",&n,&m,&k),n+m+k){ 57 A.unit(); 58 A.n=A.m=n+1; 59 while(k--){ 60 scanf("%s%d",op,&x); 61 switch(op[0]){ 62 case 'g': A.val[0][x]++; break; 63 case 'e': for(int i=0;i<=n;i++) A.val[i][x]=0; break; 64 default: scanf("%d",&y); for(int i=0;i<=n;i++) swap(A.val[i][x],A.val[i][y]); 65 } 66 } 67 A=A^m; 68 for(int i=1;i<=n;i++) printf("%I64d ",A.val[0][i]); 69 puts(""); 70 } 71 return 0; 72 }
Fast Matrix Calculation http://acm.hdu.edu.cn/showproblem.php?pid=4965
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 typedef __int64 LL; 5 class Matrix { ///矩阵 下标0开始 6 typedef int typev;///权值类型 7 static const int MV=8;///矩阵长度 8 static const int mod=6;///%mod 9 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 10 Matrix ret; 11 ret.n=a.n; 12 ret.m=b.m; 13 ret.zero(); 14 if(a.m==b.n) { 15 for(int k=0; k<a.m; k++) { 16 for(int i=0; i<a.n; i++) { 17 if(a.val[i][k]) { 18 for(int j=0; j<b.m; j++) { 19 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 20 ret.val[i][j]%=mod;///看具体需要 21 } 22 } 23 } 24 } 25 } 26 return ret; 27 } 28 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 29 Matrix ret; 30 ret.n=ret.m=a.n; 31 ret.unit(); 32 for(; b; a=a*a,b>>=1) { 33 if(b&1) { 34 ret=ret*a; 35 } 36 } 37 return ret; 38 } 39 public: 40 int n,m;///n行m列 41 typev val[MV][MV]; 42 void zero() { 43 mt(val,0); 44 } 45 void unit() { 46 zero(); 47 for(int i=0; i<MV; i++) 48 val[i][i]=1; 49 } 50 } D; 51 struct mat { 52 int n,m; 53 int data[1024][1024]; 54 }A,B,C,E; 55 class MatrixOp { 56 public: 57 int mul(mat& c,const mat& a,const mat& b) { //c=a*b 58 int i,j,k; 59 if (a.m!=b.n) return 0; 60 c.n=a.n,c.m=b.m; 61 for (i=0; i<c.n; i++) 62 for (j=0; j<c.m; j++) 63 for (c.data[i][j]=k=0; k<a.m; k++) 64 c.data[i][j]+=a.data[i][k]*b.data[k][j]; 65 return 1; 66 } 67 } gx; 68 int main(){ 69 int N,K; 70 while(~scanf("%d%d",&N,&K),N|K){ 71 A.n=N; 72 A.m=K; 73 B.n=K; 74 B.m=N; 75 for(int i=0;i<N;i++){ 76 for(int j=0;j<K;j++){ 77 scanf("%d",&A.data[i][j]); 78 } 79 } 80 for(int i=0;i<K;i++){ 81 for(int j=0;j<N;j++){ 82 scanf("%d",&B.data[i][j]); 83 } 84 } 85 gx.mul(C,B,A); 86 D.n=D.m=C.n; 87 for(int i=0;i<C.n;i++){ 88 for(int j=0;j<C.m;j++){ 89 D.val[i][j]=C.data[i][j]; 90 } 91 } 92 D=D^(N*N-1); 93 for(int i=0;i<C.n;i++){ 94 for(int j=0;j<C.m;j++){ 95 C.data[i][j]=D.val[i][j]; 96 } 97 } 98 gx.mul(E,A,C); 99 gx.mul(C,E,B); 100 LL ans=0; 101 for(int i=0;i<C.n;i++){ 102 for(int j=0;j<C.m;j++){ 103 ans+=C.data[i][j]%6; 104 } 105 } 106 printf("%I64d\n",ans); 107 } 108 return 0; 109 }
ZOJ Problem Set - 2974 Just Pour the Water http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=1973
nenucontest3,一个序列经过相同变化,变化10^9次,明显是构造矩阵,然后快速幂。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 class Matrix { ///矩阵 下标0开始 5 typedef double typev;///权值类型 6 static const int MV=128;///矩阵长度 7 // static const int mod=10000;///%mod 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法,必须a的列等于b的行 9 Matrix ret; 10 ret.n=a.n; 11 ret.m=b.m; 12 ret.zero(); 13 if(a.m==b.n) { 14 for(int k=0; k<a.m; k++) { 15 for(int i=0; i<a.n; i++) { 16 if(a.val[i][k]) { 17 for(int j=0; j<b.m; j++) { 18 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 19 // ret.val[i][j]%=mod;///看具体需要 20 } 21 } 22 } 23 } 24 } 25 return ret; 26 } 27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 28 Matrix ret; 29 ret.n=ret.m=a.n; 30 ret.unit(); 31 for(; b; a=a*a,b>>=1) { 32 if(b&1) { 33 ret=ret*a; 34 } 35 } 36 return ret; 37 } 38 public: 39 int n,m;///n行m列 40 typev val[MV][MV]; 41 void zero() { 42 mt(val,0); 43 } 44 void unit() { 45 zero(); 46 for(int i=0; i<MV; i++) 47 val[i][i]=1; 48 } 49 } A; 50 double a[128]; 51 int main(){ 52 int t,n; 53 scanf("%d",&t); 54 while(t--){ 55 scanf("%d",&n); 56 for(int i=0;i<n;i++){ 57 scanf("%lf",&a[i]); 58 } 59 A.unit(); 60 int k,id,m; 61 for(int i=0;i<n;i++){ 62 scanf("%d",&k); 63 if(k) A.val[i][i]=0; 64 for(int j=0;j<k;j++){ 65 scanf("%d",&id); 66 A.val[id-1][i]=1.0/k; 67 } 68 } 69 A.n=A.m=n; 70 scanf("%d",&m); 71 A=A^m; 72 for(int i=0;i<n;i++){ 73 if(i) printf(" "); 74 double sum=0; 75 for(int j=0;j<n;j++){ 76 sum+=A.val[i][j]*a[j]; 77 } 78 printf("%.2f",sum); 79 } 80 puts(""); 81 } 82 return 0; 83 }
end