BZOJ 3992: [SDOI2015]序列统计 NTT+快速幂

Posted 2020-09-05 Oncle_Ha

3992: [SDOI2015]序列统计

Time Limit: 30 Sec  Memory Limit: 128 MB
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

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

Round 1 感谢yts1999上传

 

想法：

  设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)

  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 }

View Code