[BZOJ]4162: shlw loves matrix II

时间:2022-05-03 16:24:59

Time Limit: 30 Sec  Memory Limit: 128 MB

Description

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

Input

  第 1 行包含两个整数 n,k,其中 n 使用二进制表示,可能含有前导零;
  余下 k 行描述了一个 k * k 的矩阵 M。

Output

  输出题目描述中要求的矩阵,格式同输入。

Sample Input

  010 3
  5 9 5
  5 4 0
  8 8 8

Sample Output

  110 121 65
  45 61 25
  144 168 104

HINT

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

Solution

  矩阵乘法特征多项式优化,矩阵M的特征多项式是$f(\lambda)=|M-\lambda I|$,随便带入k+1个$\lambda$,高斯消元求出行列式的值,然后插值就能求出M的特征多项式(这里的总复杂度为$O(k^4)$,主要在计算行列式上面)。根据Cayley-Hamilton定理,$f(M)=0$,也就是说对同一个式子减去若干个$f(M)$后值不变,我们即可用M的k-1次多项式表示一个M的次幂,两个多项式相乘后对$f(M)$取模,然后就可以快速幂了,多项式取模和乘法用$O(k^2)$暴力就够了,总复杂度$O(k^4+k^2logn)$。

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
inline int read()
{
int x;char c;
while((c=getchar())<''||c>'');
for(x=c-'';(c=getchar())>=''&&c<='';)x=x*+c-'';
return x;
}
#define ML 10000
#define MK 50
#define MN 100
#define MOD 1000000007
char s[ML+];
int k,c[MK+][MK+],x[MK+],ans[MK+][MK+];
inline int mo1(int x){return x<MOD?x:x-MOD;}
inline int mo2(int x){return x<?x+MOD:x;}
int inv(int x)
{
int r=,y=MOD-;
for(;y;y>>=,x=1LL*x*x%MOD)if(y&)r=1LL*r*x%MOD;
return r;
}
struct poly
{
int n,a[MN+];
poly(int n=,int a0=,int a1=):n(n){memset(a,,sizeof(a));a[]=a0;a[]=a1;}
friend poly operator+(poly a,const poly&b)
{
a.n=max(a.n,b.n);
for(int i=;i<=a.n;++i)a.a[i]=mo1(a.a[i]+b.a[i]);
return a;
}
friend poly operator*(const poly&a,const poly&b)
{
poly c(a.n+b.n);int i,j;ll x;
for(i=;i<=c.n;c.a[i++]=x%MOD)for(x=j=;j<=i&&j<=a.n;++j)
x+=1LL*a.a[j]*b.a[i-j],x>8e18?x%=MOD:;
return c;
}
friend poly operator*(poly a,int b)
{
for(int i=;i<=a.n;++i)a.a[i]=1LL*a.a[i]*b%MOD;
return a;
}
friend poly operator/(poly a,const poly&b)
{
poly c(a.n-b.n);int i,j;
for(i=a.n;i>=b.n;--i)
{
c.a[i-b.n]=1LL*a.a[i]*inv(b.a[b.n])%MOD;
for(j=;j<=b.n;++j)a.a[i-j]=mo2((a.a[i-j]-1LL*b.a[b.n-j]*c.a[i-b.n])%MOD);
}
return c;
}
friend poly operator%(poly a,const poly&b)
{
int i,j,x;
for(i=a.n;i>=b.n;--i)
{
x=1LL*a.a[i]*inv(b.a[b.n])%MOD;
for(j=;j<=b.n;++j)a.a[i-j]=mo2((a.a[i-j]-1LL*b.a[b.n-j]*x)%MOD);
}
a.n=b.n-;
return a;
}
}f,a,p(,);
struct mat
{
int z[MK+][MK+];
mat(){memset(z,,sizeof(z));}
friend mat operator*(const mat&a,const mat&b)
{
mat c;int i,j,k;ll x;
for(i=;i<=MK;++i)for(j=;j<=MK;c.z[i][j++]=x%MOD)
for(x=k=;k<=MK;++k)x+=1LL*a.z[i][k]*b.z[k][j],x>8e18?x%=MOD:;
return c;
}
}m,n;
int cal(int x)
{
int i,j,l,ans=;
for(i=;i<=k;++i)for(j=;j<=k;++j)c[i][j]=m.z[i][j];
for(i=;i<=k;++i)c[i][i]=mo2(c[i][i]-x);
for(i=;i<=k;++i)
{
for(j=i;j<=k;++j)if(c[j][i])break;
if(j>k)return ;
if(j>i)for(ans=MOD-ans,l=i;l<=k;++l)swap(c[i][l],c[j][l]);
ans=1LL*ans*c[i][i]%MOD;
for(j=i;++j<=k;)
{
x=1LL*c[j][i]*inv(c[i][i])%MOD;
for(l=i;l<=k;++l)c[j][l]=mo2((c[j][l]-1LL*c[i][l]*x)%MOD);
}
}
return ans;
}
int main()
{
int i,j,l,v;
scanf("%s",s+);k=read();
for(i=;i<=k;++i)for(j=;j<=k;++j)m.z[i][j]=read();
for(i=;i<=k;++i)x[i]=cal(i);
for(i=;i<=k;++i)p=p*poly(,mo2(-i),);
for(i=;i<=k;++i)
{
for(v=x[i],j=;j<=k;++j)if(i!=j)v=1LL*v*inv(mo2(i-j))%MOD;
f=f+p/poly(,mo2(-i),)*v;
}
a=poly(,);p=poly(,,);
for(i=strlen(s+);i;--i,p=p*p%f)if(s[i]>'')a=a*p%f;
for(i=;i<=k;++i)n.z[i][i]=;
for(l=;l<k;++l,n=n*m)for(i=;i<=k;++i)for(j=;j<=k;++j)
ans[i][j]=(ans[i][j]+1LL*a.a[l]*n.z[i][j])%MOD;
for(i=;i<=k;printf("%d\n",ans[i++][j]))for(j=;j<k;++j)printf("%d ",ans[i][j]);
}