[SDOI2015]序列统计

Posted 2020-10-26 hyghb

[SDOI2015]序列统计

 

很有趣的一道题目，很巧妙。

显然是一个dp，考虑最朴素的dp，f[i][j]表示选i个乘起来，%m为j的方案数为多少。转移也很简单。

然而乘法的转移并不能进行什么优化，于是考虑设法将其转为加法。

我们可以通过求出m的原根，因为原根G，G^i %m(1<=i<m)是一一对应1<=x<m的所以我们将原本S集合中的元素由"x"替换为"i",这样就完成了由乘法向加法的转换，但是需要注意的是，此时的模数并不再是m而是m-1,。下文中的m均为输入的m-1后的值。

那么问题就转换为了从一堆数中选n个加起来，最终%m为j的方案数为多少，如果我们用g数组来表示某一个元素是否存在，我们发现f[i+1][j]是由f[i][k]*g[(j-k+m)%m](0<=k<m)转移而来，我们发现他很像卷积的形式，但是%m怎么处理呢？

这里get了一个新技能，好像叫循环卷积。

就是计算f和g的卷积后，对于m<=i<l这一段，将他的值加到i%m上去。

f数组自然是可以重复使用的。

那么我们重复卷积n次，得到的即为最终答案。

然而n太大了。

卷积是具有交换律和结合律的，所以可以用快速幂来加速这一过程。

整体复杂度mlogmlogn

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 const int inf=4e4+10;
 4 const int mod=1004535809;
 5 int n,m,x,s;
 6 int a[inf];
 7 int fast(int x,int y,int z){
 8     int ans=1;
 9     while(y){
10         if(y&1)ans=1ll*ans*x%z;
11         x=1ll*x*x%z;
12         y>>=1;
13     }
14     return ans;
15 }
16 namespace RT{
17     int tmp[inf],cnt;
18     bool check(int x,int y){
19         for(int i=1;i<=cnt;i++){
20             if(::fast(x,(y-1)/tmp[i],y)==1)return 0;
21         }
22         return 1;
23     }
24     int get_rt(int x){
25         int u=x;
26         x--;
27         for(int i=2;i*i<=x;i++){
28             if(x%i)continue;
29             tmp[++cnt]=i;
30             while(x%i==0)x/=i;
31         }
32         if(x>1)tmp[++cnt]=x;
33         int now=2;
34         while(1){
35             if(check(now,u))return now;
36             now++;
37         }
38     }
39 }
40 int hs[inf];
41 int g[inf],f[inf];
42 int r[inf];
43 void ntt(int *a,int l,int type){
44     for(int i=0;i<l;i++)
45         if(i<r[i])swap(a[i],a[r[i]]);
46     for(int i=2;i<=l;i<<=1){
47         int wn=fast(3,(type*(mod-1)/i+(mod-1))%(mod-1),mod);
48         for(int j=0;j<l;j+=i){
49             int w=1;
50             for(int k=j;k<j+i/2;k++,w=1ll*w*wn%mod){
51                 int u=a[k],v=1ll*w*a[k+i/2]%mod;
52                 a[k]=(u+v)%mod;
53                 a[k+i/2]=(u-v+mod)%mod;
54             }
55         }
56     }
57 }
58 int main()
59 {
60     freopen("sdoi2015_sequence.in","r",stdin);
61     freopen("sdoi2015_sequence.out","w",stdout);
62     scanf("%d%d%d%d",&n,&m,&x,&s);
63     for(int i=1;i<=s;i++)scanf("%d",&a[i]);
64     int G=RT::get_rt(m);
65     for(int i=1;i<m;i++)hs[fast(G,i,m)]=i;
66     m--;
67     for(int i=1;i<=s;i++)
68         if(a[i])g[hs[a[i]]%m]=1;
69     int l=1,h=0;
70     while(l<m*2-1)l<<=1,h++;
71     for(int i=0;i<l;i++)
72         r[i]=(r[i>>1]>>1)+((i&1)<<(h-1));
73     int inv=fast(l,mod-2,mod);
74     f[0]=1;
75     while(n){
76         ntt(g,l,1);
77         if(n&1){
78             ntt(f,l,1);
79             for(int i=0;i<l;i++)f[i]=1ll*f[i]*g[i]%mod;
80             ntt(f,l,-1);
81             for(int i=0;i<l;i++)f[i]=1ll*f[i]*inv%mod;
82             for(int i=m;i<l;i++)f[i%m]=(f[i%m]+f[i])%mod,f[i]=0;
83         }
84         for(int i=0;i<l;i++)g[i]=1ll*g[i]*g[i]%mod;
85         ntt(g,l,-1);
86         for(int i=0;i<l;i++)g[i]=1ll*g[i]*inv%mod;
87         for(int i=m;i<l;i++)g[i%m]=(g[i%m]+g[i])%mod,g[i]=0;
88         n>>=1;
89     }
90     printf("%d\n",f[hs[x]%m]);
91     return 0;
92 }

View Code