poj Matrix Power Series 矩阵幂求和

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef struct node
{
    int matrix[55][55];
}Matrix;
Matrix a,sa,unit;
int n,m,k,i,j;
Matrix add(Matrix a,Matrix b)//矩阵加
{
    Matrix c;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            c.matrix[i][j]=a.matrix[i][j]+b.matrix[i][j];
            c.matrix[i][j]%=m;
        }
    }
    return c;
}
Matrix mul(Matrix a,Matrix b)//矩阵乘
{
    Matrix c;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            c.matrix[i][j]=0;
            for(int h=0;h<n;h++)
            {
                c.matrix[i][j]=c.matrix[i][j]+a.matrix[i][h]*b.matrix[h][j];
            }
            c.matrix[i][j]%=m;
        }
    }
    return c;
}
Matrix cal(int e)//矩阵幂
{
    Matrix p,q;
    p=a;q=unit;
    while(e!=1)
    {
        if(e&1)
        {
            e--;
            q=mul(p,q);
        }
        else
        {
            e/=2;
            p=mul(p,p);
        }
    }
    p=mul(p,q);
    return p;
}
Matrix sum(int k)//幂之和
{
    if(k==1) return a;
    Matrix temp,tnow;
    temp=sum(k/2);
    if(k&1)
    {
        tnow=cal(k/2+1);
        temp=add(temp,mul(temp,tnow));
        temp=add(tnow,temp);
    }
    else
    {
        tnow=cal(k/2);
        temp=add(temp,mul(temp,tnow));
    }
    return temp;
}
int main()
{
    int i,j,l;
    while(scanf("%d%d%d",&n,&k,&m)!=EOF)
    {
      for(i=0;i<n;i++)
      {
        for(j=0;j<n;j++)
        {
            scanf("%d",&a.matrix[i][j]);
            a.matrix[i][j]%=m;
            unit.matrix[i][j]=(i==j);//单位矩阵
        }
      }
      sa=sum(k);//幂为k的结果
      for(i=0;i<n;i++)
      {
        for(j=0;j<n-1;j++)
        {
            printf("%d ",sa.matrix[i][j]%m);//对m取模
        }
        printf("%d\n",sa.matrix[i][n-1]%m);
      }
    }
    return 0;
}