[BZOJ3992] [SDOI2015] 序列统计 [NTT&原根&指标&多项式快速幂]

Link
https://www.lydsy.com/JudgeOnline/problem.php?id=3992

---

题意
给定一个集合 $S$ ， $S$ 中元素均为小于 $M$ 的非负整数。（$M$ 为质数）
问：用集合中的数能够构造出多少个不同的数列（允许多次使用同一个数），
满足数列中所有数的乘积 $\equiv x(\mathrm{m}\mathrm{o}\mathrm{d}M)$？ $x\in [1,M-1]$
答案模 $1004535809$ 。

---

设计：填了前 $i$ 个元素，乘积 $\equiv j(\mathrm{m}\mathrm{o}\mathrm{d}M)$ 的方案数 $f(i,j)$
$ans=f(|S|,x)(\mathrm{m}\mathrm{o}\mathrm{d}p)$
转移 $f(i+1,j)=\sum _k\sum _{t=1}^{|S|}f(i,k)[S_t\equiv k^{-1}j(\mathrm{m}\mathrm{o}\mathrm{d}M)]$

---

注意到 $M$ 是质数就算了， $M\ge 3$ ……干什么用的？奇素数？？
那先不管具体怎么样，我们思考用用原根
$f(i,j)$ 表示填前 $i$ 个元素，乘积 $\equiv g^j(\mathrm{m}\mathrm{o}\mathrm{d}M)$
$f(i+1,j)\equiv \sum _{k=0}^j\sum _{t=1}^{|S|}f(i,k)[S_t\equiv g^{j-k}(\mathrm{m}\mathrm{o}\mathrm{d}M)](\mathrm{m}\mathrm{o}\mathrm{d}p)$

---

设 $G(x)=\sum _{t=1}^{|S|}[S_t\equiv g^x(\mathrm{m}\mathrm{o}\mathrm{d}M)](\mathrm{m}\mathrm{o}\mathrm{d}p)$
（实际上 G ( x ) ∈ { 0 , 1 } G(x) \in \{0,1\} G(x)∈{0,1} ）
但是这样的话 $g^0\equiv g^{M-1}(\mathrm{m}\mathrm{o}\mathrm{d}M)$
所以我们应该设 $G(x)=\sum _{i=1}^{|S|}[I_g(S_t)=x(\mathrm{m}\mathrm{o}\mathrm{d}\varphi (M))](\mathrm{m}\mathrm{o}\mathrm{d}p)$

---

$f(i+1,j)=\sum _{k=0}^{j}f(i,k)G(j-k)(\mathrm{m}\mathrm{o}\mathrm{d}p)$
即 $f(i+1)\equiv f(i)\otimes G(\mathrm{m}\mathrm{o}\mathrm{d}p)$
那么 $f(N)\equiv G^N(\mathrm{m}\mathrm{o}\mathrm{d}p)$
注意在卷 $G$ 卷卷的时候要不停把下标馍馍 $(\mathrm{m}\mathrm{o}\mathrm{d}\varphi (M))$

