2 2 4
0 1
1 1

1 2
2 3

    解題報告： 求矩陣的n次方和。

    如果接觸過矩陣，都知道矩陣快速冪。一開始我也這麼做的，二分計算前半段和後半段，即(A^1+A^2+...+A^(n/2))*(A^(n/2)+I)。題目是A了，但是耗時1600MS+。

    看了一下Status，發現一堆0MS的……又想到整數的n次方和是有公式的，那麼矩陣有嗎？

把問題轉化以加速，令

B = A  I

      0  I

則B^(k + 1) = A^(k + 1)      I + A + A² + A³ + … + A^k

                            0                          I

用二分法求B^(k + 1)

    根據該性質，複雜度大大減小。另外矩陣相乘時每次最後再取模，會快上很多。

    程式碼如下：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;

typedef long long LL;
const int SIZE = 60;
int n, mod;

struct Matrix
{
    int a[SIZE][SIZE];
    Matrix(int t=0)
    {
        memset(a, 0, sizeof(a));
        for(int i=0;i<SIZE;i++) a[i][i]=t;
    }

    Matrix operator*(const Matrix& b) const
    {
        Matrix c;
        for(int i=0;i<n;i++) for(int j=0;j<n;j++)
        {
            LL sum=0;
            for(int k=0;k<n;k++) sum+=a[i][k]*b.a[k][j];
            c.a[i][j]=sum%mod;
        }
        return c;
    }
};

Matrix pow(Matrix a, int b)
{
    Matrix res(1);
    while(b)
    {
        if(b&1)
            res=res*a;
        a=a*a;
        b>>=1;
    }
    return res;
}

int main()
{
    int k;
    while(~scanf("%d%d%d", &n, &k, &mod))
    {
        Matrix a;
        for(int i=0;i<n;i++) for(int j=0;j<n;j++)
            scanf("%d", &a.a[i][j]), a.a[i][j]%=mod;

        for(int i=0;i<n;i++)
            a.a[i][i+n]=a.a[i+n][i+n]=1;

        n<<=1;
        a = pow(a, k+1);
        n>>=1;

        for(int i=0;i<n;i++)
            if(a.a[i][i+n]==0)
                a.a[i][i+n]=mod-1;
            else
                a.a[i][i+n]--;

        for(int i=0;i<n;i++, puts("")) for(int j=0;j<n;j++)
            printf("%d ", a.a[i][j+n]);
    }
}

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;

int n, mod;

struct Matrix
{
    int a[30][30];
    Matrix(int t=0)
    {
        memset(a, 0, sizeof(a));
        for(int i=0;i<30;i++) a[i][i]=t;
    }

    Matrix operator*(const Matrix& b) const
    {
        Matrix c;
        for(int i=0;i<n;i++) for(int j=0;j<n;j++)
            for(int k=0;k<n;k++) c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod;
        return c;
    }

    Matrix operator+(const Matrix& b) const
    {
        Matrix c;
        for(int i=0;i<n;i++) for(int j=0;j<n;j++)
            c.a[i][j]=(a[i][j]+b.a[i][j])%mod;
        return c;
    }
} one(1);

Matrix pow(Matrix a, int b)
{
    Matrix res(1);
    while(b)
    {
        if(b&1)
            res=res*a;
        a=a*a;
        b>>=1;
    }
    return res;
}

Matrix powSum(Matrix a, int b)
{
    if(b==0)
        return one;
    else if(b==1)
        return a;
    else if(b%2==0)
        return powSum(a, b/2)*(pow(a, b/2)+one);
    else
        return powSum(a, b/2)*(pow(a, b/2)+one)+pow(a, b);
}

int main()
{
    int k;
    while(~scanf("%d%d%d", &n, &k, &mod))
    {
        Matrix a;
        for(int i=0;i<n;i++) for(int j=0;j<n;j++)
            scanf("%d", &a.a[i][j]), a.a[i][j]%=mod;

        a = powSum(a, k);
        for(int i=0;i<n;i++) for(int j=0;j<n;j++)
            printf("%d", a.a[i][j]%mod), printf(j==n-1?"\n":" ");
    }
}