来自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 = 0; char ch = getchar(); while(ch < '0' || ch > '9')ch = getchar(); while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0';ch = getchar();} return x; } char st[10001]; int y[51],k; int pow(int x,int k) { int sum=1; for(;k;k>>=1,x=1LL*x*x%mod) if(k&1) sum=1LL*sum*x%mod; return sum; } struct Matrix { int s[51][51]; Matrix operator*(Matrix b) { Matrix c;memset(c.s,0,sizeof(c.s)); for(int i=1;i<=k;i++) for(int l=1;l<=k;l++) if(s[i][l]) for(int j=1;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=1,j; for(int i=1;i<=k;i++) { for(j=i;j<=k;j++) if(s[j][i]) { if(j!=i) { sum=mod-sum; for(int l=1;l<=k;l++) swap(s[j][l],s[i][l]); } break; } if(j==k+1) return 0; for(int j=i+1;j<=k;j++) { int inv=1LL*pow(s[i][i],mod-2)*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=1;i<=k;i++) sum=1LL*sum*s[i][i]%mod; return sum; } }a,b,ans; struct poly { int s[101]; poly(int x=0){memset(s,0,sizeof(s));s[0]=x;} poly operator^(int x) { poly c; for(int i=k<<1;i;i--) c.s[i]=(s[i-1]+1LL*x*s[i])%mod; c.s[0]=(1LL*s[0]*x)%mod; return c; } poly operator*(int x) { poly c(0); for(int i=0;i<=k<<1;i++) c.s[i]=1LL*s[i]*x%mod; return c; } poly operator+(poly b) { poly c(0); for(int i=0;i<=k<<1;i++) c.s[i]=(s[i]+b.s[i])%mod; return c; } poly operator*(poly b) { poly c(0); for(int i=0;i<=k;i++) for(int j=0;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>=0;i--) { int t=1LL*a.s[i+k]*pow(b.s[k],mod-2)%mod; for(int j=0;j<=k;j++) a.s[i+j]=(a.s[i+j]-1LL*b.s[j]*t%mod+mod)%mod; } return a; } }F(0),Ans(1); int main() { scanf("%s",st+1);k=read(); for(register int i=1;i<=k;i++) for(register int j=1;j<=k;j++) b.s[i][j]=a.s[i][j]=read(); for(register int t=0;t<=k;t++) { memcpy(b.s,a.s,sizeof(b.s)); for(int i=1;i<=k;i++) b.s[i][i]=(b.s[i][i]-t+mod)%mod; y[t]=b.calc(); } for(register int t=0;t<=k;t++) { poly tmp(1); for(register int i=0;i<=k;i++) if(i!=t) tmp=tmp^(mod-i),tmp=tmp*pow((t-i+mod)%mod,mod-2); tmp=tmp*y[t];F=F+tmp; } poly x(0);x.s[1]=1; for(int i=strlen(st+1);i;i--) { if(st[i]=='1') Ans=Ans*x%F; x=x*x%F; } memset(b.s,0,sizeof(b.s)); for(int i=1;i<=k;i++) b.s[i][i]=1; for(int t=0;t<k;t++,b=a*b) for(int i=1;i<=k;i++) for(int j=1;j<=k;j++) ans.s[i][j]=(ans.s[i][j]+1LL*Ans.s[t]*b.s[i][j])%mod; for(register int i=1;i<=k;i++) for(register int j=1;j<=k;j++) printf("%d%c",ans.s[i][j],j!=k?' ':' '); return 0; }