[BZOJ3992][SDOI2015]序列统计(DP+原根+NTT)

时间：2022-11-26 08:00:17

3992: [SDOI2015]序列统计

Time Limit: 30 Sec  Memory Limit: 128 MB
Submit: 1888  Solved: 898
[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中的所有元素。
1<=N<=10^9，3<=M<=8000，M为质数
0<=x<=M-1，输入数据保证集合S中元素不重复x∈[1,m-1]

集合中的数∈[0,m-1]

Output

一行，一个整数，表示你求出的种类数mod 1004535809的值。

Sample Input

4 3 1 2
1 2

Sample Output

8
【样例说明】
可以生成的满足要求的不同的数列有(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)

HINT

Source

[Submit][Status][Discuss]

好久没有调过这么痛苦了。。

首先有一个简单的DP方程:$dp[i+j][(x*y)\%p]+=dp[i][x]*dp[j][y]$，dp[i][j]表示前i个数凑成余数为j的方案数。

$n^2$转移很简单，然后可以矩阵优化，但M的范围仍然过不了。

这时候就要敢往FFT方面去想。

现在这个方程之所以不能用FFT的原因在于x,y转移到的不是x+y而是x*y%mod，我们可以想到，乘法运算在作为指数的时候就是加法，又发现M是质数，于是我们考虑用原根的次幂替代S中的每个数。

设$g^{x'} \equiv x (mod\ p)$，$g^{y'} \equiv y (mod\ p)$，这样方程就变为$dp[i+j][(x'+y')\% \phi(p)]+=dp[i][x']*dp[j][y']$，这个就是标准的循环卷积了。

注意循环卷积次数界要再放大一倍！！还有S中要忽略0！！两个问题调到崩溃。

 #include<cstdio>

 #include<algorithm>

 #define rep(i,l,r) for (int i=l; i<=r; i++)

 typedef long long ll;

 using namespace std;

 const int N=,mod=,G=;

 int k,p,x,n,s,g,inv,cnt[N],f[N],pw[N],ind[N],rev[N],c[N];

 int ksm(int a,int b,int p){

     int res;

     for (res=; b; a=(1ll*a*a)%p,b>>=)

         if (b & ) res=(1ll*res*a)%p;

     return res;

 }

 bool chk(){

     for (int i=; i*i<=p; i++) if ((p-)%i== && ksm(g,(p-)/i,p)==) return ;

     return ;

 }

 void getroot(){

     if (p==) g=; else for (g=; !chk(); g++);

     ind[]=; pw[]=;

     for (int i=; i<p-; i++) pw[i]=pw[i-]*g%p,ind[pw[i]]=i;

 }

 namespace NTT{

     int n,L,rev[N];

     void init(int m){

         n=; L=;

         for (; n<=m; n<<=) L++;

         n<<=; L++; inv=ksm(n,mod-,mod);

         for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));

     }

     void DFT(int a[],int n,int f){

         for (int i=; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);

         for (int i=; i<n; i<<=){

             int wn=ksm(G,(f==) ? (mod-)/(i<<) : (mod-)-(mod-)/(i<<),mod);

             for (int p=i<<,j=; j<n; j+=p){

                 int w=;

                 for (int k=; k<i; k++,w=1ll*w*wn%mod){

                     int x=a[j+k],y=1ll*w*a[i+j+k]%mod;

                     a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;

                 }

             }

         }

         if (f==) return;

         for (int i=; i<n; i++) a[i]=1ll*a[i]*inv%mod;

     }

     void mul(int a[],int b[]){

         for (int i=; i<n; i++) c[i]=b[i];

         DFT(a,n,); DFT(c,n,);

         for (int i=; i<n; i++) a[i]=1ll*a[i]*c[i]%mod;

         DFT(a,n,-);

         for (int i=n-; i>=p-; i--) a[i-p+]=(a[i-p+]+a[i])%mod,a[i]=;

     }

 };

 int main(){

     freopen("bzoj3992.in","r",stdin);

     freopen("bzoj3992.out","w",stdout);

     scanf("%d%d%d%d",&k,&p,&x,&n); getroot(); NTT::init(p);

     rep(i,,n) { scanf("%d",&s); if (s) cnt[ind[s]]++; }

     for (f[]=; k; k>>=,NTT::mul(cnt,cnt)) if (k & ) NTT::mul(f,cnt);

     printf("%d\n",f[ind[x]]);

     return ;

 }

相关文章