传送门
利用Cayley-Hamilton定理，用插值法求出特征多项式 \(P(x)\)
然后 \(M^n\equiv M^n(mod~P(x))(mod~P(x))\)
然后就多项式快速幂+取模
最后得到了一个关于 \(M\) 的多项式，代入 \(M^i\) 即可
# include <bits/stdc++.h>
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;
}