一道清真的数论题

LOJ #3058

Luogu P5293

---

题解

考虑$ n=1$的时候怎么做

设$ s$为转移的方案数

设答案多项式为$\sum\limits_{i=0}^L (sx)^i\binom{L}{i}=(sx+1)^L$

答案相当于这个多项式模$ k$的各项系数的和

发现这和LJJ学二项式定理几乎一模一样

我上一题的题解

然而直接搞是$ k^2$的，无法直接通过本题

 

以下都用$ w$表示$ k$次单位根

设$ F_i$为次数模$ k$为$ i$的项的系数和

单位根反演一下得到$F(i)=\sum\limits_{j=0}^{k-1}w^{-ij}(sw^j+1)^L$

组合数有个非常优美的性质

$$ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}$$

推波式子

$$
\begin{aligned}
F(i)&=\sum_{j=0}^{k-1}w^{-ij}(sw^j+1)^L\\
&=\sum_{j=0}^{k-1}w^{\binom{i}{2}+\binom{j}{2}-\binom{i+j}{2}}(sw^j+1)^L\\
&=w^\binom{i}{2}\sum_{j=0}^{k-1}w^{-\binom{i+j}{2}}w^\binom{j}{2}(sw^j+1)^L\\
\end{aligned}
$$

发现这是一个差卷积的形式

扔一个$ MTT$上去就能拿那$ 40$分了

 

考虑$ n>1$的情况

相当于把$ s$从一个数变成了矩阵，把$ 1$变成单位矩阵

$ (sw^j+1)^L$这个矩阵我们只需要关注一个位置上的值

因此可以乘出来然后取[x][y]这个位置上的数即可

这样就可以通过此题

复杂度：大常数单 $\log$

---

 

代码

