[bzoj4162]shlw loves matrix II

题目大意

给你$n*n$的矩阵$A$，求$A^m$。

特征多项式

这是一个特征多项式IDE练习题。
矩阵A的特征多项式f(x)为det(A-Ix)。
可以发现f(A)=0。
如何求f(x)？代入n+1个点值求行列式，再插值插出f(x)。
设g(x)=x^m mod f(x)。
可以发现g(A)=A^m mod f(A)=A^m。
因此我们可以多项式快速幂+取模求出g(x)，再代入A即可。

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const int maxn=50+10,mo=1000000007;
int a[maxn][maxn],b[maxn][maxn],ans[maxn][maxn],o[maxn][maxn],x[maxn],y[maxn];
int f[maxn*2],g[maxn*2],h[maxn*2],sta[10000+10];
int i,j,k,l,t,n,m,top;
char ch;
bool czy;
int qsm(int x,int y){
    if (!y) return 1;
    int t=qsm(x,y/2);
    t=(ll)t*t%mo;
    if (y%2) t=(ll)t*x%mo;
    return t;
}
int det(){
    int i,j,k,t,cnt=1;
    fo(i,1,n){
        fo(j,i+1,n)
            if (b[j][i]){
                t=(ll)b[i][i]*qsm(b[j][i],mo-2)%mo;
                fo(k,i,n) (b[i][k]-=(ll)b[j][k]*t%mo)%=mo;
                fo(k,i,n) swap(b[i][k],b[j][k]);
                cnt=-cnt;
            }
        if (!b[i][i]) return 0;
        cnt=(ll)cnt*b[i][i]%mo;
    }
    return cnt;
}
void work(){
    int i,j,k,t=qsm(f[n],mo-2);
    fd(i,2*n,n){
        k=(ll)h[i]*t%mo;
        fo(j,0,n) (h[i-j]-=(ll)f[n-j]*k%mo)%=mo;
    }
    fo(i,0,n) g[i]=h[i];
}
int main(){
    freopen("data.in","r",stdin);freopen("data.out","w",stdout);
    top=0;
    ch=getchar();
    while (ch>='0'&&ch<='1'){
        if (ch=='1') czy=1;
        if (czy) sta[++top]=ch-'0';
        ch=getchar();
    }
    scanf("%d",&n);
    if (!czy){
        fo(i,1,n){
            fo(j,1,n) printf("1 ");
            printf("\n");
        }
        return 0;
    }
    reverse(sta+1,sta+top+1);
    fo(i,1,n)
        fo(j,1,n)
            scanf("%d",&a[i][j]);
    fo(t,0,n){
        x[t]=t;
        fo(i,1,n)
            fo(j,1,n)
                b[i][j]=a[i][j];
        fo(i,1,n) b[i][i]-=t;
        y[t]=det();
    }
    fo(i,0,n){
        t=y[i];
        fo(j,1,n) g[j]=0;
        g[0]=1;
        fo(j,0,n)
            if (i!=j){
                t=(ll)t*qsm(x[i]-x[j],mo-2)%mo;
                fo(k,1,n) h[k]=g[k-1];
                h[0]=0;
                fo(k,0,n) (h[k]-=(ll)g[k]*x[j]%mo)%=mo;
                fo(k,0,n) g[k]=h[k];
            }
        fo(j,0,n) g[j]=(ll)g[j]*t%mo;
        fo(j,0,n) (f[j]+=g[j])%=mo;
    }
    fo(i,1,n) g[i]=0;
    g[0]=1;
    while (top){
        fo(i,0,2*n) h[i]=0;
        fo(i,0,n)
            fo(j,0,n)
                (h[i+j]+=(ll)g[i]*g[j]%mo)%=mo;
        work();
        if (sta[top]){
            fo(i,0,2*n) h[i]=0;
            fo(i,0,n) h[i+1]=g[i];
            work();
        }
        top--;
    }
    fo(i,1,n)
        fo(j,1,n) b[i][j]=0;
    fo(i,1,n) b[i][i]=1;
    fo(l,0,n-1){
        fo(i,1,n)
            fo(j,1,n)
                (ans[i][j]+=(ll)g[l]*b[i][j]%mo)%=mo;
        fo(i,1,n)
            fo(j,1,n) o[i][j]=0;
        fo(k,1,n)
            fo(i,1,n)
                fo(j,1,n)
                    (o[i][j]+=(ll)b[i][k]*a[k][j]%mo)%=mo;
        fo(i,1,n)
            fo(j,1,n) b[i][j]=o[i][j];
    }
    fo(i,1,n){
        fo(j,1,n){
            (ans[i][j]+=mo)%=mo;
            printf("%d ",ans[i][j]);
        }
        printf("\n");
    }
}