[BZOJ3160]万径人踪灭-快速傅里叶变换-Manacher算法

万径人踪灭

Description

Input

Output

Sample Input

Sample Output

HINT

---

千山鸟飞绝
思路特别好的一道题~学习一个

---

思路:
可以发现这个“不连续的回文子序列”不好求。
那么考虑容斥，用“所有的回文子序列”减去“所有的回文子串”即可。

后者很显然是个Manacher。
前者的话，考虑设$f[i]$表示以$i$为对称中心的对称字符对数。
则因为只有选或不选两种情况，答案为$\sum 2^{f[i]}-1$。

考虑如何求$f[i]$:
有一个很明显的性质是，所有关于同一位置对称的字符对的下标和均为相同的值。
所以在进行一次卷积后，关于同一位置对称的字符对的贡献会落入同一个下标，也就是同一个$f[i]$中
先令所有$a$为1,$b$为0，做一遍FFT。
再令所有$b$为1,$a$为0，做一遍FFT。
两边FFT的和+1后除以2就是 f <script type="math/tex" id="MathJax-Element-829">f</script>数组。

那这就做完了~

#include<iostream>
#include<cmath>
#include<cstdio>
#include<complex>
#include<cstring>
#include<cstdlib>
#include<algorithm>

using namespace std;

typedef long long ll;
typedef double db;
typedef complex<db> cp;
const int N=400009;
const ll md=1e9+7;
const db pi=3.1415926535897;

int n,m,l;
int p[N],rev[N];
ll f[N],ans;
cp a[N],b[N];
char s[N],ch[N];

inline ll qpow(ll a,ll b)
{
    ll ret=1ll;
    while(b)
    {
        if(b&1ll)
            ret=ret*a%md;
        a=a*a%md;
        b>>=1;
    }
    return ret;
}

inline void FFT(cp *a,int n,int f)
{
    for(int i=0;i<n;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int h=1;h<n;h<<=1)
    {
        cp w(cos(pi/h),f*sin(pi/h));
        for(int j=0;j<n;j+=(h<<1))
        {
            cp wn(1,0);
            for(int k=j;k<j+h;k++)
            {
                cp x=a[k],y=wn*a[k+h];
                a[k]=x+y;
                a[k+h]=x-y;
                wn*=w;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<n;i++)
            a[i].real()/=(db)n;
}

int main()
{
    scanf("%s",s);
    n=strlen(s);

    for(int i=0;i<n;i++)
    {
        ch[i<<1]='#';
        ch[i<<1|1]=s[i];
    }
    m=n<<1;
    ch[m]='#';

    p[0]=1;
    int mx=0,mp=0;
    for(int i=1;i<=m;i++)
    {
        if(i<=mx)
            p[i]=min(p[2*mp-i],mx-i+1);
        while(i+p[i]<=m && i-p[i]>=0 && ch[i+p[i]]==ch[i-p[i]])
            p[i]++;
        if(i+p[i]-1>mx)
            mx=i+p[i]-1,mp=i;
    }

    for(int i=0;i<=m;i++)
        ans=(ans-(p[i]>>1)+md)%md;

    for(m=1,l=0;m<=(n<<1);m<<=1)l++;
    for(int i=0;i<m;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<n;i++)
        a[i]=(s[i]=='a'),b[i]=(s[i]=='b');

    FFT(a,m,1);FFT(b,m,1);
    for(int i=0;i<m;i++)
        a[i]*=a[i],b[i]*=b[i];
    FFT(a,m,-1);FFT(b,m,-1);

    for(int i=0;i<m;i++)
        f[i]=(int)(a[i].real()+0.5)+(int)(b[i].real()+0.5);
    for(int i=0;i<m;i++)
        ans=(ans+qpow(2,(f[i]+1)>>1)-1+md)%md;

    printf("%lld\n",ans);
    return 0;
}