[SDOI2015]序列统计

题意

$a_i\in S,$求$\prod_{i=1}^na_i\equiv x\mod m$的方案数

$m\in P,m\le8000,n\le10^9,x\in[1,m-1]$

---

题解

设$f[i][j]$为$\prod_{k=1}^ia_k\equiv j\mod m$的方案数

那么有

$$f[i+1][j]=\sum_{ab\equiv j\mod m}f[i][a]f[i][b]$$

取个对数$\ln a+\ln b=\ln j$这样似乎就是卷积了

考虑到$m$是一个质数$,$然后注意到质数的原根$g$满足

$g^0,g^1,\ldots g^{m-2}$在模$m$的意义下可以取遍$[1,m-1]$所有数

设$g^A\equiv a\mod m,$因为$A$和$a$是一一对应的$,$不妨用$A$来替代$a$

那么有

$$f[i+1][J]=\sum_{g^Ag^B\equiv g^J\mod m}f[i][A]f[i][B]$$

$$\Rightarrow f[i+1][J]=\sum_{g^{A+B}\equiv g^J\mod m}f[i][A]f[i][B]$$

根据费马小定理$a^b\equiv a^{b\mod P-1}\mod P$

$$\Rightarrow f[i+1][J]=\sum_{g^{A+B\mod m-1}=g^J}f[i][A]f[i][B]$$

$$\Rightarrow f[i+1][J]=\sum_{A+B\equiv J\mod m-1}f[i][A]f[i][B]$$

设$g[J]=\sum\limits_{A+B=J}f[i][A]f[i][B]$

$$\Rightarrow f[i][J]=g[J]+g[J+m-1]$$

因为这是个卷积并且每一层的转移都是一样的$,$所以可以直接用多项式快速幂算

注意集合中有$0,$直接判掉即可

#include<bits/stdc++.h>
#define fp(i,a,b) for(register int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(register int i=a,I=b-1;i>I;--i)
#define go(u) for(register int i=fi[u],v=e[i].to;i;v=e[i=e[i].nx].to)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
char ss[1<<17],*A=ss,*B=ss;
inline char gc(){return A==B&&(B=(A=ss)+fread(ss,1,1<<17,stdin),A==B)?-1:*A++;}
template<class T>inline void sd(T&x){
    char c;T y=1;while(c=gc(),(c<48||57<c)&&c!=-1)if(c==45)y=-1;x=c-48;
    while(c=gc(),47<c&&c<58)x=x*10+c-48;x*=y;
}
char sr[1<<21],z[20];int C=-1,Z;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
template<class T>inline void we(T x){
    if(C>1<<20)Ot();if(x<0)sr[++C]=45,x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]='\n';
}
const int N=1e5+5,P=1004535809,G=3;
typedef int arr[N];
typedef long long ll;
int n,m,M,K,X,pw;arr a,b,c,pos,R,wn;
inline int pls(int a,int b){return a+=b,a<P?a:a-P;}
inline int sub(int a,int b){return a-=b,a<0?a+P:a;}
inline int fpm(int a,int b,const int P){int x=1;for(;b;b>>=1,a=(ll)a*a%P)if(b&1)x=(ll)x*a%P;return x;}
inline void init(){for(int i=0,j=0;i<K;++i){R[i]=j;for(int l=K>>1;(j^=l)<l;l>>=1);}}
inline void dft(int*a,int ty=1){
    fp(i,0,K-1)if(i>R[i])swap(a[i],a[R[i]]);
    for(int i=2;i<=K;i<<=1){
        int d=i>>1,o;wn[1]=ty==1?fpm(G,(P-1)/i,P):fpm(G,P-1-(P-1)/i,P);
        fp(i,2,d-1)wn[i]=(ll)wn[i-1]*wn[1]%P;
        for(int*x=a;x!=a+K;x+=i)fp(k,0,d-1)
            o=(ll)x[k+d]*wn[k]%P,x[k+d]=sub(x[k],o),x[k]=pls(x[k],o);
    }if(ty^1){
        const int inv=fpm(K,P-2,P);
        fp(i,0,K-1)a[i]=(ll)a[i]*inv%P;
    }
}
inline void calc(int*a){fp(i,0,m-1)a[i]=pls(a[i],a[i+m]);}
inline void Sqr(int*a){
    dft(a);fp(i,0,K-1)a[i]=(ll)a[i]*a[i]%P;
    dft(a,-1);calc(a);fp(i,m,K-1)a[i]=0;
}
inline void Mul(int*a,int*b){
    fp(i,0,m)c[i]=b[i];dft(a),dft(c);
    fp(i,0,K-1)a[i]=(ll)a[i]*c[i]%P;
    dft(a,-1);calc(a);fp(i,m,K-1)a[i]=c[i]=0;
}
inline void Fpm(int*x,int*a,int b){for(;b;b>>=1,Sqr(a))if(b&1)Mul(x,a);}
inline int Gg(const int P){
    static int s[N],t=0;bool fg;int Phi=P-1;
    fp(i,2,sqrt(Phi))if(Phi%i==0)s[++t]=i,Phi/i!=i?s[++t]=Phi/i:0;
    fp(i,2,Phi){fg=1;
        fp(j,1,t)if(fpm(i,s[j],P)==1){fg=0;break;}
        if(fg)return i;
    }return -1;
}
int main(){
    #ifndef ONLINE_JUDGE
        file("s");
    #endif
    sd(pw),sd(m),sd(X),sd(n);wn[0]=1;
    K=1;while(K<m)K<<=1;K<<=1;init();
    int p=1,g=Gg(m);fp(i,1,m-2)pos[(p*=g)%=m]=i;
    while(n--)sd(p),p?++b[pos[p]]:0;M=2*(--m);
    a[0]=1;Fpm(a,b,pw);we(a[pos[X]]);
return Ot(),0;
}