【SDOI2015】【NTT】【dp/生成函数】序列统计

【描述】
小C有一个集合S，里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数列，数列中的每个数都属于集合S。小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数x，求所有可以生成出的，且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为，两个数列{Ai}和{Bi}不同，当且仅当至少存在一个整数i，满足Ai≠Bi。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案mod 1004535809的值就可以了。

【思路】

前置知识：
1.NTT应用拓展
这也是一道NTT好题。我们知道NTT和FFT可以快速计算如下形式的式子：
$$f[k]=\sum _{i+j=k}f[i]\ast f[j]$$
而在模数较小时，我们还可以支持快速计算如下形式的式子：
$$f[k]=\sum _{i\ast j=k}f[i]\ast f[j]$$
怎么把乘法变为加法呢？取对数。所以设g为模数的一个原根，我们知道：$g^0,g^1......g^{mod-2}$是模意义下[1,mod-1]的一个排列。我们考虑把下标i的含义变为$g^i$就可以把乘法变为加法了。
2.找原根
我们先把mod-1质因数分解。对于mod-1的每个质因数p，我们验证一个数a是否在模意义下满足$a^{\frac{mod-1}{p}}=1$即可。

下面进入正题。我们设生成函数$f(x)={\sum }_{a_iϵS}{x}^{{\log }_g{a}_i}$，设题目中要求的数列的积为need，那么我们要求的就是这个生成函数模意义下的n次方$x^{{\log }_{g}need}$的系数，注意这里的多项式乘法是循环卷积，而且是模mod-1意义下的循环卷积。当然这样的方法也可以理解为dp倍增优化，即设$f[i][j]$是使用了i个数，乘积为j的方案数，倍增即可，这里倍增时也需要使用前置知识第一条优化为$O(m\log m\log n)$，否则$O(m^2\log n)$会超时。其实这两种方法本质是一样的，代码写起来都一个样子。
代码：

#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define re register
using namespace std;
const int N=2e4+5;
const int mod=1004535809;
map<int,int>mp;
inline int red(){
    int data=0;int w=1; char ch=0;
    ch=getchar();
    while(ch!='-' && (ch<'0' || ch>'9')) ch=getchar();
    if(ch=='-') w=-1,ch=getchar();
    while(ch>='0' && ch<='9') data=(data<<3)+(data<<1)+ch-'0',ch=getchar();
    return data*w;
}
int n,m,a,b,c;
int f[N],lim=1,l=0,rev[N],g;
inline int add(const int&a,const int&b){return (a+b)>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return (a-b)<0?a-b+mod:a-b;}
inline int mul(const int&a,const int&b){return 1ll*a*b%mod;}
inline int ksm(int a,int b,const int&mod){
	int ret=1;
	while(b){if(b&1)ret=1ll*ret*a%mod;a=1ll*a*a%mod;b>>=1;}
	return ret;
}
inline void pre(){lim=1;l=0;
	while(lim<((m-1)<<1))lim<<=1,++l;
	for(int re i=0;i<lim;++i)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
} 
inline void ntt(int *f,int type){
	for(int re i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	int now=type==1?3:(mod+1)/3;
	for(int re mid=1;mid<lim;mid<<=1){
		int tmp=ksm(now,(mod-1)/(mid<<1),mod);
		for(int re j=0;j<lim;j+=mid<<1){int w=1;
			for(int re k=0;k<mid;k++,w=mul(w,tmp)){
				int x=f[j+k],y=mul(w,f[j+k+mid]);
				f[j+k]=add(x,y);f[j+k+mid]=dec(x,y);
			} 
		}
	}
	if(type==-1){
		int inv=ksm(lim,mod-2,mod);
		for(int re i=0;i<lim;++i)f[i]=mul(f[i],inv);
	}
}
struct poly{
	int ret[17000];
	friend inline poly operator*(poly a,poly b){
		static poly c;memset(c.ret,0,sizeof(c.ret));
		ntt(a.ret,1);ntt(b.ret,1);
		for(int re i=0;i<lim;++i)c.ret[i]=mul(a.ret[i],b.ret[i]);
		ntt(c.ret,-1);
		for(int re i=0;i<m-1;i++)c.ret[i]=add(c.ret[i],c.ret[i+m-1]);
		for(int re i=m-1;i<lim;++i)c.ret[i]=0;
		return c;
	}
	friend inline poly operator^(poly a,int b){
		static poly c;memset(c.ret,0,sizeof(c.ret));c.ret[0]=1;
		for(;b;b>>=1,a=a*a)if(b&1)c=c*a;
		return c;
	}
}A; 
int p[N],cnt=0;
inline void div(int m){
	for(int re i=2;i*i<=m;++i){
		if(m%i==0){
			p[++cnt]=i;
			while(m%i==0)m/=i;
		}
	}if(m>1)p[++cnt]=m;
}
inline void getg(){bool flag;div(m-1);
	while(g=rand()%(m-1)+1){flag=0;
		for(int re i=1;i<=cnt;i++){
			if(ksm(g,(m-1)/p[i],m)==1){flag=1;break;}
		}if(flag)continue;
		else return;
	}
}
int x;
int main(){
	n=red();m=red();x=red();getg();pre();
	for(int re p=g,i=1;i<m-1;(p*=g)%=m,++i)mp[p]=i;
	for(int re i=red();i;--i)
		if(a=red())A.ret[mp[a]]=1;
	A=A^n;cout<<A.ret[mp[x]];
}