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

Round 1 感謝yts1999上傳

[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！！兩個問題調到崩潰。

 1 #include<cstdio>
 2 #include<algorithm>
 3 #define rep(i,l,r) for (int i=l; i<=r; i++)
 4 typedef long long ll;
 5 using namespace std;
 6 
 7 const int N=20100,mod=1004535809,G=3;
 8 int k,p,x,n,s,g,inv,cnt[N],f[N],pw[N],ind[N],rev[N],c[N];
 9 
10 int ksm(int a,int b,int p){
11     int res;
12     for (res=1; b; a=(1ll*a*a)%p,b>>=1)
13         if (b & 1) res=(1ll*res*a)%p;
14     return res;
15 }
16 
17 bool chk(){
18     for (int i=2; i*i<=p; i++) if ((p-1)%i==0 && ksm(g,(p-1)/i,p)==1) return 0;
19     return 1;
20 }
21 
22 void getroot(){
23     if (p==2) g=1; else for (g=2; !chk(); g++);
24     ind[1]=0; pw[0]=1;
25     for (int i=1; i<p-1; i++) pw[i]=pw[i-1]*g%p,ind[pw[i]]=i;
26 }
27 
28 namespace NTT{
29     int n,L,rev[N];
30     void init(int m){
31         n=1; L=0;
32         for (; n<=m; n<<=1) L++;
33         n<<=1; L++; inv=ksm(n,mod-2,mod);
34         for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
35     }
36     void DFT(int a[],int n,int f){
37         for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
38         for (int i=1; i<n; i<<=1){
39             int wn=ksm(G,(f==1) ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1),mod);
40             for (int p=i<<1,j=0; j<n; j+=p){
41                 int w=1;
42                 for (int k=0; k<i; k++,w=1ll*w*wn%mod){
43                     int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
44                     a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
45                 }
46             }
47         }
48         if (f==1) return;
49         for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
50     }
51     void mul(int a[],int b[]){
52         for (int i=0; i<n; i++) c[i]=b[i];
53         DFT(a,n,1); DFT(c,n,1);
54         for (int i=0; i<n; i++) a[i]=1ll*a[i]*c[i]%mod;
55         DFT(a,n,-1);
56         for (int i=n-1; i>=p-1; i--) a[i-p+1]=(a[i-p+1]+a[i])%mod,a[i]=0;
57     }
58 };
59 
60 int main(){
61     freopen("bzoj3992.in","r",stdin);
62     freopen("bzoj3992.out","w",stdout);
63     scanf("%d%d%d%d",&k,&p,&x,&n); getroot(); NTT::init(p);
64     rep(i,1,n) { scanf("%d",&s); if (s) cnt[ind[s]]++; }
65     for (f[0]=1; k; k>>=1,NTT::mul(cnt,cnt)) if (k & 1) NTT::mul(f,cnt);
66     printf("%d\n",f[ind[x]]);
67     return 0;
68 }

[BZOJ3992][SDOI2015]序列統計(DP+原根+NTT)