来自FallDream的博客，未经允许，请勿转载，谢谢
给定矩阵k*k的矩阵M，请计算 M^n，并将其中每一个元素对 1000000007 取模输出。 k<=50 n<=2^10000
考虑求出矩阵的特征多项式，这点我们可以通过带入$\lambda=x0..xk$，求出矩阵的行列式，然后通过插值求出多项式。
然后搬出一个很厉害的定理:f(A)=0 其中f(x)是特征多项式，A是矩阵，所以我们可以把所求的$x^{n}$对f(x)取膜，从而让答案变成一个k-1次多项式。然后我们把原矩阵带进去就行了。
插值的复杂度是$O(n^{4})$,然后后面那部分的复杂度是$k^{2}logn$
然后做了这道题好像懂了怎么在O(klogklogn)内做完"常系数线性递推"，貌似还要nlogn的多项式取余 真麻烦TAT
#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
#define mod 1000000007
using namespace std;
inline int read()
{
int x = ; char ch = getchar();
while(ch < '' || ch > '')ch = getchar();
while(ch >= '' && ch <= ''){x = x * + ch - '';ch = getchar();}
return x;
} char st[];
int y[],k; int pow(int x,int k)
{
int sum=;
for(;k;k>>=,x=1LL*x*x%mod)
if(k&) sum=1LL*sum*x%mod;
return sum;
}
struct Matrix
{
int s[][];
Matrix operator*(Matrix b)
{
Matrix c;memset(c.s,,sizeof(c.s));
for(int i=;i<=k;i++)
for(int l=;l<=k;l++) if(s[i][l])
for(int j=;j<=k;j++)
c.s[i][j]=(c.s[i][j]+1LL*s[i][l]*b.s[l][j])%mod;
return c;
}
int calc()
{
int sum=,j;
for(int i=;i<=k;i++)
{
for(j=i;j<=k;j++)
if(s[j][i])
{
if(j!=i)
{
sum=mod-sum;
for(int l=;l<=k;l++)
swap(s[j][l],s[i][l]);
}
break;
}
if(j==k+) return ;
for(int j=i+;j<=k;j++)
{
int inv=1LL*pow(s[i][i],mod-)*s[j][i]%mod;
for(int l=i;l<=k;l++)
s[j][l]=(1LL*s[i][l]*inv%mod-s[j][l]+mod)%mod;
}
}
for(int i=;i<=k;i++) sum=1LL*sum*s[i][i]%mod;
return sum;
}
}a,b,ans; struct poly
{
int s[];
poly(int x=){memset(s,,sizeof(s));s[]=x;}
poly operator^(int x)
{
poly c;
for(int i=k<<;i;i--)
c.s[i]=(s[i-]+1LL*x*s[i])%mod;
c.s[]=(1LL*s[]*x)%mod;
return c;
}
poly operator*(int x)
{
poly c();
for(int i=;i<=k<<;i++)
c.s[i]=1LL*s[i]*x%mod;
return c;
}
poly operator+(poly b)
{
poly c();
for(int i=;i<=k<<;i++)
c.s[i]=(s[i]+b.s[i])%mod;
return c;
}
poly operator*(poly b)
{
poly c();
for(int i=;i<=k;i++)
for(int j=;j<=k;j++)
c.s[i+j]=(c.s[i+j]+1LL*s[i]*b.s[j])%mod;
return c;
}
friend poly operator%(poly a,poly b)
{
for(int i=k;i>=;i--)
{
int t=1LL*a.s[i+k]*pow(b.s[k],mod-)%mod;
for(int j=;j<=k;j++)
a.s[i+j]=(a.s[i+j]-1LL*b.s[j]*t%mod+mod)%mod;
}
return a;
}
}F(),Ans(); int main()
{
scanf("%s",st+);k=read();
for(register int i=;i<=k;i++)
for(register int j=;j<=k;j++)
b.s[i][j]=a.s[i][j]=read();
for(register int t=;t<=k;t++)
{
memcpy(b.s,a.s,sizeof(b.s));
for(int i=;i<=k;i++)
b.s[i][i]=(b.s[i][i]-t+mod)%mod;
y[t]=b.calc();
}
for(register int t=;t<=k;t++)
{
poly tmp();
for(register int i=;i<=k;i++) if(i!=t)
tmp=tmp^(mod-i),tmp=tmp*pow((t-i+mod)%mod,mod-);
tmp=tmp*y[t];F=F+tmp;
}
poly x();x.s[]=;
for(int i=strlen(st+);i;i--)
{
if(st[i]=='') Ans=Ans*x%F;
x=x*x%F;
}
memset(b.s,,sizeof(b.s));
for(int i=;i<=k;i++) b.s[i][i]=;
for(int t=;t<k;t++,b=a*b)
for(int i=;i<=k;i++)
for(int j=;j<=k;j++)
ans.s[i][j]=(ans.s[i][j]+1LL*Ans.s[t]*b.s[i][j])%mod;
for(register int i=;i<=k;i++)
for(register int j=;j<=k;j++)
printf("%d%c",ans.s[i][j],j!=k?' ':'\n');
return ;
}