bzoj 4162: shlw loves matrix II 拉格朗日插值法+矩阵乘法

题意

给定矩阵 M，请计算 M^n，并将其中每一个元素对 1000000007 取模输出。

对于 100% 数据，满足 n <= 2^10000;k <= 50; 0 <= Mij < 10^9 +7

分析

我们可以带入k+1个不同的λ算出所有的|A−λI|，然后用拉格朗日插值法把M的特征多项式f(x)插出来。

因为根据Cayley-Hamilton定理有f(M)=0，所以xn=xnmodf(x)，用快速幂算出来模意义下的值后把A带进去算即可。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=55;
const int MOD=1000000007;

int n,a

,b[N*2],c

,d

,r[N*2],tmp[N*2],f
,g
,mat

,ans

;
pair<int,int> pts
;
char str[10005];

int read()
{
int x=0,f=1;char ch=getchar();
while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
return x*f;
}

int ksm(int x,int y)
{
int ans=1;
while (y)
{
if (y&1) ans=(LL)ans*x%MOD;
x=(LL)x*x%MOD;y>>=1;
}
return ans;
}

int det()
{
int ans=1;
for (int i=1;i<=n;i++)
{
int l=i;
for (int j=i;j<=n;j++) if (mat[j][i]) {l=j;break;}
if (l!=i)
{
for (int j=i;j<=n;j++) swap(mat[i][j],mat[l][j]);
ans=-ans;
}
for (int j=i+1;j<=n;j++)
if (mat[j][i])
{
int w=(LL)mat[j][i]*ksm(mat[i][i],MOD-2)%MOD;
for (int k=i;k<=n;k++) (mat[j][k]-=(LL)mat[i][k]*w%MOD)%=MOD;
}
ans=(LL)ans*mat[i][i]%MOD;
}
ans+=ans<0?MOD:0;
return ans;
}

void mul(int *a,int *b,int *c)
{
for (int i=0;i<=n*2-2;i++) tmp[i]=0;
for (int i=0;i<n;i++)
for (int j=0;j<n;j++)
(tmp[i+j]+=(LL)a[i]*b[j]%MOD)%=MOD;
int ny=ksm(f
,MOD-2);
for (int i=n*2-2;i>=n;i--)
for (int j=n-1;j>=0;j--)
(tmp[i-n+j]-=(LL)f[j]*tmp[i]%MOD*ny%MOD)%=MOD;
for (int i=0;i<n;i++) c[i]=tmp[i];
}

int main()
{
scanf("%s",str);n=read();
for (int i=1;i<=n;i++)
for (int j=1;j<=n;j++)
a[i][j]=read();
for (int i=0;i<=n;i++)
{
memcpy(mat,a,sizeof(a));
for (int j=1;j<=n;j++) mat[j][j]-=i;
pts[i]=make_pair(i,det());
}
for (int i=0;i<=n;i++)
{
for (int j=0;j<=n;j++) g[j]=0;
g[0]=pts[i].second;
for (int j=0;j<=n;j++) if (i!=j) g[0]=(LL)g[0]*ksm(pts[j].first-pts[i].first,MOD-2)%MOD;
for (int j=0;j<=n;j++)
if (i!=j)
{
for (int k=n;k>=1;k--) g[k]=((LL)g[k]*pts[j].first%MOD-g[k-1])%MOD;
g[0]=(LL)g[0]*pts[j].first%MOD;
}
for (int j=0;j<=n;j++) (f[j]+=g[j])%=MOD;
}
r[0]=1;b[1]=1;
int len=strlen(str);
for (int i=len-1;i>=0;i--)
{
if (str[i]=='1') mul(r,b,r);
mul(b,b,b);
}
for (int i=1;i<=n;i++) c[i][i]=1;
for (int i=0;i<n;i++)
{
for (int j=1;j<=n;j++)
for (int k=1;k<=n;k++)
(ans[j][k]+=(LL)c[j][k]*r[i]%MOD)%=MOD;
memset(d,0,sizeof(d));
for (int j=1;j<=n;j++)
for (int k=1;k<=n;k++)
for (int l=1;l<=n;l++)
(d[j][k]+=(LL)c[j][l]*a[l][k]%MOD)%=MOD;
memcpy(c,d,sizeof(d));
}
for (int i=1;i<=n;i++)
{
for (int j=1;j<n;j++) printf("%d ",(ans[i][j]+MOD)%MOD);
printf("%d\n",(ans[i]
+MOD)%MOD);
}
return 0;
}