Bzoj3992:[SDOI2015]序列统计:NTT+DP

首先设dp[i][j]表示前i个数余j的方案数

dp[i][j]->dp[i+1][j*s[k]%m]=>O(n*m^2);

发现M是个质数，有原根

将j用j的原根表示为j'

则dp[i][j']=∑dp[i-1][k']((k'+s[p]')%(m-1)==j')

这是一个卷积，可以用NTT优化=>O(nmlogm)

对n快速幂一下就是O(mlognlogm)

#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=30010;
const int G=3;
const int mod=1004535809;
int n,m,x,sum,s[maxn],a[maxn],b[maxn],g,inv,N;
int A[maxn],B[maxn],C[maxn],tim=0,flag[maxn],D[maxn];

int quick_pow_m(int x,int y){
	int ret=1;
	while (y){
		if (y&1) ret=1ll*ret*x%m;
		x=1ll*x*x%m; y>>=1;
	}return ret;
}

int quick_pow_mod(int x,int y){
	int ret=1;
	while (y){
		if (y&1) ret=1ll*ret*x%mod;
		x=1ll*x*x%mod; y>>=1;
	}return ret;
}

bool getg(int x){
	++tim;
	for (int i=1;i<m;++i) flag[quick_pow_m(x,i)]=tim;
	for (int i=1;i<m;++i) if (flag[i]!=tim) return 0;
	return 1;
}

void NTT(int a[],int flag){
	for (int i=0,j=0;i<N;++i){
		if (i>j) swap(a[i],a[j]);
		for (int t=N>>1;(j^=t)<t;t>>=1);
	}
	for (int i=2;i<=N;i<<=1){
		int wn=quick_pow_mod(G,flag==1?(mod-1)/i:(mod-1-(mod-1)/i));
		for (int k=0;k<N;k+=i){
			int w=1;
			for (int j=k;j<k+i/2;++j){
				int x=a[j],y=1ll*a[j+i/2]*w%mod;
				a[j]=(1ll*x+1ll*y)%mod; a[j+i/2]=(x-y+mod)%mod;
				w=1ll*w*wn%mod;
			}
		}
	}
	if (flag==-1)for (int i=0;i<N;++i) a[i]=1LL*a[i]*inv%mod;
}

void solve1(){
	for (int i=0;i<m;++i) C[i]=A[i],D[i]=B[i];
	for (int i=m;i<N;++i) C[i]=D[i]=0;
	
	NTT(D,1); NTT(C,1);
	for (int i=0;i<N;++i) B[i]=1ll*C[i]*D[i]%mod;
	NTT(B,-1);
	for (int i=m-1;i<N;++i) B[i%(m-1)]+=B[i],B[i%(m-1)]%=mod,B[i]=0;
}

void solve2(){
	NTT(A,1);
	for (int i=0;i<N;++i) A[i]=1ll*A[i]*A[i]%mod;
	NTT(A,-1);
	for (int i=m-1;i<N;++i) A[i%(m-1)]+=A[i],A[i%(m-1)]%=mod,A[i]=0;
}

void polynomial_quick_pow(int y){
	for (int i=1;i<=sum;++i) if (!s[i]) continue; else A[b[i]]=1;
	B[0]=1;
	while (y){
		if (y&1) solve1();
		solve2(); y>>=1;
	}
	printf("%d\n",B[a[x]]);
}

int main(){
	scanf("%d%d%d%d",&n,&m,&x,&sum);
	for (N=1;N<=(2*m-2);N<<=1);
	inv=quick_pow_mod(N,mod-2);
	for (int i=1;i<=sum;++i) scanf("%d",&s[i]),s[i]%=m;
	for (int i=2;i;++i) if(getg(i)){g=i;break;}
	for (int i=0;i<m-1;++i) a[quick_pow_m(g,i)]=i;
	for (int i=1;i<=sum;++i) b[i]=a[s[i]];
    polynomial_quick_pow(n);
}