1. 程式人生 > >BZOJ.3160.萬徑人蹤滅(FFT Manacher)

BZOJ.3160.萬徑人蹤滅(FFT Manacher)

BZOJ
洛谷

\(Description\)

給定一個只含\(a,b\)的字串,求不連續迴文子序列個數(不連續指子序列不是連續一段,迴文要求字元和位置都關於某條對稱軸對稱)。
\(n\leq10^5\)

\(Solution\)

不連續的限制比較麻煩。連續迴文子序列個數可以用\(manacher\)求出,所以可以去掉不連續這個限制。

若對稱軸在\(x\)處,假設它的兩邊共有\(f(x)\)對對稱的字元,那麼以\(x\)為中心的非空迴文子序列數有\(2^{f(x)}-1\)個。因為\(x\)是可以取\(\frac k2\)的(字元之間的位置做對稱軸),用\(2x\)表示對稱軸為\(x\)

,答案就是\(\sum_{i=2}^{2n}(2^{f(i)}-1)\)

考慮如何求\(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;
}