[SDOI2015]序列统计

Posted 2021-02-22 troverld

XIII.[SDOI2015]序列统计

一个非常naive的想法就是多项式快速幂。

我们令一个函数(f_1(x)=[xin S])。并有(f_i(x)=sumlimits_{ijequiv xmod m}f_{i-1}(i)*f_{i-1}(j))。则答案为(f_n(x))。

后面那个奇妙的卷积不是任何我们已知的卷积，倒是有点长得像狄利克雷卷积。我们只能通过(m^2)暴力计算。再乘上多项式快速幂的(log n)，复杂度(O(m^2log n))，期望得分(60\%)。

有什么东西可以化乘为加呢？对数运算。

因为一个质数的原根的(0,1,dots,m-1)次方恰好覆盖了(0)到(m-1)内所有整数，所以你如果对一个(0)到(m-1)内的整数求关于原根的对数的话，求出来的结果是互不相同的。至于怎么求对数呢，直接将原根的各次幂打个表，求对数时直接查表即可。

我们现在求完对数得到了如下递推式：

(f_i(x)=sumlimits_{log_gi+log_gjequivlog_gxmod varphi(m)}f_{i-1}(i)*f_{i-1}(j))

是卷积的形式，直接NTT！

不过有一些注意点。因为求了对数，所以我们在NTT时实际上是在对指数瞎搞，因此模的应该是(varphi(m))，即(m-1)。

具体来说，这是卷积的代码：

void mul(int *a,int *b,int *c){
	for(int i=0;i<lim;i++)f[i]=a[i],g[i]=b[i];
	NTT(f,1),NTT(g,1);
	for(int i=0;i<lim;i++)f[i]=(1ll*f[i]*g[i])%mod;
	NTT(f,-1);
	for(int i=0;i<m-1;i++)c[i]=(f[i]+f[i+m-1])%mod;
}

原本是对(m)同余的加在一起，现在是对(m-1)同余的加在一起。

至于如何求出任意质数的原根吗……自己看其它JULAO的blog吧

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod=1004535809,G=3;
int n,m,X,S,lim=1,invlim,lg,rev[1<<16],f[1<<16],g[1<<16],bs[1<<16],qwq[1<<16];
int pov(int x,int y,int z){
	int res=1;
	for(;y;x=(1ll*x*x)%z,y>>=1)if(y&1)res=(1ll*res*x)%z;
	return res;
}
int getroot(int x){
	vector<int>v;v.clear();
	int phi=x-1;
	for(int i=2;i*i<=phi;i++){
		if(phi%i)continue;
		v.push_back(i);
		while(!(phi%i))phi/=i;
	}
	if(phi>1)v.push_back(phi);
	phi=x-1;
	for(int i=2;i<=phi;i++){
		bool ok=true;
		for(int j=0;j<v.size();j++)if(pov(i,phi/v[j],x)==1){ok=false;break;}
		if(ok)return i;
	}
	return 19260817;
}
void NTT(int *a,int tp){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int md=1;md<lim;md<<=1){
		int rt=pov(G,(mod-1)/(md<<1),mod);
		if(tp==-1)rt=pov(rt,mod-2,mod);
		for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
			int w=1;
			for(int i=0;i<md;i++,w=(1ll*w*rt)%mod){
				int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
				a[pos+i]=(x+y)%mod;
				a[pos+md+i]=(x-y+mod)%mod;
			}
		}
	}
	if(tp==-1)for(int i=0;i<lim;i++)a[i]=(1ll*a[i]*invlim)%mod;
}
void mul(int *a,int *b,int *c){
	for(int i=0;i<lim;i++)f[i]=a[i],g[i]=b[i];
	NTT(f,1),NTT(g,1);
	for(int i=0;i<lim;i++)f[i]=(1ll*f[i]*g[i])%mod;
	NTT(f,-1);
	for(int i=0;i<m-1;i++)c[i]=(f[i]+f[i+m-1])%mod;
}
void ksm(int x){
	if(x==1){for(int i=0;i<lim;i++)qwq[i]=bs[i];return;}
	ksm(x>>1);
	mul(qwq,qwq,qwq);
	if(x&1)mul(qwq,bs,qwq);
}
map<int,int>mp;
int main(){
	scanf("%d%d%d%d",&n,&m,&X,&S);
	int GG=getroot(m);
	for(int i=0;i<m-1;i++)mp[pov(GG,i,m)]=i;
	for(int i=0,x;i<S;i++){
		scanf("%d",&x);
		if(x)bs[mp[x]]=1;
	}
	while(lim<=m*2)lim<<=1,lg++;
	invlim=pov(lim,mod-2,mod);
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	ksm(n);
	printf("%d
",qwq[mp[X]]);
	return 0;
}