| A E |
| 0 E | 

| A E |^2    | A^2   A + E|
| 0 E |    = | 0     E    |

| A E |^3    | A^3   A^2 + A + E|
| 0 E |    = | 0     E          |

应该能发现了了吧,构造矩阵的 (k + 1) 次方的右上角 减去 单位矩阵E 就是答案

直接qpow算就完事

#include <bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
using namespace std;
typedef long long ll;

const int N = 31;

int mod, n, k;

struct Matrix {
    int n;
    vector<vector<ll>> a;

    Matrix(int x) {
        this->n = x;
        vector<vector<ll>>(x, vector<ll>(x, 0ll)).swap(a);
    }

    Matrix() {}

    Matrix operator*(const Matrix& y) {
        Matrix c(this->n);
        rep(i, 0, this->n - 1)
            rep(j, 0, this->n - 1)
                rep(k, 0, this->n - 1)
                    c.a[i][j] = (c.a[i][j] + this->a[i][k] * y.a[k][j] % mod) % mod;
        return c;
    }

}t;

Matrix qpow(Matrix a, int b) {
    Matrix ans(n << 1);
    rep(i, 0, ans.n - 1) ans.a[i][i] = 1;
    for (; b; a = a * a, b >>= 1)
        if (b & 1) ans = ans * a;
    return ans;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0);
    cin >> n >> k >> mod;
    t = Matrix(n << 1);
    rep(i, 0, n - 1) t.a[i + n][i + n] = t.a[i][i + n] = 1;

    rep(i, 0, n - 1)
        rep(j, 0, n - 1)
            cin >> t.a[i][j];

    t = qpow(t, k + 1);
    rep(i, 0, n - 1) t.a[i][i + n] = (t.a[i][i + n] - 1 + mod) % mod;

    rep(i, 0, n - 1) {
        rep(j, 0, n - 1) cout << t.a[i][j + n] << ' ';
        cout << '\n';
    }
    return 0;
}