矩阵快速幂

poj 3070 Fibonacci http://poj.org/problem?id=3070

模板题，矩阵都给你写好了。

 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 }

View Code

hdu 1575 Tr A http://acm.hdu.edu.cn/showproblem.php?pid=1575

模板题

 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 }

View Code

Kiki & Little Kiki 2 http://acm.hdu.edu.cn/showproblem.php?pid=2276

矩阵先不说了。

 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 }

View Code

Matrix Power Series http://poj.org/problem?id=3233

测模板，具体构造再说吧。

 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 }

View Code

Training little cats http://poj.org/problem?id=3735

 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 }

View Code

Fast Matrix Calculation http://acm.hdu.edu.cn/showproblem.php?pid=4965

  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 }

View Code

ZOJ Problem Set - 2974 Just Pour the Water http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=1973

nenucontest3,一个序列经过相同变化，变化10^9次，明显是构造矩阵，然后快速幂。

 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 }

View Code

end

posted on 2014-09-01 20:58 gaolzzxin 阅读(174) 评论(0) 编辑收藏举报