BZOJ3160:万径人踪灭(FFT,Manacher)
Solution
$ans=$回文子序列$-$回文子串的数目。
后者可以用$manacher$直接求。
前者设$f[i]$表示以$i$为中心的对称的字母对数。
那么回文子序列的数量也就是$\sum_{i=0}^{n-1}2^{f[i]-1}$
构造两个数组$a[i],b[i]$。若第$i$位为$a$，那么$a[i]=1$，否则$b[i]=1$。
可以发现$a$数组自身卷积就是$a$字母对$f$数组的贡献，$b$数组同理。
卷下$a$，卷下$b$，对应位置求和，就是$f$数组。
因为在卷积中每对对称字符被算了两次，而自己和自己关于自己对称只算了一次，所以要把答案除2向上取整。
Code
1 #include<iostream> 2 #include<cstring> 3 #include<cstdio> 4 #include<cmath> 5 #define N (400009) 6 #define LL long long 7 #define MOD (1000000007) 8 using namespace std; 9 10 int n,fn,l,tot,r[N],len[N],p[N]; 11 LL Re,fun; 12 char s[N],st[N]; 13 double pi=acos(-1.0); 14 struct complex 15 { 16 double x,y; 17 complex (double xx=0,double yy=0) 18 { 19 x=xx; y=yy; 20 } 21 }a[N],b[N]; 22 23 complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);} 24 complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);} 25 complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} 26 complex operator / (complex a,double b) {return complex(a.x/b,a.y/b);} 27 28 void FFT(int n,complex *a,int opt) 29 { 30 for (int i=0; i<n; ++i) 31 if (i<r[i]) swap(a[i],a[r[i]]); 32 for (int k=1; k<n; k<<=1) 33 { 34 complex wn=complex(cos(pi/k),opt*sin(pi/k)); 35 for (int i=0; i<n; i+=k<<1) 36 { 37 complex w=complex(1,0); 38 for (int j=0; j<k; ++j,w=w*wn) 39 { 40 complex x=a[i+j], y=w*a[i+j+k]; 41 a[i+j]=x+y; a[i+j+k]=x-y; 42 } 43 } 44 } 45 if (opt==-1) for (int i=0; i<n; ++i) a[i]=a[i]/n; 46 } 47 48 void Manacher() 49 { 50 s[++tot]='('; s[++tot]='#'; 51 for (int i=0; i<n; ++i) 52 s[++tot]=st[i], s[++tot]='#'; 53 s[++tot]=')'; 54 int maxn=0,mid=0,x; 55 for (int i=1; i<=tot; ++i) 56 { 57 if (i>maxn) x=1; 58 else x=min(maxn-i+1,len[mid*2-i]); 59 while (s[i+x]==s[i-x]) x++; 60 len[i]=x; 61 if (i+x-1>maxn) maxn=i+x-1, mid=i; 62 fun=(fun+len[i]/2)%MOD; 63 } 64 } 65 66 int main() 67 { 68 p[0]=1; 69 for (int i=1; i<=100000; ++i) 70 p[i]=p[i-1]*2%MOD; 71 scanf("%s",st); n=strlen(st); 72 Manacher(); 73 74 fn=1; 75 while (fn<=n+n) fn<<=1, l++; 76 for (int i=0; i<fn; ++i) 77 r[i]=(r[i>>1]>>1) | ((i&1)<<(l-1)); 78 for (int i=0; i<n; ++i) 79 if (st[i]=='a') a[i].x=1; 80 else b[i].x=1; 81 FFT(fn,a,1); FFT(fn,b,1); 82 for (int i=0; i<fn; ++i) 83 a[i]=a[i]*a[i], b[i]=b[i]*b[i]; 84 FFT(fn,a,-1); FFT(fn,b,-1); 85 for (int i=0; i<fn; ++i) 86 { 87 int x=(a[i].x+b[i].x+0.5); 88 x=(x+1)>>1; 89 Re=(Re+p[x]-1)%MOD; 90 } 91 printf("%lld\n",(Re-fun+MOD)%MOD); 92 }
