给出(n,m,x,S)，其中(Ssubseteq [0,m))，问有多少个长度为(n)的数列(a)使得(a_iin S)，并且数列中所有元素的乘积mod (m)为(x)。
答案对(1004535809=479 imes 2^{21}+1)取模。(m)为质数。(nle 10^9,m< 8000)。
分析
在实数意义下，一个经典的处理乘法的方法是通过取对数转化成加法，例如在计算最小乘积生成树时，我们对边权取对数，转化成一般的最小生成树问题。在模意义下也有类似的巧妙方法，利用模数的原根。
可以证明(1,2,4,p,2p,p^n)都有原根，而且只有这些数有原根（(p)表示奇素数）。原根(g)满足(g^xequiv 1)成立当且仅当(x=p-1)。这就是说，用(g)的幂可以在模意义下表示出([1,p))的所有数。也就是说，如果模数有原根，那么模意义下的乘法可以转化成原根的指数的加法。
令(x=g^b)，那么问题就变成了，在数集中，选出(n)个数，一个数可以选多次，有多少种方法能够让选出来的数的和为(b)。
运用生成函数(f)，那么问题就变成了(f^n(x))，注意这里其实是一个(mod m-1)的循环卷积，最后累加回来即可。(n)比较大用快速幂。
这题的关键在于利用原根把乘法转化成加法。
代码
写这个题主要是想试一下任意模数fft的做法，所以谢了两份代码，第一份是ntt的，第二个是拆模数的fft。
ntt
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long giant;
const int q=479<<21|1;
const int g=3;
int read() {
int x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const int maxj=16;
const int maxn=1<<maxj|1;
int a[maxn],b[maxn],c[maxn],mg,M,xj,rev[maxn],n,m,d[maxn];
int Wn[maxj][2];
int sub(int x,int y) {
int ret=((giant)x-y)%q;
if (ret<0) ret+=q;
return ret;
}
int Plus(int x,int y) {
return ((giant)x+y)%q;
}
int multi(int x,int y) {
return (giant)x*y%q;
}
int mi(int x,int y) {
int ret=1;
for (;y;y>>=1,x=multi(x,x)) if (y&1) ret=multi(ret,x);
return ret;
}
int inv(int x) {
return mi(x,q-2);
}
void ntt(int a[],int op=0) {
for (int i=0;i<M;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int the=1,i=2;i<=M;i<<=1,++the) {
for (int j=0;j<M;j+=i) {
for (int k=j,w=1;k<j+i/2;++k,w=multi(w,Wn[the][op])) {
int u=a[k],v=multi(a[k+i/2],w);
a[k]=Plus(u,v),a[k+i/2]=sub(u,v);
}
}
}
if (op) {
int iv=inv(M);
for (int i=0;i<M;++i) a[i]=multi(a[i],iv);
}
}
void self(int a[]) {
ntt(a);
for (int i=0;i<M;++i) a[i]=multi(a[i],a[i]);
ntt(a,1);
}
void mul(int c[],int a[]) {
ntt(a);
ntt(c);
for (int i=0;i<M;++i) c[i]=multi(c[i],a[i]);
ntt(a,1);
ntt(c,1);
}
void back(int a[]) {
memset(d,0,sizeof d);
for (int i=0;i<M;++i) d[i%(m-1)]=Plus(d[i%(m-1)],a[i]);
memcpy(a,d,(sizeof d[0])*M);
}
void ans(int a[],int len,int y) {
memset(c,0,sizeof c);
c[0]=1;
for (xj=0,M=1;M<=(len<<1);++xj,M<<=1);
for (int i=0;i<M;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
while (y) {
if (y&1) {
mul(c,a);
back(c);
}
y>>=1;
self(a);
back(a);
}
}
int getg(int m) {
for (int i=2;i<m;++i) {
int tmp=i;
bool flag=true;
for (int j=1;j<m;++j,tmp=tmp*i%m) if (tmp==1 && j!=m-1) {
flag=false;
break;
}
if (flag) return i;
}
return 0;
}
int S[maxn];
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
freopen("my.out","w",stdout);
#endif
Wn[0][0]=Wn[0][1]=1;
for (int i=1;i<maxj;++i) Wn[i][1]=inv(Wn[i][0]=mi(g,(q-1)>>i));
n=read(),m=read();
int x=read(),s=read();
mg=getg(m);
for (int tmp=1,i=0;i<m-1;++i,tmp=tmp*mg%m) b[tmp]=i;
for (int i=1;i<=s;++i) S[i]=read();
for (int i=1;i<=s;++i) if (S[i]) ++a[b[S[i]]];
ans(a,m,n);
printf("%d
",c[b[x]]);
return 0;
}
fft
任意模数fft的做法是，设模数为(p)，把多项式的系数拆成(asqrt p+b)，分成两个多项式((A,B))。
假设(X=Asqrt p+B,Y=Csqrt p+D)，那么进行分步计算：
[X imes Y=A imes C*(sqrt p*sqrt p) +(B imes C+A imes D)*sqrt p+B imes D mod p
]
分开取模，这样就不会受到精度问题影响了。要开long double。
#include<cstdio>
#include<cctype>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long double lb;
int read() {
int x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
typedef long long giant;
const int q=479<<21|1;
const int sq=31695;
const int sqb=(giant)sq*sq%q;
const lb pi=acos(-1.0);
struct C {
lb r,i;
C (lb r=0,lb i=0):r(r),i(i) {}
};
C operator * (C a,C b) {return C(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
C operator + (C a,C b) {return C(a.r+b.r,a.i+b.i);}
C operator - (C a,C b) {return C(a.r-b.r,a.i-b.i);}
C operator / (C a,lb x) {return C(a.r/x,a.i/x);}
const int maxj=16;
const int maxn=1<<maxj|1;
int M,a[maxn],c[maxn],xj,rev[maxn],b[maxn];
struct Array {
C a[maxn];
C& operator [] (int x) {return a[x];}
void back(int m) {
for (int i=m-1;i<M;++i) a[i%(m-1)]=a[i%(m-1)]+a[i],a[i]=0;
}
void init(int b[]) {
for (int i=0;i<M;++i) a[i]=b[i];
}
void one() {
for (int i=0;i<M;++i) a[i]=0;
a[0]=1;
}
void fft(int op=1) {
for (int i=0;i<M;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=2;i<=M;i<<=1) {
C wn(cos(2*pi/i),op*sin(2*pi/i));
for (int j=0;j<M;j+=i) {
C w(1);
for (int k=j;k<j+i/2;++k,w=w*wn) {
C u=a[k],v=a[k+i/2]*w;
a[k]=u+v,a[k+i/2]=u-v;
}
}
}
if (op==-1) {
for (int i=0;i<M;++i) a[i]=a[i]/M;
}
}
} A,B,D,E,F,G,H,TA,TC;
Array operator * (Array &a,Array &b) {
Array ret;
a.fft(),b.fft();
for (int i=0;i<M;++i) ret[i]=a[i]*b[i];
a.fft(-1),b.fft(-1),ret.fft(-1);
return ret;
}
Array operator + (Array &a,Array &b) {
Array ret;
for (int i=0;i<M;++i) ret[i]=a[i]+b[i];
return ret;
}
void operator *= (Array &a,Array &b) {
a.fft(),b.fft();
for (int i=0;i<M;++i) a[i]=a[i]*b[i];
a.fft(-1),b.fft(-1);
}
void split(Array &a,Array &c,Array &d) {
for (int i=0;i<M;++i) {
giant x=(giant)(a[i].r+0.5);
giant y=x%sq;
giant z=x/sq;
c[i]=z;
d[i]=y;
}
}
int multi(giant x,giant y) {
return (giant)x*y%q;
}
int Plus(giant x,giant y) {
return ((giant)x+y)%q;
}
Array mul(Array &a,Array &b) {
Array ret;
split(a,A,B);
split(b,D,E);
F=A*D;
G=A*E;
H=B*D;
G=G+H;
H=B*E;
for (int i=0;i<M;++i) {
giant x=(giant)(F[i].r+0.5);
giant y=(giant)(G[i].r+0.5);
giant z=(giant)(H[i].r+0.5);
x=multi(x,sqb);
y=multi(y,sq);
x=Plus(x,Plus(y,z));
ret[i]=x;
}
return ret;
}
void ans(int a[],int m,int y) {
for (xj=0,M=1;M<=(m<<1);M<<=1,++xj);
for (int i=0;i<M;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
TA.init(a);
TC.one();
while (y) {
if (y&1) TC=mul(TC,TA),TC.back(m);
y>>=1;
TA=mul(TA,TA),TA.back(m);
}
for (int i=0;i<M;++i) c[i]=((giant)(TC[i].r+0.5))%q;
}
int getm(int m) {
for (int i=2;i<m;++i) {
bool flag=true;
for (int tmp=i,j=1;j<m;++j,tmp=tmp*i%m) if (tmp==1 && j!=m-1) {
flag=false;
break;
}
if (flag) return i;
}
return 0;
}
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
freopen("my.out","w",stdout);
#endif
int n=read(),m=read(),x=read(),s=read();
int mg=getm(m);
for (int tmp=1,i=0;i<m-1;++i,tmp=tmp*mg%m) b[tmp]=i;
for (int i=1,the;i<=s;++i) if ((the=read())!=0) ++a[b[the]];
ans(a,m,n);
printf("%d
",c[b[x]]);
return 0;
}