BZOJ4162：shlw loves matrix II

传送门
利用Cayley-Hamilton定理，用插值法求出特征多项式 $P(x)$
然后 $M^n\equiv M^n(modP(x))(modP(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;
}