[BZOJ3160]万径人踪灭 FFT+manacher

所求可以转化为所有位置和字母都回文的子序列减去回文子串。后者显然可以用manacher解决。
假设字符串倍增后，关于i对称的字母对一共有f[i]对，那么以其为对称中心的回文子序列就一共有2^f[i]-1种。我们还发现f[i]=∑[s[j]==s[i-j]]{1<=j<=i-1}，s为未倍增的原串。这是一个卷积的形式，可以用fft解决，分成∑[s[j]==s[i-j]==’a’]和∑[s[j]==s[i-j]==’b’]两部分dft，结果相加在idft即可，(这样求出来的是2*f[i]-1，因为除了这一位本身被算了一次，其它的都被算了两次)。
代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<complex>
#define ll long long 
using namespace std;
typedef complex<double> cd;
const double PI=acos(-1);
const int mod=1000000007;
char cs[100010],s[200010];
int n,p[200010],r[280010];
cd f[280010],a[280010];
ll ksm(ll a,int b)
{
    ll re=1;
    while(b){if(b&1) re=re*a%mod;a=a*a%mod;b>>=1;}
    return re;
}
void manacher(int n)
{
    int rz=0,mxr=0;
    for(int i=1;i<=n;i++)
    {
        if(i>mxr) {p[i]=0;while(s[i-p[i]]==s[i+p[i]]) p[i]++;}
        else
        {
            if(p[2*rz-i]<mxr-i+1) {p[i]=p[2*rz-i];continue;}
            else {p[i]=mxr-i;while(s[i-p[i]]==s[i+p[i]]) p[i]++;}
        }
        rz=i;mxr=i+p[i]-1;
    }
}
void getrev(int bit)
{
    for(int i=0;i<(1<<bit);i++)
        r[i]=(r[i>>1]>>1)+((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft)
{
    for(int i=0;i<n;i++)
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<n;i<<=1)
    {
        cd wn=exp(cd(0,dft*PI/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            cd wnk=cd(1,0);
            for(int k=j;k<i+j;k++)
            {
                cd x=a[k];
                cd y=wnk*a[k+i];
                a[k]=x+y;
                a[k+i]=x-y;
                wnk*=wn;
            }
        }
    }
    if(dft==-1) for(int i=0;i<n;i++) a[i]/=n;
}
int main()
{
    scanf("%s",cs+1);
    n=strlen(cs+1);
    for(int i=1;i<=n;i++)
        {s[2*i-1]='#';s[2*i]=cs[i];}
    s[0]='$';s[2*n+2]='&';s[2*n+1]='#'; 
    manacher(2*n+1);
    int bit=0;
    while((1<<bit)<2*n+1) bit++;
    getrev(bit);
    for(int i=0;i<(1<<bit);i++)
        if(i<=n&&cs[i]=='a') a[i]=1; else a[i]=0;

    fft(a,(1<<bit),1);
    for(int i=0;i<(1<<bit);i++)
        f[i]=a[i]*a[i];
    for(int i=0;i<(1<<bit);i++)
        if(i<=n&&cs[i]=='b') a[i]=1;else a[i]=0;

    fft(a,(1<<bit),1);

    for(int i=0;i<(1<<bit);i++)
        f[i]+=a[i]*a[i];
    fft(f,(1<<bit),-1);
    ll ans=0;
    for(int i=0;i<(1<<bit);i++)
        ans=(ans+ksm(2,((ll)(round(f[i].real()))+1)/2)-1+mod)%mod;
    for(int i=1;i<=2*n+1;i++)
        ans=(ans-p[i]/2+mod)%mod;           
    printf("%lld",ans);
    return 0;
}