【BZOJ3992】[SDOI2015] 序列统计（原根+NTT+倍增）

点此看题面

大致题意： 给定一个所有元素都为小于(m)的非负整数的无重复元素的集合(S)。求有多少个不同的长度为(n)的序列，满足其中每个元素都属于集合(S)，且所有元素之积模(m)的值为给定值(x)。

前言

亡在了数学上。。。

前面的式子我都推出来了啊，乘法转化成加法的套路也不是没见过，可不知道原根这种用法我真的还是太弱了。

说起来原根这东西我也挺熟悉的——写(NTT)每次都要用的啊，可我这种写(FFT/NTT)从没搞懂过原理只会背板子的蒟蒻，根本未曾想过原根究竟是个啥玩意。

现在是明白了，难怪(NTT)必须要用原根，原来就是利用了原根的性质。

原根

什么是原根？

简单的说，对于一个数(p)，若(g)满足(g^0 mod p,g^1 mod p,g^2 mod p,...,g^{p-2} mod p)恰好是(1sim p-1)的一个排列，那么(g)就是(p)的原根。

那么如何求原根呢？很暴力，枚举+验证，只不过验证是需要点技巧的。

我们可以求出(phi(p))（此题为质数，因此(phi(p)=p-1)）的所有质因数(P_i)，若存在(P_i)使得(x^{frac{phi(p)}{P_i}}equiv1(mod p))，则(x)不是(p)的原根，否则(x)就是(p)的原根。

那么原根在这道题中具体有什么用呢？

等到后面就知道了。

大致想法

如果我们设(g_{i,j})为长度为(i)的序列中元素之积模(m)为(j)的方案数，则显然有转移：

[g_{x+y,a}=sum_{b imes cequiv a(mod m)}g_{x,b} imes g_{y,c} ]

显而易见这是一个多项式乘法，但转移条件是(b imes c)十分麻烦，因此我们就想到把乘法转化为加法，然后就能用(NTT)了。

而怎么把乘法转化为加法呢？

当然是用对数啦！

然后我们的原根就闪亮登场了，正因为它的性质，(log_gx=y)中(x,y)是一一对应的。

例如在上面的式子中，我们分别设(w=log_ga,u=log_gb,v=log_gc)，就可以得到：

[g_{x+y,w}=sum_{u+vequiv w(mod m-1)}g_{x,u} imes g_{y,v} ]

而加法的取模只要在(NTT)做完后把后面的给加到前面去就好了，不用多说。

这样一来，我们考虑倍增，用(f_{i,j})表示长度为(2^i)的序列中元素之积模(m)为(j)的方案数，然后就可以轻松搞定了。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define LN 30
#define M 8000
#define X 1004535809
using namespace std;
int n,m,tar,k,GR,a[M+5],cnt,v[M+5],LG[M+5],P[LN+5][M+5],ans[M+5];
I int Qpow(RI x,RI y,CI p) {RI t=1;W(y) y&1&&(t=1LL*t*x%p),x=1LL*x*x%p,y>>=1;return t;}
I bool IsP(CI x) {for(RI i=2;1LL*i*i<=x;++i) if(!(x%i)) return 0;return 1;}
I void GR_Init()//求原根
{
	RI i,j;for(i=2;i^m;++i) !((m-1)%i)&&IsP(i)&&(v[++cnt]=i);
	for(i=2;i^m;++i)//枚举
	{
		for(j=1;j<=cnt;++j) if(Qpow(i,(m-1)/v[j],m)%m==1) break;//验证
		if(j>cnt) {GR=i;break;}
	}
	for(i=0,j=1;i^(m-1);++i,j=1LL*j*GR%m) LG[j]=i;//预处理对数
}
class Poly
{
	private:
		int P,L,Inv,PR,IPR,R[4*M+5],A[4*M+5],B[4*M+5];
		I int Sum(CI x,CI y) {return x+y>=X?x+y-X:x+y;}I int Sub(CI x,CI y) {return x-y<0?x-y+X:x-y;}
		I void T(int *s,CI op)
		{
			RI i,j,k,U,S,x,y;for(i=0;i^P;++i) i<R[i]&&(s[i]^=s[R[i]]^=s[i]^=s[R[i]]);
			for(i=1;i^P;i<<=1) for(U=Qpow(~op?PR:IPR,(X-1)/(i<<1),X),j=0;j^P;j+=i<<1)
				for(S=1,k=0;k^i;++k,S=1LL*S*U%X) s[j+k]=Sum(x=s[j+k],y=1LL*S*s[i+j+k]%X),s[i+j+k]=Sub(x,y);
		}
	public:
		I Poly() {PR=3,IPR=Qpow(3,X-2,X);}
		I void Init()//因为数组大小固定，因此可以先预处理一波
		{
			P=1,L=0;W(P<=2*(m-2)) P<<=1,++L;Inv=Qpow(P,X-2,X);
			for(RI i=0;i^P;++i) R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
		}
		I void Mul(int *a,int *b,int *t)//NTT
		{
			RI i;for(i=0;i<=m-2;++i) A[i]=a[i],B[i]=b[i];for(;i^P;++i) A[i]=B[i]=0;
			for(T(A,1),T(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
			for(T(A,-1),i=0;i<=m-2;++i) t[i]=1LL*(A[i]+A[i+m-1])*Inv%X;//注意把后面的加到前面
		}
}NTT;
int main()
{
	RI i;for(scanf("%d%d%d%d",&n,&m,&tar,&k),i=1;i<=k;++i) scanf("%d",a+i);
	for(GR_Init(),ans[0]=1,i=1;i<=k;++i) a[i]&&(P[0][LG[a[i]]]=1);//初始化第一个数组
	for(NTT.Init(),i=1;i<=LN;++i) NTT.Mul(P[i-1],P[i-1],P[i]);//倍增
	for(i=0;i<=LN;++i) (n>>i)&1&&(NTT.Mul(ans,P[i],ans),0);//计算答案
	return printf("%d",ans[LG[tar]]),0;
}