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

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define N (400009)
#define LL long long
#define MOD (1000000007)
using namespace std;

int n,fn,l,tot,r[N],len[N],p[N];
LL Re,fun;
char s[N],st[N];
double pi=acos(-1.0);
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0)
    {
        x=xx; y=yy;
    }
}a[N],b[N];

complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
complex operator / (complex a,double b) {return complex(a.x/b,a.y/b);}

void FFT(int n,complex *a,int opt)
{
    for (int i=0; i<n; ++i)
        if (i<r[i]) swap(a[i],a[r[i]]);
    for (int k=1; k<n; k<<=1)
    {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=k<<1)
        {
            complex w=complex(1,0);
            for (int j=0; j<k; ++j,w=w*wn)
            {
                complex x=a[i+j], y=w*a[i+j+k];
                a[i+j]=x+y; a[i+j+k]=x-y;
            }
        }
    }
    if (opt==-1) for (int i=0; i<n; ++i) a[i]=a[i]/n;
}

void Manacher()
{
    s[++tot]='('; s[++tot]='#';
    for (int i=0; i<n; ++i)
        s[++tot]=st[i], s[++tot]='#';
    s[++tot]=')';
    int maxn=0,mid=0,x;
    for (int i=1; i<=tot; ++i)
    {
        if (i>maxn) x=1;
        else x=min(maxn-i+1,len[mid*2-i]);
        while (s[i+x]==s[i-x]) x++;
        len[i]=x;
        if (i+x-1>maxn) maxn=i+x-1, mid=i;
        fun=(fun+len[i]/2)%MOD;
    }
}

int main()
{
    p[0]=1;
    for (int i=1; i<=100000; ++i)
        p[i]=p[i-1]*2%MOD;
    scanf("%s",st);    n=strlen(st);
    Manacher();

    fn=1;
    while (fn<=n+n) fn<<=1, l++;
    for (int i=0; i<fn; ++i)
        r[i]=(r[i>>1]>>1) | ((i&1)<<(l-1));
    for (int i=0; i<n; ++i)
        if (st[i]=='a') a[i].x=1;
        else b[i].x=1;
    FFT(fn,a,1); FFT(fn,b,1);
    for (int i=0; i<fn; ++i)
        a[i]=a[i]*a[i], b[i]=b[i]*b[i];
    FFT(fn,a,-1); FFT(fn,b,-1);
    for (int i=0; i<fn; ++i)
    {
        int x=(a[i].x+b[i].x+0.5);
        x=(x+1)>>1;
        Re=(Re+p[x]-1)%MOD;
    }
    printf("%lld
",(Re-fun+MOD)%MOD);
}