题目大意
给你
n∗n
的矩阵
A
,求
特征多项式
这是一个特征多项式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");
}
}