loj2183/洛谷P3321/bzoj3992 序列统计 原根+NTT

题目分析

哎呀原根这个东西忘得差不多了…=。=

对于$P$的剩余系的原根$g$，用$g^k\mathrm{m}\mathrm{o}\mathrm{d}P$可以表示$P$的剩余系中的所有数。

原根的定义是，对于一个群$G$，单位元为$e$。将满足$a^d=e$的最小正整数$d$记作$a$的阶，也就是$ord(a)$，那么使得$ord(g)=|G|$的$g$就是$G$的原根。

由于原根通常情况下都不大，所以我们可以暴力枚举这个原根。那么$g$是$P$的剩余系的原根的条件是，$g^{P-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}P)$且对于任意$P-1$的质因子$p_i$，都有$g^{\frac{P-1}{p_i}}̸\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}P)$。后面那一个条件是为了满足$P-1$是满足$g^d\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}P)$的最小正整数，因为$1$的幂次都是$1$，所以该条件满足时，$P-1$无论去掉多少个质因子，$g$的这么多次方模$P$都不会为1。

找到$m$的原根后，若$g^k(\mathrm{m}\mathrm{o}\mathrm{d}m)$在$|S|$中出现过，则$a_k=1$，否则$a_k=0$。那么将$a$做$n$次最高次项为$m-1$的循环卷积，就能得到选出$n$个数，乘积模$m$为$g^k(\mathrm{m}\mathrm{o}\mathrm{d}m)$的方案数。用快速幂实现即可。

循环卷积？就是做一次NTT后，将$a_i+=a_{m-1+i}$，然后$a_{m-1+i}=0$即可。

代码

#include<bits/stdc++.h>
using namespace std;
#define RI register int
int read() {
	int q=0;char ch=' ';
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9') q=q*10+ch-'0',ch=getchar();
	return q;
}
const int G=3,mod=1004535809,N=16390;
int n,m,x,ss,gm,kn,len;
int A[N],k1[N],k2[N],res[N],rev[N],gk[N];

int qm(int x,int p) {return x>=p?x-p:x;}
int ksm(int x,int y,int p) {
	int re=1;
	for(;y;y>>=1,x=1LL*x*x%p) if(y&1) re=1LL*re*x%p;
	return re;
}

void NTT(int *a,int n,int x) {
	for(RI i=0;i<n;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
	for(RI i=1;i<n;i<<=1) {
		int gn=ksm(G,(mod-1)/(i<<1),mod);
		for(RI j=0;j<n;j+=i<<1) {
			int g=1,t1,t2;
			for(RI k=0;k<i;++k,g=1LL*g*gn%mod) {
				t1=a[j+k],t2=1LL*g*a[j+i+k]%mod;
				a[j+k]=qm(t1+t2,mod),a[j+i+k]=qm(t1-t2+mod,mod);
			}
		}
	}
	if(x==1) return;
	int inv=ksm(n,mod-2,mod);reverse(a+1,a+n);
	for(RI i=0;i<n;++i) a[i]=1LL*a[i]*inv%mod;
}
void mul(int *a,int *b) {
	for(RI i=0;i<kn;++i) k1[i]=a[i],k2[i]=b[i];
	NTT(k1,kn,1),NTT(k2,kn,1);
	for(RI i=0;i<kn;++i) a[i]=1LL*k1[i]*k2[i]%mod;
	NTT(a,kn,-1);
	for(RI i=m-1;i<kn;++i)
		a[i%(m-1)]=qm(a[i%(m-1)]+a[i],mod),a[i]=0;
}

void work1() {
	int flag=0;
	for(RI i=1;i<=ss;++i) if(read()==0) flag=1;
	if(flag) printf("%d\n",qm(ksm(ss,n,mod)-ksm(ss-1,n,mod)+mod,mod));
	else puts("0");
}
int pri[8005],vis[8005],js;
int getgm(int x) {
	int kx=x-1;
	for(RI i=2;i*i<=kx;++i)
		if(kx%i==0) {pri[++js]=i;while(kx%i==0) kx/=i;}
	if(kx!=1) pri[++js]=kx;
	for(RI i=2;;++i) {
		int flag=1;
		if(ksm(i,x-1,x)!=1) continue;
		for(RI j=1;j<=js;++j)
			if(ksm(i,(x-1)/pri[j],x)==1) {flag=0;break;}
		if(flag) return i;
	}
}
void work2() {
	res[0]=1;
	kn=1;while(kn<=2*m-4) kn<<=1,++len;
	for(RI i=0;i<kn;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	for(RI i=n;i;i>>=1,mul(A,A)) if(i&1) mul(res,A);
	printf("%d\n",res[gk[x]]);
}
int main()
{
	n=read(),m=read(),x=read(),ss=read();
	if(x==0) {work1();return 0;}
	gm=getgm(m);
	for(RI i=1;i<=ss;++i) vis[read()]=1;
	for(RI now=1,i=0;i<m-1;++i,now=1LL*now*gm%m) gk[now]=i,A[i]=vis[now];
	work2();
    return 0;
}