1. 程式人生 > >BZOJ3992: [SDOI2015]序列統計

BZOJ3992: [SDOI2015]序列統計

++ src ons close pla std lose == con

https://www.lydsy.com/JudgeOnline/problem.php?id=3992

有a*bΞx(mod m),則ord(a)+ord(b)Ξord(x)(mod φ(m))

對Σxord(si)循環卷積FFT即可

技術分享圖片
#include<cstdio>
#include<algorithm>
#define re register 
#define rep(i,s,t) for(re int i=s;i<=t;++i)
#define Rep(i,s,t) for(re int i=s;i<t;++i)
using
namespace std; const int N=4e5+11,mod=1004535809; int n,m,x,s,g,sum,d,lg2; int q[N],ord[N],a[N],p[N],res[N]; inline void inc(int &a,int b){ a+=b; if(a>=mod)a-=mod; if(a<0)a+=mod; } inline int fp(int a,int b,int m=mod){ if(b<0)b+=m-1; if(b>=m-1)b-=m-1;
int res=1; for(;b;b>>=1,a=1ll*a*a%m) if(b&1) res=1ll*res*a%m; return res; } inline int findg(int x){ int t=0; for(re int i=2;i*i<=x-1;++i) if((x-1)%i==0){ q[++t]=i; if(i*i!=x-1) q[++t]=(x-1)/i; } rep(i,
2,(1<<30)){ rep(j,1,t) if(fp(i,q[j],x)==1) goto end; return i; end:; } } inline void dft(int *a,int n,int f){ Rep(i,0,n) if(i<p[i]) swap(a[i],a[p[i]]); for(re int i=1,x,y;i<n;i<<=1){ for(re int w,j=0,wn=fp(3,(mod-1)/(i<<1)*f);w=1,j<n;j+=(i<<1)){ for(re int k=j;k<i+j;++k,w=1ll*w*wn%mod){ x=a[k],y=1ll*w*a[i+k]%mod; a[k]=(x+y)%mod; a[i+k]=(x-y+mod)%mod; } } } if(f==-1){ int inv=fp(n,mod-2); Rep(i,0,d) a[i]=1ll*a[i]*inv%mod; } } int main(){ //freopen("test.in","r",stdin); scanf("%d%d%d%d",&n,&m,&x,&s); g=findg(m),sum=1; //printf("%d\n",g); rep(j,0,m-1){ ord[sum]=j; sum=1ll*sum*g%m; } for(int x;s--;){ scanf("%d",&x); if(x)++a[ord[x]]; } for(d=1,lg2=-1;d<=m+m;d<<=1,++lg2); Rep(i,0,d)p[i]=(p[i>>1]>>1)^((i&1)<<lg2); Rep(i,0,d)res[i]=a[i]; for(--n;n;n>>=1){ dft(a,d,1); if(n&1){ dft(res,d,1); Rep(i,0,d)res[i]=1ll*res[i]*a[i]%mod; dft(res,d,-1); Rep(i,0,m-1)inc(res[i],res[i+m-1]); Rep(i,m-1,d)res[i]=0; } Rep(i,0,d)a[i]=1ll*a[i]*a[i]%mod; dft(a,d,-1); Rep(i,0,m-1)inc(a[i],a[i+m-1]); Rep(i,m-1,d)a[i]=0; } printf("%d\n",res[ord[x]]); return 0; }
code

BZOJ3992: [SDOI2015]序列統計