BZOJ3992: [SDOI2015]序列统计

https://www.lydsy.com/JudgeOnline/problem.php?id=3992

有a*bΞx(mod m),则ord(a)+ord(b)Ξord(x)(mod φ(m))

对Σx^ord(s_i)循环卷积FFT即可

#include<cstdio>
#include<algorithm>
#define re register 
#define rep(i,s,t) for(re int i=s;i<=t;++i)
#define Rep(i,s,t) for(re int i=s;i<t;++i)
using namespace std;
const int N=4e5+11,mod=1004535809;
int n,m,x,s,g,sum,d,lg2;
int q[N],ord[N],a[N],p[N],res[N];
inline void inc(int &a,int b){
    a+=b;
    if(a>=mod)a-=mod;
    if(a<0)a+=mod;
}
inline int fp(int a,int b,int m=mod){
    if(b<0)b+=m-1;
    if(b>=m-1)b-=m-1;
    int res=1;
    for(;b;b>>=1,a=1ll*a*a%m)
        if(b&1)
            res=1ll*res*a%m;
    return res;
}
inline int findg(int x){
    int t=0;
    for(re int i=2;i*i<=x-1;++i)
        if((x-1)%i==0){
            q[++t]=i;
            if(i*i!=x-1)
                q[++t]=(x-1)/i;
        }
    rep(i,2,(1<<30)){
        rep(j,1,t)
            if(fp(i,q[j],x)==1)
                goto end;
        return i;
        end:;
    }
}
inline void dft(int *a,int n,int f){
    Rep(i,0,n)
        if(i<p[i])
            swap(a[i],a[p[i]]);
    for(re int i=1,x,y;i<n;i<<=1){
        for(re int w,j=0,wn=fp(3,(mod-1)/(i<<1)*f);w=1,j<n;j+=(i<<1)){
            for(re int k=j;k<i+j;++k,w=1ll*w*wn%mod){
                x=a[k],y=1ll*w*a[i+k]%mod;
                a[k]=(x+y)%mod;
                a[i+k]=(x-y+mod)%mod;
            }
        }
    }
    if(f==-1){
        int inv=fp(n,mod-2);
        Rep(i,0,d)
            a[i]=1ll*a[i]*inv%mod;
    }
}
int main(){
    //freopen("test.in","r",stdin);
    scanf("%d%d%d%d",&n,&m,&x,&s);
    g=findg(m),sum=1;
    //printf("%d\n",g);
    rep(j,0,m-1){
        ord[sum]=j;
        sum=1ll*sum*g%m;
    }
    for(int x;s--;){
        scanf("%d",&x);
        if(x)++a[ord[x]];
    }
    for(d=1,lg2=-1;d<=m+m;d<<=1,++lg2);
    Rep(i,0,d)p[i]=(p[i>>1]>>1)^((i&1)<<lg2);
    Rep(i,0,d)res[i]=a[i];
    for(--n;n;n>>=1){
        dft(a,d,1);
        if(n&1){
            dft(res,d,1);
            Rep(i,0,d)res[i]=1ll*res[i]*a[i]%mod;
            dft(res,d,-1);
            Rep(i,0,m-1)inc(res[i],res[i+m-1]);
            Rep(i,m-1,d)res[i]=0;
        }
        Rep(i,0,d)a[i]=1ll*a[i]*a[i]%mod;
        dft(a,d,-1);
        Rep(i,0,m-1)inc(a[i],a[i+m-1]);
        Rep(i,m-1,d)a[i]=0;
    }
    printf("%d\n",res[ord[x]]);
    return 0;
}