我们需要求(g(P))，其中(g)是一个只有一项的多项式
暴力是(n ^ 3 log k)的，显然过不去
怎么办？
特征方程优化矩阵快速幂～
考虑方程(|xI - P| = 0)，把det展开后可以得到一个方程(f(x) = 0)，这个方程称为(P)的特征方程，(f)称为(P)的特征多项式
基于某个厉害的定理，有(f(P) = 0)，也就是说矩阵(P)是特征方程的一个解
于是有(g(P) = (g mod f)(P))，由于(f)的次数只有(n)，而且多项式取模的复杂度很容易做到(O(n ^ 2))
因此只要算出特征多项式就可以(O(n ^ 4 + n ^2 log k))计算
那么怎么求出特征多项式呢？
牛顿恒等式～
对于多项式(F(x) = sum_{i = 0}^{n} c_{n-i} x^i)
若有方程(F(x) = 0)，它显然有(n)个根 (x_1, x_2 ... x_n)，令(s_{k} = sum_{i = 1}^{n} x_{i}^k)
则有：
当(1 leq k leq n)时，有(sum_{i = 0}^{k - 1} c_{i} * s_{k - i} +c_{k} * k = 0)
当(n < k)时，有(sum_{i = 0}^{n} c_{i} * s_{k - i} = 0)
这玩意对于矩阵同样成立～
基于另一个厉害的定理，特征方程(f(x) = 0)的根(lambda_{i})满足：(s_k = sum lambda _{i}^{k} = tr(P^k))
于是可以通过牛顿恒等式解出特征方程每一项的系数～
这样就可以(O(n^4))解出特征多项式了了……
时间复杂度(O(n ^ 4 + n ^2 log k))
#include <bits/stdc++.h> using namespace std; // #define int long long const int N = 55; const int MOD = 1000000007; char st[20000]; int k; struct mat { int d[N][N]; mat() { memset(d, 0, sizeof d); } } X, P[N], I, ans; mat operator * (mat a, mat b) { mat c; for (int i = 0; i < k; ++ i) for (int j = 0; j < k; ++ j) for (int p = 0; p < k; ++ p) c.d[i][j] = (1ll * a.d[i][p] * b.d[p][j] + c.d[i][j]) % MOD; return c; } mat operator * (mat a, int b) { mat c; for (int i = 0; i < k; ++ i) for (int j = 0; j < k; ++ j) c.d[i][j] = 1ll * a.d[i][j] * b % MOD; return c; } mat operator + (mat a, mat b) { mat c; for (int i = 0; i < k; ++ i) for (int j = 0; j < k; ++ j) c.d[i][j] = (a.d[i][j] + b.d[i][j]) % MOD; return c; } int powi(int a, int b) { int c = 1; for (; b; b /= 2, a = 1ll * a * a % MOD) if (b & 1) c = 1ll * c * a % MOD; return c; } int S[N], C[N]; struct poly { int d[N * 2]; poly () { memset(d, 0, sizeof d); } } pI, pS; poly operator * (poly a, poly b) { poly c; for (int i = 0; i < k; ++ i) for (int j = 0; j < k; ++ j) c.d[i + j] = (c.d[i + j] + 1ll * a.d[i] * b.d[j]) % MOD; for (int i = k + k - 2; i >= k; -- i) { int v = c.d[i]; for (int j = 0; j <= k; ++ j) c.d[i - j] = (c.d[i - j] - 1ll * v * C[j] % MOD + MOD) % MOD; } return c; } poly powi(poly a, char b[], int l) { poly c = pI; for (int i = l - 1; ~i; -- i, a = a * a) if (b[i] == '1') c = c * a; return c; } signed main() { scanf("%s%d", st, &k); int l = strlen(st); for (int i = 0; i < k; ++ i) I.d[i][i] = 1; pI.d[0] = 1; for (int i = 0; i < k; ++ i) for (int j = 0; j < k; ++ j) scanf("%d", &X.d[i][j]); P[0] = I; C[0] = 1; for (int i = 1; i <= k; ++ i) { P[i] = P[i - 1] * X; for (int j = 0; j < k; ++ j) S[i] = (S[i] + P[i].d[j][j]) % MOD; for (int j = 1; j <= i; ++ j) C[i] = (C[i] - 1ll * C[j - 1] * S[i - j + 1] % MOD + MOD) % MOD; C[i] = 1ll * C[i] * powi(i, MOD - 2) % MOD; } pS.d[1] = 1; pS = powi(pS, st, l); for (int i = 0; i < k; ++ i) ans = ans + P[i] * pS.d[i]; for (int i = 0; i < k; ++ i) { for (int j = 0; j < k - 1; ++ j) printf("%d ", ans.d[i][j]); printf("%d", ans.d[i][k - 1]); puts(""); } }