poj3233_矩阵快速幂水题

被这水题卡了半天，这题主要思路是矩阵快速幂：

A^1+A^2+...+A^k=(A^1+..+A^(k/2))+(A^(k/2))(+...+A^k)=(E+A^(k/2))*(A^1+..+A^(k/2));而对(A^1+..+A^(k/2))可以采用相同的方法；注意是0矩阵还是e（单位矩阵）；注意要long long

#include<cstdio>
#include<cstring>
const int maxn=35;

int n,mod;

struct matrix{
    long long a[maxn][maxn];
    matrix(){
        memset(a,0,sizeof(a));
    }
};

void matrixPrint(matrix a)
{
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            printf("%lld%c",a.a[i][j],j==n-1?'\n':' ');
}

matrix matMul(matrix a,matrix b)
{
    matrix ans;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
            ans.a[i][j]=0;
            for(int k=0;k<n;k++)
                ans.a[i][j]+=a.a[i][k]*b.a[k][j];
            ans.a[i][j]%=mod;
        }
    return ans;
}

matrix matPls(matrix a, matrix b)
{
    matrix ans;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
            ans.a[i][j]=a.a[i][j]+b.a[i][j];
            ans.a[i][j]%=mod;
        }
   return ans;
}

matrix fastMod(matrix a,int k)
{
    matrix ans;
    for(int i=0;i<n;i++)ans.a[i][i]=1;
    while(k){
        if(k&1)ans=matMul(ans,a);
        a=matMul(a,a);
        k=k>>1;
    }
    return ans;
}

matrix fun(matrix a,int k)
{
    if(k==1)return a;
    matrix ans;
    memset(ans.a,0,sizeof(ans.a));
    for(int i=0;i<n;i++)ans.a[i][i]=1;
    ans=matPls(ans,fastMod(a,k>>1));
    ans=matMul(ans,fun(a,k>>1));
    if(k&1)ans=matPls(ans,fastMod(a,k));
    return ans;
}

int main()
{
    int k;
    scanf("%d%d%d",&n,&k,&mod);
    matrix a,ans;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            scanf("%lld",&a.a[i][j]);
    ans=fun(a,k);
    matrixPrint(ans);
    return 0;
}