利用Cayley-Hamilton定理，用插值法求出特征多项式 \(P(x)\)

然后 \(M^n\equiv M^n(mod~P(x))(mod~P(x))\)

然后就多项式快速幂+取模

最后得到了一个关于 \(M\) 的多项式，代入 \(M^i\) 即可

# include 
    
     using namespace std;typedef long long ll;const int mod(1e9 + 7);inline int Pow(ll x, int y) {    register ll ret = 1;    for (; y; y >>= 1, x = x * x % mod)        if (y & 1) ret = ret * x % mod;    return ret;}inline void Inc(int &x, int y) {    x = x + y >= mod ? x + y - mod : x + y;}inline int Dec(int x, int y) {    return x - y < 0 ? x - y + mod : x - y;}int n, m, a[55][55], b[55][55], mt[55][55], tmt[55][55], len, c[55], d[55], p[55], tmp[105], yi[55];char str[10005];inline int Gauss() {    register int i, j, k, inv, ans = 1;    for (i = 1; i <= n; ++i) {        for (j = i; j <= n; ++j)            if (b[j][i]) {                if (i != j) swap(b[i], b[j]), ans = mod - ans;                break;            }        for (j = i + 1; j <= n; ++j)            if (b[j][i]) {                inv = (ll)b[j][i] * Pow(b[i][i], mod - 2) % mod;                for (k = i; k <= n; ++k) Inc(b[j][k], mod - (ll)b[i][k] * inv % mod);            }        ans = (ll)ans * b[i][i] % mod;    }    return ans;}inline void Mul(int *x, int *y, int *z) {    register int i, j, inv;    memset(tmp, 0, sizeof(tmp));    for (i = 0; i <= n; ++i)        for (j = 0; j <= n; ++j) Inc(tmp[i + j], (ll)x[i] * y[j] % mod);    for (i = m; i >= n; --i) {        inv = (ll)tmp[i] * Pow(p[n], mod - 2);        for (j = 0; j <= n; ++j) Inc(tmp[i - j], mod - (ll)p[n - j] * inv % mod);    }    for (i = 0; i <= n; ++i) z[i] = tmp[i];}int main() {    register int i, j, k, l, inv;    scanf(" %s%d", str + 1, &n), len = strlen(str + 1), m = n << 1;    for (i = 1; i <= n; ++i)        for (j = 1; j <= n; ++j) scanf("%d", &a[i][j]);    for (i = 0; i <= n; ++i) {        memset(b, 0, sizeof(b));        for (j = 1; j <= n; ++j)            for (k = 1; k <= n; ++k)                b[j][k] = (j ^ k) ? mod - a[j][k] : Dec(i, a[j][k]);        yi[i] = Gauss();    }    for (i = 0; i <= n; ++i) {        memset(tmp, 0, sizeof(tmp)), tmp[0] = yi[i];        for (j = 0; j <= n; ++j)            if (j ^ i) {                for (k = n; k; --k) tmp[k] = Dec(tmp[k - 1], (ll)tmp[k] * j % mod);                tmp[0] = mod - (ll)tmp[0] * j % mod, inv = Pow(Dec(i, j), mod - 2);                for (k = 0; k <= n; ++k) tmp[k] = (ll)tmp[k] * inv % mod;            }        for (j = 0; j <= n; ++j) Inc(p[j], tmp[j]);    }    c[0] = d[1] = 1;    for (i = len; i; --i) {        if (str[i] == '1') Mul(c, d, c);        Mul(d, d, d);    }    memset(b, 0, sizeof(b));    for (i = 1; i <= n; ++i) mt[i][i] = 1;    for (l = 0; l <= n; ++l) {        for (i = 1; i <= n; ++i)            for (j = 1; j <= n; ++j)                Inc(b[i][j], (ll)c[l] * mt[i][j] % mod);        memset(tmt, 0, sizeof(tmt));        for (i = 1; i <= n; ++i)            for (j = 1; j <= n; ++j)                for (k = 1; k <= n; ++k)                    Inc(tmt[i][k], (ll)mt[i][j] * a[j][k] % mod);        memcpy(mt, tmt, sizeof(mt));    }    for (i = 1; i <= n; ++i, putchar('\n'))        for (j = 1; j <= n; ++j) printf("%d ", b[i][j]);    return 0;}