特征多项式优化矩阵快速幂

F(X)=C0Xn+C1Xn-1+……+Cn-1X1+Cn（C0≠0）为矩阵的特征多项式

牛顿恒等式

C0Sk+C1Sk-1+……+Ck-1S1+kCk=0 （当1≤k≤n）

其中Sk为矩阵的k次方的主对角线的值的和

即可n^4求出矩阵的特征多项式

一个矩阵的x次方对矩阵的特征多项式取模表示为一个矩阵的0-n-1次方乘系数后的和

这一步可以n^2*logx求出系数多项式

暴力求和即可，复杂度n^4

 BZOJ4162

#include <cstdio>
#include <cstring>
#include <algorithm>
#define LL long long 
using namespace std;

  const LL mo=1e9+7;
  LL rev_for_CK,rev_for_S0,C[101],n,s[101],ans[101][101];
  char st[10001];

  LL qpow(LL bas,int powe){
      LL ret=1;
      for (;powe;bas*=bas,bas%=mo){
        if (powe&1) ret*=bas,ret%=mo;
      powe>>=1;    
    }
    return(ret);
  }

  struct matrix{
    int n,m;
    LL a[101][101];LL tmp[101][101];
          
    void mul(matrix &b){
      for (int i=0;i<=n;i++) 
        for (int j=0;j<=b.m;j++) 
          tmp[i][j]=0;
        
      for (int i=0;i<=n;i++)
        for (int k=0;k<=m;k++)
          if (a[i][k]) 
            for (int j=0;j<=b.m;j++)
              tmp[i][j]+=a[i][k]*b.a[k][j]%mo,tmp[i][j]%=mo;
        
      for (int i=0;i<=n;i++)
        for (int j=0;j<=b.m;j++)
          a[i][j]=tmp[i][j];
    }
  }bas,a;
  
  struct poly{
    int a[3001];LL tmp[3001];
    int k;
     
    void mul(poly&b){
      for (int i=0;i<=2*k-2;i++) tmp[i]=0;
      for (int i=0;i<=k-1;i++)
        for (int j=0;j<=k-1;j++)
          tmp[i+j]+=1LL*a[i]*b.a[j]%mo;
      for (int i=0;i<=2*k-2;i++) tmp[i]%=mo;
      
      for (int i=2*k-2;i>=k;i--){
          tmp[i]%=mo;
        int bas=1LL*tmp[i]*rev_for_CK%mo;
        for (int j=k-1;j>=0;j--)
          tmp[i-(k-j)]-=1LL*C[j]*bas%mo;
      }
      for (int i=k-1;i>=0;i--) a[i]=tmp[i]%mo; 
    }
  }re,ba;
   
  int main(){      
      scanf("%s",&st);scanf("%d",&n);
      bas.n=bas.m=n;
      for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
          scanf("%d",&bas.a[i][j]);
      a.n=n;a.m=n;for (int i=1;i<=n;i++) a.a[i][i]=1;
      for (int i=0;i<=n;i++){
        for (int j=1;j<=n;j++) s[i]+=a.a[j][j],s[i]%=mo;
      a.mul(bas);    
    }
    
    C[0]=1;
    for (int i=1;i<=n;i++){
      LL tot=0;
      for (int j=0;j<i;j++)
        tot+=C[j]*s[i-j]%mo,tot%=mo;
      tot*=-1;tot%=mo;tot+=mo;tot%=mo;
      tot*=qpow(i,mo-2);
      C[i]=tot%mo;
    }
    reverse(C,C+n+1);
    rev_for_CK=qpow(C[n],mo-2);
    
    re.k=ba.k=n;re.a[0]=1;ba.a[1]=1;
    int tlen=strlen(st);
    for (int i=tlen-1;i>=0;i--){
      if (st[i]=='1') re.mul(ba);
      ba.mul(ba);        
    }
    
    memset(a.a,0,sizeof(a.a));
    a.n=n;a.m=n;for (int i=1;i<=n;i++) a.a[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]+=a.a[j][k]*re.a[i]%mo,ans[j][k]%=mo;
      a.mul(bas);    
    }
    
    for (int i=1;i<=n;i++){
      for (int j=1;j<n;j++) printf("%lld ",(ans[i][j]%mo+mo)%mo);
      printf("%lld
",(ans[i][n]%mo+mo)%mo);
    }
  }