poj 3233 S = A + A^2 + A^3 + … + A^k A是一个n X n矩阵 （矩阵快速幂）

S = A + A^2 + A^3 + … + A^k A是一个n*n矩阵

Sample Input

2 2 4 //n k MOD
0 1
1 1
Sample Output

1 2
2 3

 

先求 I + A + A^2 + A^3 + … + A^k  I为单位矩阵

 

我们来设置这样一个矩阵 B=
A I
O I
其中O是零矩阵，I是单位矩阵

将它乘方，得到
A^2 I+A
O I
乘三方，得到
A^3 I+A+A^2
O I
乘四方，得到
A^4 I+A+A^2+A^3
O I

然后B^(k+1) 的右上小矩阵 就是 I + A + A^2 + A^3 + … + A^k

记得 数组开大一点...因为要扩展成2n*2n的矩阵， 第一次开小了，数组越界=.=

 

 1 # include <iostream>
 2 # include <cstdio>
 3 # include <algorithm>
 4 # include <cmath>
 5 # define LL long long
 6 using namespace std ;
 7 
 8 int MOD ;
 9 int n ;
10 
11 struct Matrix
12 {
13     LL mat[70][70];
14 };
15 
16 Matrix mul(Matrix a,Matrix b)
17 {
18     Matrix c;
19     for(int i=0;i<2*n;i++)
20         for(int j=0;j<2*n;j++)
21         {
22             c.mat[i][j]=0;
23             for(int k=0;k<2*n;k++)
24             {
25                 c.mat[i][j]=(c.mat[i][j] + a.mat[i][k]*b.mat[k][j])%MOD;
26             }
27         }
28     return c;
29 }
30 Matrix pow_M(Matrix a,int k )
31 {
32     Matrix ans;
33     memset(ans.mat,0,sizeof(ans.mat));
34     for (int i=0;i<2*n;i++)
35         ans.mat[i][i]=1;
36     Matrix temp=a;
37     while(k)
38     {
39         if(k&1)ans=mul(ans,temp);
40         temp=mul(temp,temp);
41         k>>=1;
42     }
43     return ans;
44 }
45 
46 
47 
48 int main ()
49 {
50     //freopen("in.txt","r",stdin) ;
51     int  k ;
52     while(cin>>n>>k>>MOD)
53     {
54         Matrix A ;
55         int i , j ;
56         for (i = 0 ; i < n ; i++)
57             for (j = 0 ; j < n ; j++)
58                cin>>A.mat[i][j] ;
59         Matrix B ;
60         memset(B.mat,0,sizeof(B.mat));
61         for (i = 0 ; i < n ; i++) //扩展成2n * 2n的矩阵
62         {
63             for (j = 0 ; j < n ; j++)
64             {
65                 B.mat[i][j] = A.mat[i][j] ;
66             }
67             B.mat[n+i][n+i] = 1 ;
68             B.mat[i][n+i] = 1 ;
69 
70         }
71         B = pow_M(B,k+1) ;
72         LL t ;
73         for (i = 0 ; i < n ; i++)
74             for (j = 0 ; j < n ; j++)
75             {
76                 t = B.mat[i][n+j] % MOD ;
77                 if (i == j)
78                     t = (t + MOD - 1)%MOD ;
79                 if (j+1 != n)
80                     cout<<t<<" " ;
81                 else
82                     cout<<t<<endl ;
83             }
84 
85     }
86 
87 
88 
89     return 0 ;
90 }
View Code

 

posted @ 2015-05-29 18:40  __Meng  阅读(381)  评论(0)  编辑  收藏  举报