\[ 注:W是一个矩阵,表示题中w[i][j],下列式子中的W做加减法时表示W_{xy}\\ ans_t=\Sigma_{i\ mod\ k=t}^L{L \choose i}W^i\\ =\Sigma_{i=0}^L[k|(i-t)]{L \choose i}W^i\\ =\Sigma_{i=0}^L\frac1k\Sigma_{j=0}^{k-1}(\omega_k^{i-t})^j{L \choose i}W^i\\ =\Sigma_i^L\frac1k\Sigma_{j=0}^{k-1}\omega_k^{ij-tj}{L \choose i}W^i\\ =\frac1k\Sigma_{j=0}^{k-1}\omega_k^{

{t \choose 2}+{j \choose 2}-{t+j \choose 2}}\Sigma_{i=0}^L\omega_k^{ij}{L \choose i}W^i\\ 设C(j)=\Sigma_{i=0}^L\omega_k^{ij}{L \choose i}W^i\\ 把它变成二项式的形式:\\ A(j)=\Sigma_{i=0}^L{L \choose i}(W\omega_k^j)^i*1^{L-i}\\ \therefore A(j)=(W\omega_k^j+1)^L\\ \therefore ans_t=\frac{\omega_k^{t \choose 2}}k\Sigma_{j=0}^{k-1}\omega_k^{j \choose 2}A(j)\omega_k^{-{t+j \choose 2}}\\ 这好像是还不是卷积,但是我们可以把它变成卷积\\ 令f(x)=\omega_k^{x \choose 2}\\ 令g(x)=A(x)\omega_k^{-{t+x \choose 2}}\\ 那么按照常规操作,将f倍增,再将g翻转.\\ \therefore ans_t=\frac{\omega_k^{t \choose 2}}k(f*g)(k+t)\\ 然后卷积用MTT算,就没了. \]

#pragma GCC optimize(3)#include
    
     #define N (65536*4)using namespace std;const double pi=acos(-1);int n,k,l,x,y,p;struct Matrix{    int a[3][3];    void format(){        memset(a,0,sizeof(a));    }    friend Matrix operator + (const Matrix &a,const Matrix &b){        Matrix re;        for(int i=0;i<3;i++)        for(int j=0;j<3;j++){            re.a[i][j]=(a.a[i][j]+b.a[i][j])%p;        }        return re;    }    friend Matrix operator * (const Matrix &a,const Matrix &b){        Matrix re;        re.format();        for(int i=0;i<3;i++)        for(int j=0;j<3;j++)        for(int k=0;k<3;k++){            re.a[i][j]=(re.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%p;        }        return re;    }    friend Matrix operator * (const int &a,const Matrix &b){        Matrix re;        for(int i=0;i<3;i++)        for(int j=0;j<3;j++){            re.a[i][j]=(1ll*a*b.a[i][j])%p;        }        return re;    }    friend Matrix operator * (const Matrix &a,const int &b){return b*a;}}I,W;Matrix Pow(Matrix x,int y){    Matrix re=I;    while(y){        if(y&1)re=re*x;        y>>=1;        x=x*x;    }    return re;}struct Complex{double x,y;};Complex operator + (Complex a,Complex b){return (Complex){a.x+b.x,a.y+b.y};}Complex operator - (Complex a,Complex b){return (Complex){a.x-b.x,a.y-b.y};}Complex operator * (Complex a,Complex b){return (Complex){a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};}Complex a[N],b[N],c[N];int Pow(int x,int y){    int re=1;    while(y){        if(y&1)re=(1ll*re*x)%p;        y>>=1;        x=(1ll*x*x)%p;    }    return re;}int rev[N];int pre(int n,int m){    int high=0,cnt=1;;    while(cnt
     
      <<=1,high++;    for(int i=0;i
      
       >1]>>1)|((i&1)<<(high-1));    return cnt;}int root(int num){    int div[128],sum=0;    int phi=num-1;    for(int i=2;i<=phi;i++){        if(phi%i==0){            div[++sum]=i;            while(phi%i==0)phi/=i;        }    }    if(phi!=1)div[++sum]=phi;    int re=2;    while(1){        int flag=1;        for(int i=1;i<=sum;i++){            int pw=Pow(re,(num-1)/div[i]);            if(pw==1){flag=0;break;}        }        if(flag)return re;        re++;    }}Complex w[N];void fft(int lim,Complex *buf,int dft=1){    for(int i=0;i
       
        
         >15,B1[i].x=A[i]&32767;    for(int i=0;i
         
          >15,B2[i].x=B[i]&32767; for(int i=0;i