首页 > 代码库 > Poj3233(Matrix Power Series)

Poj3233(Matrix Power Series)

题目大意:

给一个n*n的方阵A,令S=A+A^2+...+A^k 给定k,m

输出这样的S方阵,里面的每一项模m后的结果。

方阵中的元素可以是方阵.

令S(i)=A+A^2+...+A^i 

那么观察相邻的两个S

S(i-1)=A+A^2+...+A^(i-1)

S(i)=A+A^2+...+A^(i-1)+A^i

S(i)=S(i-1)*A+A

有了这个递推关系.构造一个初始矩阵

(A E) ps:A为刚开始给定的方阵,E为单位矩阵.这个矩阵是n行2*n列

再构造一个快速幂的方阵

(A E

 A 0) ps:0代表全是0的矩阵

这是一个2*n行2*n列的方阵

然后快速幂k-1次再与初始矩阵相乘就是结果(别忘了取模)

1A代码

  1 #include <iostream>  2 #include <cstdio>  3 #include <string.h>  4 #include <stdlib.h>  5 #define maxn 65  6 using namespace std;  7 int MOD;  8 __int64 num[maxn][maxn];  9 int N; 10 class Matrix{ 11 public: 12     int row,col; 13     __int64 mapp[maxn][maxn]; 14     Matrix(int r,int c):row(r),col(c){} 15     Matrix(){} 16     void ini() 17     { 18         for(int i=0;i<2;i++) 19          for(int j=0;j<2;j++) 20          { 21            for(int q=0;q<N;q++) 22             for(int k=0;k<N;k++) 23             { 24              if((i==0&&j==0)||(i==1&&j==0)) 25               mapp[i*N+q][j*N+k]=num[q][k]; 26              else 27               if(i==1&&j==1) 28                if(q==k) mapp[i*N+q][j*N+k]=1; 29                else mapp[i*N+q][j*N+k]=0; 30               else 31                mapp[i*N+q][j*N+k]=0; 32             } 33          } 34     } 35     void unit() //单位矩阵的初始化 36     { 37         memset(mapp,0,sizeof(mapp)); 38         for(int i=0;i<row;i++) 39          mapp[i][i]=1; 40     } 41     void zero() 42     { 43         memset(mapp,0,sizeof(mapp)); 44     } 45     void result() 46     { 47         for(int i=0;i<row;i++) 48         { 49          for(int j=0;j<row;j++) 50          { 51           if(j!=0) cout<<" "; 52           cout<<mapp[i][j]; 53          } 54          cout<<endl; 55         } 56     } 57     friend Matrix operator*(const Matrix&,const Matrix&); 58     void Pow(int); 59 }; 60 Matrix operator*(const Matrix& M1,const Matrix& M2) 61 { 62     Matrix M(M1.row,M2.col); 63     for(int i=0;i<M1.row;i++) 64      for(int j=0;j<M2.col;j++) 65      { 66         M.mapp[i][j]=0; 67         for(int k=0;k<M2.col;k++) 68          M.mapp[i][j]=(M.mapp[i][j]+M1.mapp[i][k]*M2.mapp[k][j])%MOD; 69      } 70     return M; 71 } 72 Matrix M; 73 void Matrix::Pow(int n)//矩阵快速幂 74 { 75     //cout<<n<<endl; 76     Matrix temp(N*2,N*2); 77     temp.ini(); 78     while(n) 79     { 80         if(n&1) 81          M=M*temp; 82         temp=temp*temp; 83         n>>=1; 84     } 85 } 86 int main() 87 { 88     int m,k; 89     scanf("%d%d%d",&N,&k,&m); 90     MOD=m; 91     Matrix resM(N,2*N); 92     M.row=N*2;M.col=N*2; 93     M.unit(); 94     for(int i=0;i<N;i++) 95      for(int j=0;j<N;j++) 96      { 97       scanf("%I64d",&resM.mapp[i][j]); 98       num[i][j]=resM.mapp[i][j]; 99      }100     for(int i=0;i<N;i++)101      for(int j=N;j<2*N;j++)102       if(i==j-N) resM.mapp[i][j]=1;103       else resM.mapp[i][j]=0;104     if(k>1)105     {106      M.Pow(k-1);107      resM=resM*M;108     }109     resM.result();110     return 0;111 }

还可以有其他构造方法

比如

S(i)=S(i-1)+A^i

按照这样的递推式构造

前一个矩阵(A A^2)

后一个方阵(E 0

               E A)

这样应该也是可以的.

结论:不同的递推式有不用的构造方法...

 

Poj3233(Matrix Power Series)