Description：

求用所有小于 $M$ 的非负整数构成长为 $N$ 的数列乘积 $\text{mod }M=x$ 的数列的个数 $\text{mod }1004535809$ 。

$1\leqslant N\leqslant10^9,3\leqslant M\leqslant 8000,M$ 为质数

Solution：

首先乘积这个东西不是很可做，如果不是模意义下我们可以两边取对数，然后统计和，而模意义下我们就可以使用原根来求模意义下的离散对数，由于 $M$ 是质数所以一定存在原根，且原根一定是 $\varphi(M)=M-1$ 的约数，证明： $$ \begin{align} &由欧拉定理可知：g^{\varphi(m)}\equiv1(\text{mod M})\\ &用反证法：设g^b\equiv1(\text{mod M})\\ &设b的最小值为c_1，g^{c_1}\equiv(\text{mod M})\\ &设c_2=\varphi(\text{M})-c_1，那么g^{c_2}\equiv1(\text{mod M})\\ &c_1\not|c_2，否则就会有c_1|\varphi(\text{M})\\ &因为c_2>c_1，所以设k=c_2-c_1，一定有g^k\equiv1(\text{mod M})\\ &而k<c_1，与原命题矛盾 \end{align} $$ 参考

所以我们只要枚举 $M-1$ 的所有约数，然后暴力判断对于所有 $i\in[1,\varphi(\text{M}))$ ， $g^i$ 是否都不为 $1$ 即可。

这样每个 $k$ 都会对应一个 $g^i$ ，那么我们就成功地把乘积变成了求和，

然后问题就变成了给出很多数他们的和是某个值的方案数，生成函数还是比较显然的，答案是某个多项式的 $N$ 次方，于是多项式 $+$ 快速幂即可，注意处理取模，这个可以在每次计算完一次后从高往低加。

注意所有值都是在指数上计算的，每个值还有模数都会改变，当用这个方法的时候要特别小心。

还要注意必须要跳过 $0$ 。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

ll n,m,x,s;

#define MAXN 32010

ll val[MAXN];

ll gcd(ll a,ll b){return (b == 0 ? a : gcd(b,a % b));}

bool check(ll k)

{

	if(gcd(k,m) != 1)return false;

	for(ll i = 1,v = k;i < m - 1;++i,v = v * k % m)

	{

		if(v == 1)return false;

	}

	return true;

}

ll G()

{

	for(ll i = 1;i <= m - 1;++i)

	{

		if(check(i))return i;

	}

}

ll my_lg[MAXN],my_pow[MAXN];

const ll P = 1004535809;

int rev[MAXN];

ll ww[MAXN << 1],*g = ww + MAXN;

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % P;

		a = a * a % P;

		b = b >> 1;

	}

	return res;

}

ll res[MAXN],a[MAXN];

ll inv;

void NTT(ll f[],int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		ll wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			ll w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				ll u = f[k],t = w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		for(int i = 0;i < l;++i)

		{

			f[i] = f[i] * inv % P;

		}

	}

	return;

}

void KSM(int len)

{

	res[0] = 1;

	while(n > 0)

	{

		if(n & 1)

		{

			NTT(res,len,1);NTT(a,len,1);

			for(int i = 0;i < len;++i)res[i] = res[i] * a[i] % P;

			NTT(res,len,-1);NTT(a,len,-1);

			for(int i = len - 1;i >= m - 1;--i)

			{

				res[i - m + 1] = (res[i - m + 1] + res[i]) % P;

				res[i] = 0;

			}

		}

		NTT(a,len,1);

		for(int i = 0;i < len;++i)a[i] = a[i] * a[i] % P;

		NTT(a,len,-1);

		for(int i = len - 1;i >= m - 1;--i)

		{

			a[i - m + 1] = (a[i - m + 1] + a[i]) % P;

			a[i] = 0;

		}

		n = n >> 1;

	}

	return;

}

int main()

{

	scanf("%lld%lld%lld%lld",&n,&m,&x,&s);

	for(int i = 1;i <= s;++i)scanf("%lld",&val[i]);

	ll m_g = G();

	for(ll i = 0,cur = 1;i < m - 1;++i,cur = cur * m_g % m)

	{

		my_lg[cur] = i;my_pow[i] = cur;

	}

	int l = 0,len = 1;

	for(;len <= 2 * m;++l,len = len << 1);

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	for(int i = 1;i <= s;++i)

	{

		if(val[i] == 0)continue;

		a[my_lg[val[i]]] = 1;

	}

	inv = power(len,P - 2);

	KSM(len);

	printf("%lld\n",res[my_lg[x]]);

	return 0;

}