求一个 $m \times 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;
- }
- }
- View Code
来源: http://www.bubuko.com/infodetail-2983818.html