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

来源：互联网发布：手机截取音乐软件编辑：程序博客网时间：2024/04/24 17:32

首先设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);}

2 0