文章标题 POJ : 3233 Matrix Power Series (矩阵快速幂+求等比矩阵的和）

Matrix Power Series

//http://www.cnblogs.com/jiangjing/archive/2013/05/28/3103336.html

/*
分析：求a^1+..a^n这是矩阵乘法中关于等比矩阵的求法：

|A  E|

|0  E|

其中的A为m阶矩阵，E是单位矩阵，0是零矩阵。而我们要求的是：

A^1+A^2+..A^L，由等比矩阵的性质

|A  ,  1|                 |A^n , 1+A^1+A^2+....+A^(n-1)|

|0  ,  1| 的n次方等于     |0   ,         1             |

所以我们只需要将A矩阵扩大四倍，变成如上形式的矩阵B，
然后开L+1次方就可以得到1+A^1+A^2+....+A^L。
由于多了一个1，所以最后得到的答案我们还要减去1。同理我们把矩阵A变成B：

|A  E|

|0  E|

然后我们就是求B的n+1次幂之后得到的矩阵为
|x1   x2|

|x3   x4|

右上角的矩阵x2减去单位矩阵E，得到就是要求的矩阵了！

*/

#include<iostream>
#include<string>
#include<cstdio>
#include<cstring>
#include<vector>
#include<math.h>
#include<map>
#include<queue>
#include<algorithm>
using namespace std;
const int inf = 0x3f3f3f3f;
typedef pair<int,int> pii;
typedef long long ll;

const int N = 100;
int mod;
int n,M,K;

struct Matrix {
ll mat

;
Matrix operator * (const Matrix m)const {
Matrix tmp;
for (int i=0;i<n;i++){
for (int j=0;j<n;j++){
tmp.mat[i][j]=0;
for (int k=0;k<n;k++){
tmp.mat[i][j]+=mat[i][k]*m.mat[k][j]%mod;
tmp.mat[i][j]%=mod;
}
}
}
return tmp;
}
};

Matrix Pow(Matrix &m,int k){
Matrix ans;
memset (ans.mat,0,sizeof (ans.mat));
for (int i=0;i<n;i++)ans.mat[i][i]=1;
while (k){
if (k&1)ans=ans*m;
k>>=1;
m=m*m;
}
return ans;
}

int main ()
{

while (scanf ("%d%d%d",&n,&K,&M)!=EOF){
Matrix ans;
memset (ans.mat,0,sizeof (ans.mat));
for (int i=0;i<n;i++){
for (int j=0;j<n;j++){
scanf ("%d",&ans.mat[i][j]);
}
}
for (int i=0;i<n;i++)ans.mat[i][i+n]=1;//右上角的单位矩阵
for (int i=0;i<n;i++)ans.mat[i+n][i+n]=1;//右下角的单位矩阵

n=n*2;//扩大两倍
mod = M;
ans=Pow(ans,K+1);

for (int i=0;i<n/2;i++)//右上角的矩阵减去单位矩阵
ans.mat[i][i+n/2]=(ans.mat[i][i+n/2]-1+mod)%mod;

for (int i=0;i<n/2;i++){//输出
for (int j=n/2;j<n;j++){
if (j==n/2)printf ("%lld",ans.mat[i][j]);
else printf (" %lld",ans.mat[i][j]);
}
printf ("\n");
}
}
return 0;
}