[bzoj3160] 万径人踪灭 [FFT+manacher]

题面

传送门

思路

首先我们明确一个想法：答案等于“位置对称的回文子序列数”-“回文子串数”，而后面的回文子串数可以通过$manacher$做出来

所以我们的问题转化为：求“位置对称的会文字序列个数”

我们考虑以字符串$S$的第$i$个字符$S_i$作为一个子序列的对称中心

那么当$S_{i-j}=S_{i+j}left(1leq jleq i-1 ight)$时，我们的答案就加一种选择

如果对于$S_i$，有$x$个这样的相等字符对，那么最终以第$i$个位置为中心的答案就是$2^{x+1}-1$

为什么呢？

因为算上对称中心本身，一共有$x+1$组可以选或者不选的字符，除去全都不选的一种就是答案

然而这是i为某个特定字符的情况，或者说i为整数的情况，但是i还有在两个字符中间的情况，此时答案就是$2^x-1$，因为没有中间的自己那个字符

接下来，我们考虑如何计算每个位置上这样的字符组数量

我们构建两个多项式$A$与$B$，令$A$的所有原串中为'a'的位置系数为1，'b'系数为零，$B$则反过来

然后用$A$自乘，得到的新多项式的第i位的值，就是对称中心为第$frac i2$的，两个字符都是'a'的字符组数，$B$同理

最后我们把$A$和$B$的每一位系数加起来，减去((i&1)^1)（代表对称中心是否在某个字符上），求2的次幂再减去1，最后减去manacher求出的回文子串数目就可以了

Code：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
const ll MOD=1e9+7;
struct complex{
    double x,y;
    complex(double xx=0,double yy=0){x=xx;y=yy;}
    complex operator +(const complex &b){return complex(x+b.x,y+b.y);}
    complex operator -(const complex &b){return complex(x-b.x,y-b.y);}
    complex operator *(const complex &b){return complex(x*b.x-y*b.y,x*b.y+y*b.x);}
}A[400010],B[400010];
const double pi=acos(-1.0);
ll n,m,limit=1,cnt=0,r[400010];
void fft(complex *a,double type){
    ll i,j,k,mid;complex x,y,w,wn;
    for(i=0;i<limit;i++) if(i<r[i]) swap(a[i],a[r[i]]);
    for(mid=1;mid<limit;mid<<=1){
        wn=complex(cos(pi/mid),type*sin(pi/mid));
        for(j=0;j<limit;j+=(mid<<1)){
            w=complex(1,0);
            for(k=0;k<mid;k++,w=w*wn){
                x=a[j+k];y=w*a[j+k+mid];
                a[j+k]=x+y;a[j+k+mid]=x-y;
            }
        }
    }
}
ll s[100010],ans[200010],x[200010],p[200010],sum;char ss[100010];
void manacher(){
    ll maxn=-1,id,mx=0,i;
    for(i=1;i<=(n<<1)+1;i++){
        if(i<mx) p[i]=min(p[2*id-i],mx-i);
        else p[i]=1;
        while(x[i-p[i]]==x[i+p[i]]) p[i]++;
        if(mx<i+p[i]){
            mx=p[i]+i;id=i;
        }
    }
}
ll qpow(ll x,ll y){
    ll re=1;
    while(y){
        if(y&1) re=re*x%MOD;
        x=x*x%MOD;y>>=1;
    }
    return re;
}
int main(){
    ll i;
    scanf("%s",ss);n=strlen(ss);
    for(i=0;i<n;i++) s[i+1]=(ss[i]=='a');
    for(i=1;i<=(n<<1)+1;i++){
        if(i%2) x[i]=2;
        else x[i]=s[i>>1];
    }x[0]=-1,x[(n+1)<<1]=-2;
    
    while(limit<=(n<<1)) limit<<=1,cnt++;
    for(i=0;i<limit;i++) r[i]=((r[i>>1]>>1)|((i&1)<<(cnt-1)));
    
    for(i=1;i<=n;i++) A[i].x=B[i].x=s[i];
    fft(A,1);fft(B,1);
    for(i=0;i<=limit;i++) A[i]=A[i]*B[i];
    fft(A,-1);
    for(i=1;i<=(n<<1)+1;i++) ans[i]+=((ll)(A[i].x/limit+0.5)-((i&1)^1));
    
    memset(A,0,sizeof(A));memset(B,0,sizeof(B));
    for(i=1;i<=n;i++) A[i].x=B[i].x=(s[i]^1);
    fft(A,1);fft(B,1);
    for(i=0;i<=limit;i++) A[i]=A[i]*B[i];
    fft(A,-1);
    for(i=1;i<=(n<<1)+1;i++) ans[i]+=((ll)(A[i].x/limit+0.5)-((i&1)^1));
    for(i=1;i<=(n<<1)+1;i++) ans[i]=((ans[i]+((i&1)^1))>>1)+((i&1)^1);
    
    for(i=1;i<=(n<<1)+1;i++) ans[i]=qpow(2,ans[i])-1,ans[i]%=MOD;
    
    manacher();
    
    for(i=1;i<=(n<<1)+1;i++) ans[i]-=(p[i]>>1),ans[i]%=MOD;
    for(i=1;i<=(n<<1)+1;i++) sum+=ans[i],sum%=MOD;
    printf("%lld",sum);
}