BZOJ 3992: [SDOI2015]序列统计 NTT+快速幂
3992: [SDOI2015]序列统计
Submit: 1155 Solved: 532
[Submit][Status][Discuss]
Description
小C有一个集合S，里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数列，数列中的每个数都属于集合S。
小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数x，求所有可以生成出的，且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为，两个数列{Ai}和{Bi}不同，当且仅当至少存在一个整数i，满足Ai≠Bi。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案mod 1004535809的值就可以了。
Input
一行，四个整数，N、M、x、|S|，其中|S|为集合S中元素个数。第二行，|S|个整数，表示集合S中的所有元素。
Output
一行，一个整数，表示你求出的种类数mod 1004535809的值。
Sample Input
4 3 1 2
1 2
1 2
Sample Output
8
HINT
【样例说明】
可以生成的满足要求的不同的数列有(1,1,1,1)、(1,1,2,2)、(1,2,1,2)、(1,2,2,1)、(2,1,1,2)、(2,1,2,1)、(2,2,1,1)、(2,2,2,2)。
【数据规模和约定】
对于10%的数据，1<=N<=1000；
对于30%的数据，3<=M<=100；
对于60%的数据，3<=M<=800；
对于全部的数据，1<=N<=109，3<=M<=8000，M为质数，1<=x<=M-1，输入数据保证集合S中元素不重复
Source
想法：
设a[i]表示数字i是否属于集合S，C[i]表示数列之积%M=i的方案数。 当n=2时：c[(i*j)%M]=∑a[i]*a[j]
令A[i]=a[g^i],C[i]=c[g^i]//∵g为M原根，遍历0~M-1，而将数组映射到另一个数组，并不影响答案，只要改变运算规则。
由 c[g^(i+j)%M]=∑a[g^i]*a[g^j] 得到：
C[(i+j)%(M-1)]=∑A[i]*A[j]//费马小定理：g^(M-1)
∵j+i≤m*2
∴每次FFT后将后面的累加到前面来就行了。
当n=y时,C=A^y，找到g^j=x,输出C[j] 于是NTT+快速幂O(nlog^2n)
∴每次FFT后将后面的累加到前面来就行了。
当n=y时,C=A^y，找到g^j=x,输出C[j] 于是NTT+快速幂O(nlog^2n)
1 #include<cstdio> 2 #define ll long long 3 const int MP(1004535809),lem(16400),g(3); 4 int n,m,x,size,gm; 5 int a[lem+10],num,p[30],tp; 6 struct data{int a[lem+10];}A,C,B; 7 int power(int a,int b,int MP) 8 { 9 ll t=1,y=a;b+=b<0?MP-1:0; 10 for(;b;b>>=1,y=(y*y)%MP)if(b&1)t=(t*y)%MP; 11 return (int)t; 12 } 13 bool check(int y) 14 { 15 for(int j=1;j<=tp;j++) 16 if(power(y,(m-1)/p[j],m)==1)return false; 17 return true; 18 } 19 void Get_g(int x) 20 { 21 if(!(x&1)) 22 { 23 p[++tp]=2; 24 while(!(x&1))x>>=1; 25 } 26 for(int i=3;i*i<=x;i+=2) 27 { 28 if(x%i==0) 29 { 30 p[++tp]=i; 31 while(x%i==0)x/=i; 32 } 33 } 34 if(x>1)p[++tp]=x; 35 for(int i=2;i<=m-1;i++) 36 if(check(i)){gm=i;break;} 37 } 38 int R[lem+10],w[lem+10],wn,l,il,h; 39 void deal() 40 { 41 l=1;w[0]=1; 42 while(l<=m+m)l<<=1,h++; 43 for(int i=0;i<l;i++)R[i]=(R[i>>1]>>1|(i&1)<<(h-1)); 44 il=power(l,MP-2,MP); 45 } 46 void swap(int &a,int &b){if(a==b)return;a^=b;b^=a;a^=b;} 47 void NTT(int *a,int l,int ty) 48 { 49 for(int i=0;i<l;i++)if(i<R[i])swap(a[i],a[R[i]]); 50 for(int leng=2;leng<=l;leng<<=1) 51 { 52 int M=leng>>1; 53 wn=power(g,ty*(MP-1)/leng,MP); 54 for(int i=1;i<M;i++)w[i]=(1ll*w[i-1]*wn)%MP; 55 for(int i=0;i<l;i+=leng) 56 { 57 for(int j=0;j<M;j++) 58 { 59 int x=a[i+j],y=(1ll*w[j]*a[i+j+M])%MP; 60 a[i+j]=x+y;a[i+j+M]=x-y; 61 a[i+j]-=a[i+j]>=MP?MP:0; 62 a[i+j+M]+=a[i+j+M]<0?MP:0; 63 } 64 } 65 } 66 if(ty==-1) 67 for(int i=0;i<l;i++)a[i]=(1ll*a[i]*il)%MP; 68 } 69 void three(data &A,data &B) 70 { 71 NTT(A.a,l,1);NTT(B.a,l,1); 72 for(int i=0;i<l;i++)B.a[i]=(1ll*B.a[i]*A.a[i])%MP; 73 NTT(B.a,l,-1); 74 for(int i=m-1;i<l;i++)B.a[i%(m-1)]=(B.a[i%(m-1)]+B.a[i])%MP,B.a[i]=0; 75 } 76 void run() 77 { 78 C.a[0]=1; 79 while(n) 80 { 81 if(n&1) 82 { 83 B=A; 84 three(B,C); 85 } 86 B=A; 87 three(B,A); 88 n>>=1; 89 } 90 } 91 int main() 92 { 93 scanf("%d%d%d%d",&n,&m,&x,&size); 94 for(int i=1;i<=size;i++){scanf("%d",&num);a[num]=1;} 95 Get_g(m-1);deal(); 96 for(int i=0,j=1;i<m-1;i++,j=(j*gm)%m) 97 { 98 A.a[i]=a[j]; 99 if(j==x)num=i; 100 } 101 run(); 102 printf("%d",C.a[num]); 103 return 0; 104 }
