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的矩阵， 第一次开小了，数组越界=.=

# include <iostream>
# include <cstdio>
# include <algorithm>
# include <cmath>
# define LL long long
using namespace std ;

int MOD ;
int n ;

struct Matrix
{
    LL mat[70][70];
};

Matrix mul(Matrix a,Matrix b)
{
    Matrix c;
    for(int i=0;i<2*n;i++)
        for(int j=0;j<2*n;j++)
        {
            c.mat[i][j]=0;
            for(int k=0;k<2*n;k++)
            {
                c.mat[i][j]=(c.mat[i][j] + a.mat[i][k]*b.mat[k][j])%MOD;
            }
        }
    return c;
}
Matrix pow_M(Matrix a,int k )
{
    Matrix ans;
    memset(ans.mat,0,sizeof(ans.mat));
    for (int i=0;i<2*n;i++)
        ans.mat[i][i]=1;
    Matrix temp=a;
    while(k)
    {
        if(k&1)ans=mul(ans,temp);
        temp=mul(temp,temp);
        k>>=1;
    }
    return ans;
}



int main ()
{
    //freopen("in.txt","r",stdin) ;
    int  k ;
    while(cin>>n>>k>>MOD)
    {
        Matrix A ;
        int i , j ;
        for (i = 0 ; i < n ; i++)
            for (j = 0 ; j < n ; j++)
               cin>>A.mat[i][j] ;
        Matrix B ;
        memset(B.mat,0,sizeof(B.mat));
        for (i = 0 ; i < n ; i++) //扩展成2n * 2n的矩阵
        {
            for (j = 0 ; j < n ; j++)
            {
                B.mat[i][j] = A.mat[i][j] ;
            }
            B.mat[n+i][n+i] = 1 ;
            B.mat[i][n+i] = 1 ;

        }
        B = pow_M(B,k+1) ;
        LL t ;
        for (i = 0 ; i < n ; i++)
            for (j = 0 ; j < n ; j++)
            {
                t = B.mat[i][n+j] % MOD ;
                if (i == j)
                    t = (t + MOD - 1)%MOD ;
                if (j+1 != n)
                    cout<<t<<" " ;
                else
                    cout<<t<<endl ;
            }

    }



    return 0 ;
}