「SDOI2015」序列统计

Source and Judge

SDOI2015
luogu3321

Record

1h

Analysis

对M求原根g，对每个s以及x求g下的离散对数，化乘法为加法
然后相当于对生成函数求循环卷积（循环长度为phi(M)，需要手动），ntt即可
注意si可能为0

//Zory-2018
#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<map>
#include<set>
#include<queue>
#include<deque>
#include<stack>
#include<bitset>
#include<vector>
#include<algorithm>
#include<iostream>
using namespace std;

namespace mine
{
	typedef long long ll;
	const int INF=0x3f3f3f3f;
	//#define pr pair<double,int>
	#define FR first
	#define SE second
	#define MP make_pair

	const int MAX_NUM=1e4;
	bool no[MAX_NUM];
	int pr=0,prime[MAX_NUM],phi[MAX_NUM];
	void pre()
	{
		for(int i=2;i<MAX_NUM;i++)
		{
			if(!no[i]) prime[++pr]=i,phi[i]=i-1;
			for(int j=1;j<=pr and (ll)i*prime[j]<MAX_NUM;j++)
			{
				int t=i*prime[j];
				no[t]=1;
				if(i%prime[j]==0) {phi[t]=phi[i]*prime[j];break;}
				phi[t]=phi[i]*(prime[j]-1);
			}
		}
	}

	int MOD=1004535809,M;
	ll qpower(ll x,ll e,ll MOD)
	{
		ll ans=1;
		while(e)
		{
			if(e&1) ans=ans*x%MOD;
			x=x*x%MOD;e>>=1;
		}
		return ans;
	}
	ll inv(ll x,ll MOD) {return qpower(x,MOD-2,MOD);}
	ll add(ll x,ll y)
	{
		x+=y;
		if(x>=MOD) x-=MOD;
		if(x<=-MOD) x+=MOD;
		return x;
	}

	bool check(int g,int num)
	{
		int tmp=phi[num];
		for(int i=1;i<=pr and (ll)prime[i]*prime[i]<=tmp;i++) if(tmp%prime[i]==0)
		{
			if(qpower(g,phi[num]/prime[i],num)==1) return 0;
			while(tmp%prime[i]==0) tmp/=prime[i];
		}
		return qpower(g,tmp,num)!=1;
	}
	int getyg(int num)
	{
		for(int g=2;g<num;g++) if(check(g,num)) return g;
		return -1;
	}

	const int MAX_N=MAX_NUM*4;
	struct Ntt
	{
		int R[MAX_N];
		void preR(int n)
		{
			R[0]=0;for(int i=1;i<n;i++) R[i]=(R[i>>1]>>1)|((i&1)<<((int)log2(n)-1));
		}
		ll w[2][MAX_N];
		void prew(int n)
		{
			ll g=qpower(3,(MOD-1)/n,MOD),g2=inv(g,MOD);
			w[0][0]=w[1][0]=1;
			for(int i=1;i<n;i++) w[0][i]=w[0][i-1]*g%MOD,w[1][i]=w[1][i-1]*g2%MOD;
		}
		void solve(ll A[],int n,int f)
		{
			for(int i=0;i<n;i++) if(i<R[i]) swap(A[i],A[R[i]]);
			for(int ln=1;ln<=n/2;ln<<=1)
				for(int st=0;st<n;st+=ln*2)
					for(int k=0;k<ln;k++)
					{
						ll a=A[st+k],b=w[f][k*(n/(ln*2))]*A[st+ln+k]%MOD;
						A[st+k]=add(a,b);A[st+ln+k]=add(a,-b);
					}
		}
	}ntt;
	ll ans[MAX_N],A[MAX_N],C[MAX_N];
	void cheng(int ln,int n,ll X[],ll Y[])
	{
		memcpy(C,Y,sizeof C);
		ntt.solve(X,ln,0);ntt.solve(C,ln,0);
		for(int i=0;i<ln;i++) X[i]=X[i]*C[i]%MOD;
		ntt.solve(X,ln,1);
		for(int i=0;i<n;i++) X[i]=(X[i]+X[n+i])*inv(ln,MOD)%MOD,X[n+i]=0;
	}
	void power(int n,int e)
	{
		int ln=1;while(ln<4*n) ln<<=1;

		ntt.preR(ln);ntt.prew(ln);
		ans[0]=1;
		while(e)
		{
			if(e&1) cheng(ln,n,ans,A);
			e>>=1;cheng(ln,n,A,A);
		}
	}

	bool app[MAX_N];
	void main()
	{
		pre();
		int k,wt,n;scanf("%d%d%d%d",&k,&M,&wt,&n);
		for(int i=1;i<=n;i++) {int t;scanf("%d",&t);app[t]=1;}
		int g=getyg(M);
		bool bk=1;//debug
		for(int now=1,t=0;t<=M-2;t++,now=now*g%M)
		{
			if(app[now]) A[t]=1;
			if(bk==1 and now==wt) wt=t,bk=0;
		}
		power(M-1,k);
		printf("%lld",(ans[wt]+MOD)%MOD);
	}
};
int main()
{
	srand(time(0));
	mine::main();
}