#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
    ll x=0;char zf=1;char ch=getchar();
    while(ch!='-'&&!isdigit(ch))ch=getchar();
    if(ch=='-')zf=-1,ch=getchar();
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();return x*zf;
}
void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
void writeln(const ll y){write(y);putchar('\n');}
int k,m,n,x,y,cnt,p,L;
int w[65539],val[65539];
int f[65539],g[65539],F[65539];
int ksm(int x,int y=p-2){
    int ans=1;
    for(;y;y>>=1,x=1ll*x*x%p)if(y&1)ans=1ll*ans*x%p;
    return ans;
}
struct mat{
    ll a[3][3];    
    void print(){
        for(rt i=0;i<n;i++)for(rt j=0;j<n;j++)cout<<a[i][j]<<" \n"[j==n-1];
    }
    inline mat operator*(const mat s)const{
        mat ret={0};
        for(rt i=0;i<n;i++)
        for(rt k=0;k<n;k++)
        for(rt j=0;j<n;j++)
        (ret.a[i][j]+=a[i][k]*s.a[k][j]);
        for(rt i=0;i<n;i++)for(rt j=0;j<n;j++)ret.a[i][j]%=p;
        return ret;
    }
    mat operator*(const int s)const{
        mat ret;
        for(rt i=0;i<n;i++)
        for(rt j=0;j<n;j++)
        ret.a[i][j]=a[i][j]*s%p;
        return ret;
    }
    mat operator+(const mat s)const{
        mat ret;memset(ret.a,0,sizeof(ret.a));
        for(rt i=0;i<n;i++)
        for(rt j=0;j<n;j++)
        ret.a[i][j]=(a[i][j]+s.a[i][j])%p;
        return ret;
    }

}I,zy;
mat ksm(mat x,int y){
    mat ans=I;
    for(;y;y>>=1,x=x*x)if(y&1)ans=ans*x;
    return ans;
}
int V(int x){
    return 1ll*x*(x-1)/2%k;
}
const double PI=acos(-1.0);
struct cp{
    double x,y;
    cp operator +(const cp &s)const{return {x+s.x,y+s.y};}
    cp operator -(const cp &s)const{return {x-s.x,y-s.y};}
    cp operator *(const cp &s)const{return {x*s.x-y*s.y,x*s.y+y*s.x};}
}W[18][1<<17],A[270010],B[270010],C[270010],D[270010];
int R[270010];
void FFT(const int n,cp *A){
    for(rt i=0;i<n;i++)if(i>R[i])swap(A[i],A[R[i]]);
    for(rt j=0;j<n;j+=2){
        const cp x=A[j],y=A[j+1];
        A[j]=x+y,A[j+1]=x-y;
    }
    for(rt i=2,s=1;i<n;i<<=1,s++)
    for(rt j=0;j<n;j+=i<<1)
    for(rt k=0;k<i;k+=2){
        cp x=A[j+k],y=W[s][k]*A[i+j+k];
        A[j+k]=x+y;A[i+j+k]=x-y;
        x=A[j+k+1],y=W[s][k+1]*A[i+j+k+1];
        A[j+k+1]=x+y;A[i+j+k+1]=x-y;
    }
}
int yg(int n){
    for(rt i=1;i<=n;i++){
        for(rt j=2;j*j<=n;j++)if((n-1)%j==0)
        if(ksm(i,(n-1)/j)==1)goto end;
        return i;end:;
    }
}
void subtask(){
    w[0]=1;w[1]=ksm(yg(p),(p-1)/k);
    for(rt i=2;i<k;i++)w[i]=1ll*w[i-1]*w[1]%p;mat s=zy;
    for(rt i=0;i<k;i++){
        g[i]=ksm(s*w[i]+I,L).a[x-1][y-1]*w[V(i)]%p;
    }    
    for(rt i=0;i<k+k;i++)f[i]=w[(k-V(i))%k];
    for(rt i=0;i<k/2;i++)swap(g[i],g[k-i]);

    n=k+k-1;m=k;
    int lim=1;while(lim<=n+m)lim<<=1;
    for(rt i=0;(1<<i)<lim;i++){
        W[i][0]={1,0};
        W[i][1]={cos(PI/(1<<i)),sin(PI/(1<<i))};   
        for(rt j=2;j<1<<i;j++)
        if(j%32==0)W[i][j]={cos(PI*j/(1<<i)),sin(PI*j/(1<<i))};
        else W[i][j]=W[i][j-1]*W[i][1];
    }
 
    for(rt i=0;i<=n;i++){
        A[i].x=f[i]>>15;
        A[i].y=f[i]&32767;
    }
    for(rt i=0;i<=m;i++){
        B[i].x=g[i]>>15;
        B[i].y=g[i]&32767;
    }
    for(rt i=1;i<lim;i++)R[i]=(R[i>>1]>>1)|(i&1)*(lim>>1);
    FFT(lim,A);FFT(lim,B);
    for(rt i=0;i<lim;i++){
        const int pl=(lim-1)&(lim-i);
        const cp ca={A[pl].x,-A[pl].y},cb={B[pl].x,-B[pl].y};
        const cp a=(A[i]+ca)*(cp){0.5,0},b=(A[i]-ca)*(cp){0,-0.5},
                  c=(B[i]+cb)*(cp){0.5,0},d=(B[i]-cb)*(cp){0,-0.5};
        C[pl]=a*c+a*d*(cp){0,1};D[pl]=b*c+b*d*(cp){0,1};
    }
    FFT(lim,C);FFT(lim,D);
    for(rt i=k;i<k+k;i++){
        const int vv=1ll*w[V(i-k)]*ksm(k,p-2)%p;
        ll a=C[i].x/lim+0.5,b=C[i].y/lim+0.5,c=D[i].x/lim+0.5,d=D[i].y/lim+0.5;
        a=((a%p)<<30)+(((b+c)%p)<<15)+d;a=(a%p+p)%p;
        writeln(1ll*a*vv%p);
    }
    exit(0);
}
int jc[100010],njc[100010],inv[100010],ff[100010][3],ans[100010];
int c(int x,int y){
    return 1ll*jc[x]*njc[y]%p*njc[x-y]%p;
}
void bl(){
    for(rt i=0;i<2;i++)jc[i]=njc[i]=inv[i]=1;
    for(rt i=2;i<=L;i++){
        jc[i]=1ll*jc[i-1]*i%p;
        inv[i]=1ll*inv[p%i]*(p-p/i)%p;
        njc[i]=1ll*njc[i-1]*inv[i]%p;
    }
    ff[0][x-1]=1;if(x==y)ans[0]=1;
    for(rt i=1;i<=L;i++)
    for(rt j=0;j<n;j++){
        for(rt k=0;k<n;k++)
        (ff[i][j]+=1ll*zy.a[k][j]*ff[i-1][k]%p)%=p;
        if(j==y-1)(ans[i%k]+=1ll*ff[i][j]*c(L,i)%p)%=p;
    }
    for(rt i=0;i<k;i++)writeln(ans[i]);exit(0);
}
int main(){
    cin>>n>>k>>L>>x>>y>>p;
    for(rt i=0;i<n;i++)
    for(rt j=0;j<n;j++)
    cin>>zy.a[i][j];
    if(L<=100000)bl();
    for(rt i=0;i<n;i++)I.a[i][i]=1;
    subtask();
    return 0;
}