BZOJ.3160.萬徑人蹤滅(FFT Manacher)
阿新 • • 發佈:2018-11-28
\(Description\)
給定一個只含\(a,b\)的字串,求不連續迴文子序列個數(不連續指子序列不是連續一段,迴文要求字元和位置都關於某條對稱軸對稱)。
\(n\leq10^5\)。
\(Solution\)
不連續的限制比較麻煩。連續迴文子序列個數可以用\(manacher\)求出,所以可以去掉不連續這個限制。
若對稱軸在\(x\)處,假設它的兩邊共有\(f(x)\)對對稱的字元,那麼以\(x\)為中心的非空迴文子序列數有\(2^{f(x)}-1\)個。因為\(x\)是可以取\(\frac k2\)的(字元之間的位置做對稱軸),用\(2x\)表示對稱軸為\(x\)
考慮如何求\(f(x)\)。
一對關於\(x\)對稱的字元,就是滿足\(s_{x-i}=s_{x+i}\)。因為\(x\)可以取\(\frac12\),不妨換種表示,即\(s_{i}=s_{2x-i}\)。
那麼\(f(x)=\frac{\sum_{i=1}^{2x-1} [s_{i}=s_{2x-i}]+1}{2}\)(每對對稱字元算了兩次,而自己關於自己對稱的字元算了一次,所以加\(1\)再除以\(2\))。而那個\(\sum\)就是卷積的形式。
列舉一種字元\(c\),令\(F_i=[s_i=c]\),那麼\(f(x)=\frac{\sum_{i=1}^{2x-1}F_i*F_{2x-i}+1}{2}\)
所以列舉一下字元\(a,b\)然後\(FFT\),加起來就可以求出\(\sum_{i=1}^{2x-1}F_i*F_{2x-i}\)了。
兩次計算\(\sum\)的時候中間可以不\(IDFT\)。
//12340kb 1524ms #include <cmath> #include <cstdio> #include <cstring> #include <algorithm> #define gc() getchar() #define mod 1000000007 typedef long long LL; const int N=(1<<18)+7; const double PI=acos(-1); int pw[N],rev[N]; char s[N]; struct Complex { double x,y; Complex(double x=0,double y=0):x(x),y(y) {} Complex operator +(const Complex &a)const {return Complex(x+a.x, y+a.y);} Complex operator -(const Complex &a) {return Complex(x-a.x, y-a.y);} Complex operator *(const Complex &a) {return Complex(x*a.x-y*a.y, x*a.y+y*a.x);} }A[N],f[N]; void FFT(Complex *a,int lim,int opt) { for(int i=1; i<lim; ++i) if(i<rev[i]) std::swap(a[i],a[rev[i]]); for(int i=2; i<=lim; i<<=1) { int mid=i>>1; Complex Wn(cos(PI/mid),opt*sin(PI/mid)); for(int j=0; j<lim; j+=i) { Complex w(1,0),t; for(int k=j; k<j+mid; ++k,w=w*Wn) a[k+mid]=a[k]-(t=a[k+mid]*w), a[k]=a[k]+t; } } if(opt==-1) for(int i=0; i<lim; ++i) a[i].x/=lim; } void Solve(const char c,int n,int m,int lim) { for(int i=0; i<lim; ++i) A[i]=Complex(0,0); for(int i=1; i<=n; ++i) A[i]=Complex(s[i]==c,0); FFT(A,lim,1); for(int i=0; i<lim; ++i) f[i]=f[i]+A[i]*A[i];//0~lim! // for(int i=0; i<lim; ++i) A[i]=A[i]*A[i]; // FFT(A,lim,-1); // for(int i=1; i<=m; ++i) f[i]+=(int)(A[i].x+0.5); } int Manacher(int n) { static int ext[N]; s[0]='$', s[n+n+1]='#', s[n+n+2]='%'; for(int i=n; i; --i) s[i<<1]=s[i], s[(i<<1)-1]='#'; n<<=1; for(int i=1,mx=0,p; i<=n; ++i) { ext[i]=mx>i?std::min(mx-i,ext[2*p-i]):1; while(s[i-ext[i]]==s[i+ext[i]]) ++ext[i]; if(i+ext[i]>mx) mx=i+ext[i], p=i; } LL ans=0; for(int i=1; i<=n; ++i) ans+=ext[i]>>1; return ans%mod; } int main() { scanf("%s",s+1); int n=strlen(s+1),m=n<<1,lim=1,l=-1; while(lim<=m) lim<<=1,++l; for(int i=1; i<lim; ++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<l); Solve('a',n,m,lim), Solve('b',n,m,lim); FFT(f,lim,-1); LL ans=0; pw[0]=1; for(int i=1; i<=m; ++i) pw[i]=pw[i-1]<<1, pw[i]>=mod&&(pw[i]-=mod); for(int i=1; i<=m; ++i) ans+=pw[(int)(f[i].x+1.5)>>1]-1; // for(int i=1; i<=m; ++i) f[i]=f[i]+1>>1, ans+=pw[f[i]]-1; printf("%lld\n",(ans%mod-Manacher(n)+mod)%mod); return 0; }