求一个 $m imes m$ 矩阵的 $n$ 次方
$m leq 50,n leq 2^{10000}$
sol:
特征多项式是 $f(x) = |det(Ix - A)|$,插出来
然后答案就是 $A^n space mod space f(A)$
$A^n$ 是一个多项式,$f(A)$ 是一个多项式,可以用多项式快速幂+多项式取模在 $O(m^2logn)$
多项式取模就暴力吧
假设模出来是 $g(x)$,我们算 $g(A)$ 就可以了
#include <bits/stdc++.h> #define LL long long #define rep(i,s,t) for(register int i = (s),i##end = (t); i <= i##end; ++i) #define dwn(i,s,t) for(register int i = (s),i##end = (t); i >= i##end; --i) using namespace std; inline int read() { int x=0,f=1;char ch; for(ch=getchar();!isdigit(ch);ch=getchar())if(ch=='-')f=-f; for(;isdigit(ch);ch=getchar())x=10*x+ch-'0'; return x*f; } char n[10010]; int k; const int mod = 1000000007; inline int ksm(int x, int t) { int res = 1; while(t) { if(t & 1) res = 1LL * res * x % mod; x = 1LL * x * x % mod; t = t >> 1; }return res; } int a[555][555], y[555], b[555][555], c[555][555], ans[555][555], d[555][555]; int g[555], f[555], h[555], prod[555], x[555], cur[555]; int det(int x) { int ans = 1; rep(i, 1, k) rep(j, 1, k) b[i][j] = (i == j) ? (a[i][j] - x) : a[i][j]; rep(i, 1, k) rep(j, 1, k) if(b[i][j] < 0) b[i][j] += mod; rep(i, 1, k) { rep(j, i+1, k) { int t = 1LL * b[j][i] * ksm(b[i][i], mod - 2) % mod; rep(l, i, k) { (b[j][l] -= (1LL * t * b[i][l] % mod)) %= mod; if(b[j][l] < 0) b[j][l] += mod; } } } rep(i, 1, k) ans = 1LL * ans * b[i][i] % mod; ans = ((ans % mod) + mod) % mod; return ans; } void mul(int *res, int *a, int *b) { rep(i, 0, k) rep(j, 0, k) (res[i + j] += (1LL * a[i] * b[j] % mod)) %= mod; } void Mod(int *a, int *res) { int t = ksm(f[k], mod - 2); dwn(i, k+k, k) { int nk = 1LL * a[i] * t % mod; rep(j, 0, k) (a[i - j] -= (1LL * nk * f[k - j] % mod)) %= mod; } rep(i, 0, k) res[i] = (a[i] + mod) % mod; } int main() { // freopen("a.in","r",stdin); // freopen("a.out","w",stdout); scanf("%s", n+1); k = read(); rep(i, 1, k) rep(j, 1, k) a[i][j] = read(); rep(i, 0, k) y[i] = det(i); rep(i, 0, k) { int t = y[i]; rep(j, 0, k) g[j] = 0; g[0] = 1; rep(j, 0, k) { if(i != j) { t = 1LL * ksm(i - j, mod - 2) * t % mod; rep(l, 1, k) h[l] = g[l - 1]; h[0] = 0; rep(l, 0, k) (h[l] -= 1LL * g[l] * j % mod) %= mod; rep(l, 0, k) g[l] = h[l]; } } rep(j, 0, k) g[j] = 1LL * t * g[j] % mod; rep(j, 0, k) (f[j] += g[j]) %= mod; } int len = strlen(n + 1), gw = 0; rep(i, 1, len) if(n[i] == '1') gw = 1; if(!gw) { rep(i, 1, k) { rep(j, 1, k) cout << 1 << " "; cout << endl; } return 0; } //rep(i, 0, k) cout << f[i] << " "; // cout << endl; //prod[0] = 1; x[0] = 1; rep(i, 1, len) { rep(i, 0, k+k) cur[i] = 0; mul(cur, x, x); Mod(cur, x); if(n[i] == '1') { rep(i, 0, k+k) cur[i] = 0; rep(i, 0, k) cur[i+1] = x[i]; Mod(cur, x); } //rep(i, 1, k+k) cur[i] = 0; } //rep(i, 0, k) cout << x[i] << " "; //cout << endl; rep(i, 1, k) c[i][i] = 1; rep(l, 0, k-1) { rep(i, 1, k) rep(j, 1, k) (ans[i][j] += 1LL * x[l] * c[i][j] % mod) %= mod; rep(i, 1, k) rep(j, 1, k) d[i][j] = 0; rep(i, 1, k) rep(j, 1, k) rep(l, 1, k) (d[i][j] += (1LL * c[i][l] * a[l][j] % mod)) %= mod; rep(i, 1, k) rep(j, 1, k) c[i][j] = d[i][j]; } rep(i, 1, k) { rep(j, 1, k) cout << ((ans[i][j] % mod) + mod) % mod << " "; cout << endl; } }
坑:BM 算法
可以用 $O(n^2)$ 时间求一个数列的最短线性递推式
也就是对于 ${a_i}$,你可以求长度为 $m$ 的递推式 ${r_i}$ 使得 $a_n = sumlimits_{i=1}^{m} r_i imes a_{n-i}$
这个东西就是矩阵乘法,也就是上面那一道特征多项式
也就是:如果一道题是线性递推式且递推式不是很长(考试时间内跑得出来),那我们有理由在考场上 $A$ 掉它
现场可能要保持一下 想其他题 / 想这道题 / 打 BM 三者的均衡
但网络赛可能就是:
看出 BM -> 打表 -> BM -> 过题 -> 涨 rating
但还没学,有时间整一个
另外,还要学习的就是利用上海森堡矩阵 $O(m^3)$ 求化零多项式,这样就很强了,有板子=打完表一秒切掉一道很 nb 的题