---

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<cctype>
#include<ctime>
using namespace std;
#define getchar() (frS==frT&&(frT=(frS=frBB)+fread(frBB,1,1<<12,stdin),frS==frT)?EOF:*frS++)
char frBB[1<<12], *frS=frBB, *frT=frBB;
inline void read(int& x)
{
	x = 0; char ch = getchar(); bool w = 0;
	while (!isdigit(ch)) w |= (ch == '-'), ch = getchar();
	while (isdigit(ch)) x = x * 10 + (ch ^ 48), ch = getchar();
	w ? (x = -x) : 0;
}
const int MAXN = 16384+16;
const long long p = 1004535809;
inline long long qpow(long long a, long long b, const long long& p)
{
	long long ret = 1;
	a %= p;
	while(b)
	{
		if (b & 1)
		{
			ret *= a;
			ret %= p;
		}
		a *= a;
		a %= p;
		b >>= 1;
	}
	return ret;
}
const long long g = 3;
int n, m, x, s, Lim = 1, Log;
long long q;
long long Wn[MAXN][2];
int Rev[MAXN];
inline long long Adjust(long long x)
{
	x = (x % p) + p;
	return (x >= p) ? (x - p) : x;
}
void NTT(long long *a, const bool& Type = 0)
{
	for (register int i = 0; i < Lim; ++i) if (Rev[i] > i) swap(a[Rev[i]], a[i]);
	register long long x;
	for (register int Len, Mid = 1, qwq; Mid < Lim; Mid <<= 1)
	{
		Len = Mid << 1;
		qwq = Lim / Len;
		for (register int Pos = 0; Pos < Lim; Pos += Len)
		{
			for (register int Sub = 0; Sub < Mid; ++Sub)
			{
				x = Wn[qwq * Sub][Type] * a[Pos + Mid + Sub] % p;
				a[Pos + Mid + Sub] = Adjust(a[Pos + Sub] - x);
				a[Pos + Sub] = (a[Pos + Sub] + x) % p;
			}
		}
	}
}
int gM;
int PMStack[8005],  Ind[8005];
inline void PMSplit(int x)
{
	for (register int i = 2; i * i <= x; ++i)
	{
		if (x % i == 0)
		{
			PMStack[++PMStack[0]] = i;
			while (x % i == 0) x /= i;
		}
	}
}
inline int GetPrimeRoot(int x, int phi)
{
	PMSplit(phi);
	register bool flag;
	for (register int i = 2; i < x; ++i)
	{
		flag = 0;
		for (register int j = 1; j <= PMStack[0]; ++j)
		{
			if (qpow(i, phi/PMStack[j], m) == 1)
			{
				flag = 1;
				break;
			}
		}
		if (!flag) return i;
	}
}
long long a[MAXN];
long long ta[MAXN], tb[MAXN], inv;
inline void PolyMul(long long *x, long long *y)
{
	for (register int i = 0; i < Lim; ++i) ta[i] = x[i] % p;
	for (register int i = 0; i < Lim; ++i) tb[i] = y[i] % p;
	NTT(ta); NTT(tb);
	for (register int i = 0; i < Lim; ++i) ta[i] = ta[i] * tb[i] % p;
	NTT(ta, 1);
	for (register int i = 0; i < Lim; ++i) ta[i] = ta[i] * inv % p;
	for (register int i = 0; i < m-1; ++i) x[i] = (ta[i] + ta[i+m-1]) % p;
	//同样要小心自乘的情况
}
long long Ans[MAXN];
void PolyPow(long long *a, int Index)
{
	Ans[0] = 1;
	while (Index)
	{
		if (Index & 1)
		{
			PolyMul(Ans, a);
		}
		PolyMul(a, a);
		Index >>= 1;
	}
}
int main()
{
	read(n), read(m), read(x), read(s);
	gM = GetPrimeRoot(m, m-1);
	for (register int i = 0, cur = 1; i < m - 1; ++i) Ind[cur] = i, cur *= gM, cur %= m;
	for (register int t, i = 1; i <= s; ++i) read(t), t ? (a[Ind[t]] = 1) : 0;
	int mm = m << 1;
	while (Lim <= mm) Lim <<= 1, ++Log;
	for (register int i = 0; i < Lim; ++i) Rev[i] = (Rev[i>>1]>>1)|((i&1)<<(Log-1));
	inv = qpow(Lim, p-2, p);
	q = (p - 1) / Lim;
	register long long fafa = qpow(g, q, p);
	Wn[0][0] = Wn[Lim][0] = Wn[0][1] = Wn[Lim][1] = 1;
	for (register int i = 1; i < Lim; ++i) Wn[i][0] = Wn[i-1][0] * fafa % p;
	for (register int i = 1; i < Lim; ++i) Wn[i][1] = Wn[Lim-i][0];
	PolyPow(a, n);
	printf("%lld", Ans[Ind[x]]);
	return 0;
}