文章标题 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[N][N];
    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;
}