DP一下,设$f_{i,j}$表示生成$i$个数且乘积$\%M=j$的方案数,则$f_{i+1,l}=\sum\limits_{jk\%M=l}[k\in S]f_{i,j}$
我们很不希望DP式中下标的位置出现乘法,因为这样不好转移,考虑把乘法换成加法
因为模数$M$是质数,所以它有原根,于是对于每个$1\leq x\leq M-1$都有唯一的$0\leq x'\leq M-2$使得$x=g^{x'}$,所以DP式可以写成$f_{i+1,g^{l'}}=\sum\limits_{g^{j'+k'}\%M=g^{l'}}[g^{k'}\in S]f_{i,g^{j'}}$,令$h_{i,j}=f_{i,g^j},S'=\{x|g^x\in S,0\leq x\leq M-2\}$,则$h_{i+1,l}=\sum_{(j+k)\%(M-1)=l}\limits[k\in S']h_{i,j}$,这是卷积的形式($j+k\geq M-1$的项最后加到低次的项中),所以我们可以用FFT做到在$O(n\log_2n)$的时间内从$h_i$转移到$h_{i+1}$
$n$很大,但是每次转移都是一样的,所以只需要快速幂就好了,总时间复杂度$O(M\log_2M\log_2n)$
注意题目给的$S$中是有可能出现$0$的,忽略就好
#include<stdio.h>
#include<string.h>
const int mod=1004535809;
typedef long long ll;
int mul(int a,int b,int p=mod){return a*(ll)b%p;}
int ad(int a,int b){return(a+b)%mod;}
int de(int a,int b){return(a-b)%mod;}
void swap(int&a,int&b){a^=b^=a^=b;}
int rev[40010],N,iN,m;
int pow(int a,int b){int s=1;while(b){if(b&1)s=mul(s,a);a=mul(a,a);b>>=1;}return s;
}
void pre(int n){int i,k;for(N=1,k=0;N<n;N<<=1)k++;for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));iN=pow(N,mod-2);
}
void ntt(int*a,int on){int i,j,k,t,w,wn;for(i=0;i<N;i++){if(i<rev[i])swap(a[i],a[rev[i]]);}for(i=2;i<=N;i<<=1){wn=pow(3,(on==1)?(mod-1)/i:(mod-1-(mod-1)/i));for(j=0;j<N;j+=i){w=1;for(k=0;k<i>>1;k++){t=mul(w,a[i/2+j+k]);a[i/2+j+k]=de(a[j+k],t);a[j+k]=ad(a[j+k],t);w=mul(w,wn);}}}if(on==-1){for(i=0;i<N;i++)a[i]=mul(a[i],iN);}
}
struct poly{int n,a[40010];poly(){n=0;memset(a,0,sizeof(a));}
};
poly operator*(poly x,poly y){int i;poly c;pre(x.n+y.n+1);ntt(x.a,1);ntt(y.a,1);for(i=0;i<N;i++)c.a[i]=mul(x.a[i],y.a[i]);ntt(c.a,-1);for(i=m-1;i<N;i++){c.a[i%(m-1)]=ad(c.a[i%(m-1)],c.a[i]);c.a[i]=0;}c.n=m-2;return c;
}
int getg(int p){int i,j,t;for(i=2;i<p;i++){t=1;for(j=1;j<p-1;j++){t=mul(t,i,p);if(t==1)break;}if(j==p-1)break;}return i;
}
int s[8010],r[8010];
poly pow(poly a,int b){poly s;s.a[0]=1;while(b){if(b&1)s=s*a;a=a*a;b>>=1;}return s;
}
int main(){int n,x,ss,i,g,b;poly a;scanf("%d%d%d%d",&n,&m,&x,&ss);for(i=1;i<=ss;i++)scanf("%d",s+i);g=getg(m);b=1;for(i=0;i<m-1;i++){r[b]=i;b=mul(b,g,m);}for(i=1;i<=ss;i++){if(s[i])a.a[r[s[i]]]=1;}a.n=m-2;while(a.a[a.n]==0)a.n--;a=pow(a,n);printf("%d",(a.a[r[x]]+mod)%mod);
}