[BZOJ 3992][SDOI2015]序列统计

  1 #include <bits/stdc++.h>
  2 
  3 const int G=3;
  4 const int DFT=1;
  5 const int IDFT=-1;
  6 const int MAXN=550000;
  7 const int MOD=1004535809;
  8 const int INV2=(MOD+1)>>1;
  9 const int PHI=MOD-1;
 10 
 11 typedef std::vector<int> Poly;
 12 
 13 Poly Sqrt(Poly);
 14 void Read(Poly&);
 15 Poly Inverse(Poly);
 16 Poly Ln(const Poly&);
 17 Poly Exp(const Poly&);
 18 void Print(const Poly&);
 19 void NTT(Poly&,int,int);
 20 Poly Pow(const Poly&,int);
 21 Poly Integral(const Poly&);
 22 Poly Derivative(const Poly&);
 23 Poly operator*(Poly,Poly);
 24 Poly operator/(Poly,Poly);
 25 Poly operator%(Poly,Poly);
 26 Poly operator<<(const Poly&,int);
 27 Poly operator>>(const Poly&,int);
 28 Poly operator+(const Poly&,const Poly&);
 29 Poly operator-(const Poly&,const Poly&);
 30 
 31 int rev[MAXN];
 32 
 33 int FindG(int);
 34 int NTTPre(int);
 35 int Sqrt(int,int);
 36 int Pow(int,int,int);
 37 int Log(int,int,int);
 38 int ExGCD(int,int,int&,int&);
 39 
 40 int k,m,x,n;
 41 int lg[MAXN];
 42 
 43 int main(){
 44     scanf("%d%d%d%d",&k,&m,&x,&n);
 45     Poly a(m-1);
 46     int g=FindG(m);
 47     for(int i=0,pw=1;i<m-1;i++,pw=pw*g%m)
 48         lg[pw]=i;
 49     for(int i=0;i<n;i++){
 50         int x;
 51         scanf("%d",&x);
 52         if(x)
 53             a[lg[x]]=1;
 54     }
 55     --m;
 56     Poly ans(a);
 57     --k;
 58     while(k>0){
 59         if(k&1){
 60             ans=a*ans;
 61             ans=ans+(ans>>m);
 62             ans.resize(m);
 63         }
 64         a=a*a;
 65         a=a+(a>>m);
 66         a.resize(m);
 67         k>>=1;
 68     }
 69     printf("%d\n",ans[lg[x]]);
 70     return 0;
 71 }
 72 
 73 Poly operator<<(const Poly& a,int x){
 74     Poly ans(a.size()+x);
 75     for(size_t i=0;i<a.size();i++)
 76         ans[i+x]=a[i];
 77     return ans;
 78 }
 79 
 80 Poly operator>>(const Poly& a,int x){
 81     Poly ans(a.size()-x);
 82     for(size_t i=x;i<a.size();i++)
 83         ans[i-x]=a[i];
 84     return ans;
 85 }
 86 
 87 int FindG(int mod){
 88     for(int g=2;g<mod;g++){
 89         bool flag=true;
 90         for(int i=1,pw=g;i<mod-1;i++,pw=pw*g%mod){
 91             if(pw==1){
 92                 flag=false;
 93                 break;
 94             }
 95         }
 96         if(flag)
 97             return g;
 98     }
 99     return -1;
100 }
101 
102 void Read(Poly& a){
103     for(auto& i:a)
104         scanf("%d",&i);
105 }
106 
107 void Print(const Poly& a){
108     for(auto i:a)
109         printf("%d ",i);
110     puts("");
111 }
112 
113 Poly Pow(const Poly& a,int k){
114     Poly log=Ln(a);
115     for(auto& i:log)
116         i=1ll*i*k%MOD;
117     return Exp(log);
118 }
119 
120 Poly Sqrt(Poly a){
121     int len=a.size();
122     if(len==1)
123         return Poly(1,Sqrt(a[0],MOD));
124     else{
125         Poly b=a;
126         b.resize((len+1)>>1);
127         b=Sqrt(b);
128         b.resize(len);
129         Poly inv=Inverse(b);
130         int bln=NTTPre(inv.size()+a.size());
131         NTT(a,bln,DFT);
132         NTT(inv,bln,DFT);
133         for(int i=0;i<bln;i++)
134             a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
135         NTT(a,bln,IDFT);
136         for(int i=0;i<len;i++)
137             b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
138         return b;
139     }
140 }
141 
142 Poly Exp(const Poly& a){
143     size_t len=1;
144     Poly ans(1,1),one(1,1);
145     while(len<(a.size()<<1)){
146         len<<=1;
147         Poly b=a;
148         b.resize(len);
149         ans=ans*(one-Ln(ans)+b);
150         ans.resize(len);
151     }
152     ans.resize(a.size());
153     return ans;
154 }
155 
156 Poly Ln(const Poly& a){
157     Poly ans=Integral(Derivative(a)*Inverse(a));
158     ans.resize(a.size());
159     return ans;
160 }
161 
162 Poly Integral(const Poly& a){
163     int len=a.size();
164     Poly ans(len+1);
165     for(int i=1;i<len;i++)
166         ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
167     return ans;
168 }
169 
170 Poly Derivative(const Poly& a){
171     int len=a.size();
172     Poly ans(len-1);
173     for(int i=1;i<len;i++)
174         ans[i-1]=1ll*a[i]*i%MOD;
175     return ans;
176 }
177 
178 Poly operator/(Poly a,Poly b){
179     int n=a.size()-1,m=b.size()-1;
180     Poly ans(1);
181     if(n>=m){
182         std::reverse(a.begin(),a.end());
183         std::reverse(b.begin(),b.end());
184         b.resize(n-m+1);
185         ans=Inverse(b)*a;
186         ans.resize(n-m+1);
187         std::reverse(ans.begin(),ans.end());
188     }
189     return ans;
190 }
191 
192 Poly operator%(Poly a,Poly b){
193     int n=a.size()-1,m=b.size()-1;
194     Poly ans;
195     if(n<m)
196         ans=a;
197     else
198         ans=a-(a/b)*b;
199     ans.resize(m);
200     return ans;
201 }
202 
203 Poly operator*(Poly a,Poly b){
204     int len=a.size()+b.size()-1;
205     int bln=NTTPre(len);
206     NTT(a,bln,DFT);
207     NTT(b,bln,DFT);
208     for(int i=0;i<bln;i++)
209         a[i]=1ll*a[i]*b[i]%MOD;
210     NTT(a,bln,IDFT);
211     a.resize(len);
212     return a;
213 }
214 
215 Poly operator+(const Poly& a,const Poly& b){
216     Poly ans(std::max(a.size(),b.size()));
217     std::copy(a.begin(),a.end(),ans.begin());
218     for(size_t i=0;i<b.size();i++)
219         ans[i]=(ans[i]+b[i])%MOD;
220     return ans;
221 }
222 
223 Poly operator-(const Poly& a,const Poly& b){
224     Poly ans(std::max(a.size(),b.size()));
225     std::copy(a.begin(),a.end(),ans.begin());
226     for(size_t i=0;i<b.size();i++)
227         ans[i]=(ans[i]+MOD-b[i])%MOD;
228     return ans;
229 }
230 
231 Poly Inverse(Poly a){
232     int len=a.size();
233     if(len==1)
234         return Poly(1,Pow(a[0],MOD-2,MOD));
235     else{
236         Poly b(a);
237         b.resize((len+1)>>1);
238         b=Inverse(b);
239         int bln=NTTPre(b.size()*2+a.size());
240         NTT(a,bln,DFT);
241         NTT(b,bln,DFT);
242         for(int i=0;i<bln;i++)
243             b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
244         NTT(b,bln,IDFT);
245         b.resize(len);
246         return b;
247     }
248 }
249 
250 void NTT(Poly& a,int len,int opt){
251     a.resize(len);
252     for(int i=0;i<len;i++)
253         if(rev[i]>i)
254             std::swap(a[i],a[rev[i]]);
255     for(int i=1;i<len;i<<=1){
256         int step=i<<1;
257         int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
258         for(int j=0;j<len;j+=step){
259             int w=1;
260             for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
261                 int x=a[j+k];
262                 int y=1ll*a[j+k+i]*w%MOD;
263                 a[j+k]=(x+y)%MOD;
264                 a[j+k+i]=(x-y+MOD)%MOD;
265             }
266         }
267     }
268     if(opt==IDFT){
269         int inv=Pow(len,MOD-2,MOD);
270         for(int i=0;i<len;i++)
271             a[i]=1ll*a[i]*inv%MOD;
272     }
273 }
274 
275 int NTTPre(int n){
276     int bln=1,bct=0;
277     while(bln<n){
278         bln<<=1;
279         ++bct;
280     }
281     for(int i=0;i<bln;i++)
282         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
283     return bln;
284 }
285 
286 inline int Pow(int a,int n,int p){
287     int ans=1;
288     while(n>0){
289         if(n&1)
290             ans=1ll*a*ans%p;
291         a=1ll*a*a%p;
292         n>>=1;
293     }
294     return ans;
295 }
296 
297 int ExGCD(int a,int b,int& x,int& y){
298     if(b==0){
299         x=1,y=0;
300         return a;
301     }
302     else{
303         int gcd=ExGCD(b,a%b,y,x);
304         y-=x*(a/b);
305         return gcd;
306     }
307 }
308 
309 inline int Sqrt(int a,int p){
310     int s=Pow(G,Log(G,a,p)>>1,p);
311     return std::min(s,MOD-s);
312 }
313 
314 inline int Log(int a,int x,int p){
315     int s=sqrt(p+0.5);
316     int inv=Pow(Pow(a,s,p),p-2,p);
317     std::unordered_map<int,int> m;
318     m[1]=0;
319     int pow=1;
320     for(int i=1;i<s;i++){
321         pow=1ll*a*pow%p;
322         if(!m.count(pow))
323             m[pow]=i;
324     }
325     for(int i=0;i<s;i++){
326         if(m.count(x))
327             return i*s+m[x];
328         x=1ll*x*inv%MOD;
329     }
330     return -1;
331 }