洛谷P4199 万径人踪灭（manacher+FFT）

传送门

 

题目所求为所有的不连续回文子序列个数，可以转化为回文子序列数-回文子串数

回文子串manacher跑一跑就行了，考虑怎么求回文子序列数

我们考虑，如果$S_i$是回文子序列的对称中心，那么只要$S_{i-j}$和$S_{i+j}$相等，我们就多了一种选择

设共有$x$组相等的，那么以$S_i$为对称中心的回文子序列个数就是$2^{x+1}-1$，表示这$x$组包括对称中心选或不选，除去全都不选的都能算入答案

然而对称中心不一定在字符上可能在两个字符中间，那么这种时候回文子序列数就是$2^x-1$（因为没有中间的字符所以无所谓选不选）

然后考虑如何计算每一个位置上的$x$

我们考虑构造多项式$A$，原串上为$a$的位置设为$1$，$b$的位置设为$0$，如果$s[i]==s[j]$，那么他们的贡献会加到$(i+j)/2$上

然后发现这玩意儿和卷积很像，于是我们把除以二去掉，那么每一对$s[i]==s[j]$都会把贡献加到$i+j$上，所以只要把$A$自乘一下就可以了

然后构造$B$，原串上为$b$的位置设为$0$，也自乘一下就行了

然后把$A$和$B$对应系数加起来，再减去((i&1)^1)（表示这一位对称中心是否在字符上），再求一下$2$的多少次幂减一，最后减掉回文串个数就行了

 1 //minamoto
 2 #include<iostream>
 3 #include<cstdio>
 4 #include<cstring>
 5 #include<cmath>
 6 #include<algorithm>
 7 #define ll long long
 8 #define add(x,y) (x+y>=P?x+y-P:x+y)
 9 using namespace std;
10 const int N=5e5+5,P=1e9+7;const double Pi=acos(-1.0);
11 struct complex{
12     double x,y;
13     complex(double xx=0,double yy=0){x=xx,y=yy;}
14     inline complex operator +(complex b){return complex(x+b.x,y+b.y);}
15     inline complex operator -(complex b){return complex(x-b.x,y-b.y);}
16     inline complex operator *(complex b){return complex(x*b.x-y*b.y,x*b.y+y*b.x);}
17 }A[N],B[N];
18 int n,m,l,r[N],limit=1;
19 void FFT(complex *A,int type){
20     for(int i=0;i<limit;++i)
21     if(i<r[i]) swap(A[i],A[r[i]]);
22     for(int mid=1;mid<limit;mid<<=1){
23         complex Wn(cos(Pi/mid),type*sin(Pi/mid));
24         for(int R=mid<<1,j=0;j<limit;j+=R){
25             complex w(1,0);
26             for(int k=0;k<mid;++k,w=w*Wn){
27                 complex x=A[j+k],y=w*A[j+k+mid];
28                 A[j+k]=x+y,A[j+k+mid]=x-y;
29             }
30         }
31     }
32 }
33 ll s[N],ans[N],x[N],p[N],sum;char ss[N];
34 void manacher(){
35     int mx=0,id;
36     for(int i=1,l=(n<<1)+1;i<=l;++i){
37         p[i]=i<mx?min(p[2*id-i],(ll)mx-i):1ll;
38         while(x[i-p[i]]==x[i+p[i]]) ++p[i];
39         if(mx<i+p[i]) mx=i+p[i],id=i;
40     }
41 }
42 ll ksm(ll a,ll b){
43     ll res=1;
44     while(b){
45         if(b&1) res=res*a%P;
46         a=a*a%P,b>>=1;
47     }
48     return res;
49 }
50 int main(){
51 //    freopen("testdata.in","r",stdin);
52     scanf("%s",ss+1);n=strlen(ss+1);
53     for(int i=1;i<=n;++i) s[i]=ss[i]=='a';
54     for(int i=1;i<=n;++i)
55     x[(i<<1)-1]=2,x[i<<1]=s[i];
56     x[0]=-1,x[(n+1)<<1]=-2,x[(n<<1)+1]=2;
57     while(limit<=(n<<1)+1) limit<<=1,++l;
58     for(int i=0;i<limit;++i)
59     r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
60     for(int i=1;i<=n;++i) A[i].x=B[i].x=s[i];
61     FFT(A,1),FFT(B,1);
62     for(int i=0;i<limit;++i) A[i]=A[i]*B[i];
63     FFT(A,-1);
64     for(int i=1,l=(n<<1)+1;i<=l;++i) ans[i]+=((ll)(A[i].x/limit+0.5)-((i&1)^1));
65     memset(A,0,sizeof(A)),memset(B,0,sizeof(B));
66     for(int i=1;i<=n;++i) A[i].x=B[i].x=(s[i]^1);
67     FFT(A,1),FFT(B,1);
68     for(int i=0;i<limit;++i) A[i]=A[i]*B[i];
69     FFT(A,-1);
70     for(int i=1,l=(n<<1)+1;i<=l;++i) ans[i]+=((ll)(A[i].x/limit+0.5)-((i&1)^1));
71     for(int i=1,l=(n<<1)+1;i<=l;++i) ans[i]=((ans[i]+((i&1)^1))>>1)+((i&1)^1);
72     for(int i=1,l=(n<<1)+1;i<=l;++i) ans[i]=ksm(2,ans[i])-1;
73     manacher();
74     for(int i=1,l=(n<<1)+1;i<=l;++i) ans[i]-=(p[i]>>1);
75     for(int i=1,l=(n<<1)+1;i<=l;++i) sum=add(sum,ans[i]);
76     printf("%lld\n",sum);
77     return 0;
78 }

 

转载于:https://www.cnblogs.com/bztMinamoto/p/9759777.html