我们可以对这个等比数列进行分治求和，通过函数 calc(int x) 计算 $A^1 + A^2 \cdots + A^x$的和
当x 为奇数时:
$$ A^1 + A^2 \cdots + A^x = \left(A^1 \cdots + A^{\frac{x}{2}}\right) + A^{\frac{x}{2}}\left( A^1 + \cdots + A^{\frac{x}{2}}\right) = (E + A^{\frac{x}{2}})\left( A^1 + \cdots + A^{\frac{x}{2}}\right) $$
同理当x 为偶数时：
$$ A^1 + A^2 + A^x = \left(A^1 \cdots + A^{\frac{x - 1}{2}}\right) + A^{\frac{x - 1}{2}}\left( A^1 + \cdots + A^{\frac{x - 1}{2}}\right) + A^k = (E + A^{\frac{x - 1}{2}})\left( A^1 + \cdots + A^{\frac{x - 1}{2}}\right) + A^k $$

通过递归和分治实现这个过程, 时间复杂度$O(n^3log^2k)$, 因为矩阵乘法的时间复杂度是n的三次方,每次递归种还要用一次快速幂.

C++ 代码

// #pragma GCC optimize(2)
#include<bits/stdc++.h>
#define xx first
#define yy second
using namespace std;
using ll = long long;
using ld = long double;
using PII = pair<ll, ll>;
const ll N = 1e5+7;
const ll inf = 1e18;
const ld eps = 1e-8;
// const ll mod = 1e9+7;

int n, k, mod;

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

    Matrix() : n(0) {}
    Matrix(int len) : n(len), a(len, vector<int>(len)){}

    void read(){
        for (int i = 0; i < n; i ++){
            for (int j = 0; j < n; j ++) cin >> a[i][j];
        }
    }

    Matrix operator +(const Matrix &y) const{
        Matrix res(n);
        for (int i = 0; i < n; i ++){
            for (int j = 0; j < n; j ++){
                res.a[i][j] = (this->a[i][j] + y.a[i][j]) % mod;
            }
        }
        return res;
    }

    Matrix operator *(const Matrix &y) const{
        Matrix res(n);
        for (int i = 0; i < n; i ++){
            for (int j = 0; j < n; j ++){
                for (int k = 0; k < n; k ++){
                    res.a[i][j] = (res.a[i][j] + this->a[i][k] * y.a[k][j]) % mod;
                }
            }
        }
        return res;
    }
};

Matrix E, A;

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

Matrix calc(int x){
    if (x == 1) return A;
    if (x & 1) return (E + qpow(A, (x - 1) / 2)) * calc((x - 1) / 2) + qpow(A, x);
    else return (E + qpow(A, x / 2)) * calc(x / 2);
}

void solve()
{ 
    cin >> n >> k >> mod;
    A = E = Matrix(n);
    A.read();
    for (int i = 0; i < n; i ++) E.a[i][i] = 1;
    Matrix ans = calc(k);
    for (int i = 0; i < n; i ++){
        for (int j = 0; j < n; j ++)
            cout << ans.a[i][j] << ' ';
        cout << '\n';
    }
}

signed main()
{
    ios::sync_with_stdio(0), cin.tie(0);
    cout.tie(0);
    ll tt = 1;
    // cin >> tt;
    while(tt --) solve();
    return 0;
}