求一個 $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;
    }
}

bzoj 4162 shlw loves matrix II