個人認為這這種解法是這道題比較妙的解法, 因為它是從矩陣加速遞推的板子中合理衍生得到的.
設f[n]=A^1+A^2+A^3+...+A^n
則f[n]=A+f[n-1]*A
轉化成矩陣加速遞推(矩陣套矩陣, 1為單位矩陣):
[f[n],1] = [f[n-1],1] * [A 0]
...............................[A 1]
所以:
[f[n],1] = [A,1] * [A 0]^(k-1)
.........................[A 1]
成為板子題

#include<cstdio>
#include<cstring>
#define in inline
#define re register
#define ll long long
#define db double
int n,k,yyb;
struct jz{
    int a[31][31];
    jz operator * (const jz &b)const
        {
            jz c;
            memset(c.a,0,sizeof(c.a));
            for(re int i=1;i<=n;++i)
                for(re int j=1;j<=n;++j)
                    for(re int k=1;k<=n;++k)
                        c.a[i][j]+=(ll)a[i][k]*(ll)b.a[k][j]%(ll)yyb;
            return c;
        }
    void operator += (const jz &b)
        {
            for(re int i=1;i<=n;++i)
                for(re int j=1;j<=n;++j)
                    a[i][j]=(a[i][j]+b.a[i][j])%yyb;
        }
}A,I;
struct jzjz{
    jz a[3][3];//這樣實現矩陣套矩陣
    jzjz operator * (const jzjz &b)const
        {
            jzjz c;
            memset(c.a,0,sizeof(c.a));
            for(re int i=1;i<=2;++i)
                for(re int j=1;j<=2;++j)
                    for(re int k=1;k<=2;++k)
                        c.a[i][j]+=a[i][k]*b.a[k][j];
            return c;
        }
}cs,zy,II;
jzjz ksm(int k)
{
    if(k==0) return II;
    if(k&1) return zy*ksm(k-1);
    jzjz res=ksm(k>>1);
    return res*res;
}
signed main()
{
    scanf("%d%d%d",&n,&k,&yyb);
    for(re int i=1;i<=n;++i)
        for(re int j=1;j<=n;++j)
            scanf("%d",&A.a[i][j]);
    for(re int i=1;i<=n;++i) I.a[i][i]=1;
    II.a[1][1]=I,II.a[2][2]=I;//構造單位矩陣
        cs.a[1][1]=A,cs.a[1][2]=I;//構造初始矩陣
    zy.a[1][1]=A,zy.a[2][1]=A,zy.a[2][2]=I;//構造轉移矩陣
    jzjz ans=cs*ksm(k-1);
    for(re int i=1;i<=n;++i)
    {
        for(re int j=1;j<=n;++j) printf("%d ",ans.a[1][1].a[i][j]);
        printf("\n");
    }
    return 0;
}
 

Matrix Power Series[矩陣加速]