BZOJ3160: 万径人踪灭 FFT+manacher

题意：给一个字符串，求（不能是连续的一段的）回文子序列的数量。要求回文子序列的字符和位置都必须以某一对称轴对称，就是说“a空格b空格空格a”这样的是不算的。
长度<=100000，只有a、b两种字符。
先考虑允许连续的情况下怎么求。不妨假设原串每两个字符间已经插入了‘#’。先研究每一个对称中心，对其有贡献的字符对满足在它一左一右，距离它距离相等，且字符值相等。如果已经统计出了一个中心两边符合条件的点对数x（含自身），则每个点对可以任意选或不选，一共有2^x-1种选法（不能一个都不选）。将每个中心的方法数累加，再用manacher减去连续的方案数即可。
所以现在要求的就是每个中心的点对数。显然一个中心所对应的点对坐标和为定值，因此可以用卷积来求。对于“字符值相等”，我们可以分别处理每种字符的贡献，先让所有’a’处为1，求自乘卷积，再让’b’处为1，求自乘卷积，二者相加。由于每个点对算了两遍（自身算了一遍），所以结果为（二者之和+1）/2。
manacher中，若最长延展长度为k，则有k/2个连续子串。

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#define gm 262200
using namespace std;
typedef long long ll;
const int mob=1000000007;
struct cpx
{
    double a,b;
    cpx(const double &a=0,const double &b=0):a(a),b(b){}
    #define c ano.a
    #define d ano.b
    inline cpx operator + (const cpx& ano) const
    {
        return cpx(a+c,b+d);
    }
    inline cpx operator - (const cpx& ano) const
    {
        return cpx(a-c,b-d);
    }
    inline cpx operator * (const cpx& ano) const
    {
        return cpx(a*c-b*d,b*c+a*d);
    }
    #undef c
    #undef d
}x[gm],res[gm];
char s[gm];
int pos[gm];
const double DFT=2.0,IDFT=-2.0;
inline void trans(cpx a[],const int len,const double& mode)
{
    static const double pi=acos(-1);
    for(int i=0;i<len;++i)
    if(i<pos[i]) swap(a[i],a[pos[i]]);
    for(int i=2;i<=len;i<<=1)
    {
        cpx wi(cos(2*pi/i),sin(mode*pi/i));
        int step=i>>1;
        for(int j=0;j<len;j+=i)
        {
            cpx w(1.0,0.0);
            int lim=j+step;
            for(int k=j;k<lim;++k)
            {
                cpx l=a[k],r=a[k+step]*w;
                a[k]=l+r,a[k+step]=l-r;
                w=w*wi;
            }
        }
    }
    if(mode==IDFT)
    for(int i=0;i<len;++i)
    a[i].a/=len;
}
inline void manacher(const char *__s,ll& ans)
{
    static int p[gm];
    static char s[gm]="$#";
    int top=1;
    while(*__s)
    {
        s[++top]=*__s;
        s[++top]='#';
        ++__s;
    }
    s[++top]='\0';
    int mx=0,id=0;
    for(int i=1;i<top;++i)
    {
        if(i<mx) p[i]=min(p[(id<<1)-i],mx-i);
        else p[i]=1;
        while(s[i-p[i]]==s[i+p[i]]) ++p[i];
        if(i+p[i]>mx) mx=i+p[i],id=i;
        ans-=p[i]>>1;
        if(ans<0) ans+=mob;
    }
}
#define pow __POW
int pow[gm]={1};
ll ans=0;
int main()
{
    scanf("%s",s);
    int len=strlen(s);
    int n=1;
    while(n<(len<<1)) n<<=1;
    for(int i=1;i<n;++i)
    {
        pos[i]=pos[i>>1]>>1;
        if(i&1) pos[i]|=n>>1;
    }
    for(int i=0;i<n;++i)
    x[i]=s[i]=='a'?cpx(1.0,0.0):cpx();
    trans(x,n,DFT);
    for(int i=0;i<n;++i)
    res[i]=x[i]*x[i];
    for(int i=0;i<n;++i)
    x[i]=s[i]=='b'?cpx(1.0,0.0):cpx();
    trans(x,n,DFT);
    for(int i=0;i<n;++i)
    res[i]=res[i]+x[i]*x[i];
    trans(res,n,IDFT);
    for(int i=1;i<n;++i)
    {
        pow[i]=pow[i-1]<<1;
        if(pow[i]>=mob) pow[i]-=mob;
    }
    for(int i=0;i<n;++i)
    {
        ans+=pow[ll(res[i].a+0.5)+1>>1]-1;
        if(ans>=mob) ans%=mob;
    }
    manacher(s,ans);
    printf("%lld\n",ans);
    return 0;
}