https://www.acwing.com/problem/content/submission/227/

需要构造一种新的矩阵，受到前几天xy的求和的启发，但是还是不知道矩阵的求和怎么搞。事实上矩阵的求和是一样的。

构造一个矩阵：其中E是单位矩阵，O是零矩阵，那么这个东西转移n次就得到需要的Sn，而A在此过程中自动转移。

\[ \left[ \begin{matrix} S_1 \\ A^2 \\ \end{matrix} \right] = \left[ \begin{matrix} E&E \\ O&A \\ \end{matrix} \right] * \left[ \begin{matrix} S_0 \\ A^1 \\ \end{matrix} \right] \]

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

int mod;
struct Matrix30 {
    static const int MAXN = 30;
    int ma[MAXN][MAXN];
    Matrix30 () {
        init();
    }
    void init() {
        memset(ma, 0, sizeof(ma));
    }
    void setE() {
        init();
        for(int i = 0; i < MAXN; ++i)
            ma[i][i] = 1;
    }
    Matrix30 operator+(const Matrix30 &m)const {
        Matrix30 Tmp;
        Tmp.init();
        for(int i = 0; i < MAXN; ++i) {
            for(int j = 0; j < MAXN; ++j)
                Tmp.ma[i][j] += (ma[i][j] + m.ma[i][j]) % mod;
        }
        for(int i = 0; i < MAXN; ++i) {
            for(int j = 0; j < MAXN; ++j)
                if(Tmp.ma[i][j] >= mod)
                    Tmp.ma[i][j] %= mod;
        }
        return Tmp;
    }

    Matrix30 operator*(const Matrix30 &m)const {
        Matrix30 Tmp;
        Tmp.init();
        for(int k = 0; k < MAXN; ++k) {
            for(int i = 0; i < MAXN; ++i) {
                register int r = ma[i][k];
                for(int j = 0; j < MAXN; ++j)
                    Tmp.ma[i][j] += (r * m.ma[k][j]) % mod;
            }
        }
        for(int i = 0; i < MAXN; ++i) {
            for(int j = 0; j < MAXN; ++j)
                if(Tmp.ma[i][j] >= mod)
                    Tmp.ma[i][j] %= mod;
        }
        return Tmp;
    }
    void show(int n) {
        for(int i = 0; i < n; ++i) {
            for(int j = 0; j < n; ++j)
                printf("%d%c", ma[i][j], " \n"[j == n - 1]);
        }
    }
} E, O, A;

struct Matrix2 {
    static const int MAXN = 2;
    Matrix30 ma[MAXN][MAXN];
    Matrix2 () {
        init();
    }
    void init() {
        ma[0][0] = O;
        ma[0][1] = O;
        ma[1][0] = O;
        ma[1][1] = O;
    }
    void setE() {
        ma[0][0] = E;
        ma[0][1] = O;
        ma[1][0] = O;
        ma[1][1] = E;
    }
    Matrix2 operator*(const Matrix2 &m)const {
        Matrix2 Tmp;
        Tmp.init();
        for(int k = 0; k < MAXN; ++k) {
            for(int i = 0; i < MAXN; ++i)
                for(int j = 0; j < MAXN; ++j)
                    Tmp.ma[i][j] = Tmp.ma[i][j] + (ma[i][k] * m.ma[k][j]);
        }
        return Tmp;
    }
};

Matrix2 qpow(Matrix2 x, int n) {
    Matrix2 res;
    res.setE();
    while(n) {
        if(n & 1)
            res = res * x;
        x = x * x;
        n >>= 1;
    }
    return res;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in", "r", stdin);
#endif // Yinku
    int n, k;
    scanf("%d%d%d", &n, &k, &mod);
    E.setE();
    for(int i = 0; i < n; ++i) {
        for(int j = 0; j < n; ++j) {
            scanf("%d", &A.ma[i][j]);
            A.ma[i][j] %= mod;
        }
    }

    Matrix2 NXT;
    NXT.ma[0][0] = E;
    NXT.ma[0][1] = E;
    NXT.ma[1][0] = O;
    NXT.ma[1][1] = A;

    Matrix2 I;
    I.ma[0][0] = O;
    I.ma[0][1] = O;
    I.ma[1][0] = A;
    I.ma[1][1] = O;

    Matrix2 RES = qpow(NXT, k);
    RES = RES * I;
    RES.ma[0][0].show(n);
}