BZOJ 3160: 万径人踪灭 FFT+快速幂+manacher

BZOJ 3160: 万径人踪灭

题目传送门
【题目大意】
　　给定一个长度为n的01串，求有多少个回文子序列？
　　回文子序列是指从原串中找出任意个，使得构成一个回文串，并且位置也是沿某一对称轴对称。

　　假如x是对称轴，若 i 和 j 是对称且di=dj，i,j可以视为可行的一组。可行组数记为f[x]。
　　$f[x]=\sum_{i=1}^{x-1}[d[x-i]==d[x+i]]$
　　以x为对称轴的答案是2^(f[x])-1。
　　可以观察发现将d[i]=1的A[i]标为1，A与A做一次卷积，即可得出d[i]=1的f[]。（因为$p^{x-i}*p^{x+i}=p^x$）
　　对d[i]=1的做一次卷积，对d[i]=0的做一次卷积，加起来就是完整的f[]。
　　由于单个连续一段的回文串不算，做一次manacher，减掉就行了。（ps：感觉这是强行加上去的条件，增加码量……）
　　正解：FFT+快速幂+manacher

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>

#define imax(a,b) ((a>b)?(a):(b))
#define imin(a,b) ((a<b)?(a):(b))

using namespace std;

typedef long long ll;

const ll mods=1000000007;
const int N=101000;
const double pi=acos(-1.0);
char st[N];
int n,d[N<<1],dn;
int nn,k,F[N<<1];
ll sum,tp[N<<2],power[N<<2];
struct Complex
{
    double real,image;
    Complex() {}
    Complex(double _real,double _image)
    {
        real=_real; image=_image;
    }
    friend Complex operator + (Complex A,Complex B) { return Complex(A.real+B.real,A.image+B.image); }
    friend Complex operator - (Complex A,Complex B) { return Complex(A.real-B.real,A.image-B.image); }
    friend Complex operator * (Complex A,Complex B) { return Complex(A.real*B.real-A.image*B.image,A.image*B.real+A.real*B.image); }
} A[N<<2],B[N<<2];
int rev[N<<2];

void FFT(Complex *A,int n,int DFT)
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]);
    for(int s=1;(1<<s)<=n;s++)
    {
        int mi=(1<<s);
        Complex wn=Complex(cos(2*pi/mi),DFT*sin(2*pi/mi));
        for(int t=0;t<n;t+=mi)
        {
            Complex w=Complex(1,0);
            for(int j=0;j<(mi>>1);j++)
            {
                Complex u=A[t+j],v=w*A[t+j+(mi>>1)];
                A[t+j]=u+v; A[t+j+(mi>>1)]=u-v;
                w=w*wn;
            }
        }
    }
    if(DFT==-1) for(int i=0;i<n;i++) A[i].real/=n,A[i].image/=n;
}

ll manacher()
{
    F[1]=0; int p=1; ll re=0ll;
    for(int i=2;i<=dn;i++)
    {
        F[i]=imax(0,imin(F[2*p-i],F[p]+p-i));
        while(d[i+F[i]+1]==d[i-F[i]-1]) F[i]++;
        if(F[i]+i>F[p]+p) p=i;
        re=(re+((F[i]+1)>>1))%mods;
    }
    return re;
}

ll Pow(ll x,ll y)
{
    ll s=1ll;
    for(;y;y>>=1,x=x*x%mods) s=s*x%mods;
    return s;
}

int main()
{
    scanf("%s",st); n=strlen(st);
    dn=1; d[1]=4;
    for(int i=1;i<=n;i++)
    {
        d[++dn]=3;
        d[++dn]=st[i-1]-'a';
    } d[++dn]=3; d[++dn]=5;

    sum=(mods-manacher())%mods;

    int m=n<<1; k=0;
    for(nn=1;nn<=m;nn<<=1) k++;
    for(int i=0;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));

    power[0]=1ll;
    for(int i=1;i<=nn;i++) power[i]=(power[i-1]<<1)%mods;

    for(int i=1;i<=n;i++) if(st[i-1]=='a') A[i].real=1;
    FFT(A,nn,1); for(int i=0;i<=nn;i++) B[i]=A[i]*A[i];

    memset(A,0,sizeof(A));
    for(int i=1;i<=n;i++) if(st[i-1]=='b') A[i].real=1;
    FFT(A,nn,1); for(int i=0;i<=nn;i++) B[i]=B[i]+A[i]*A[i]; FFT(B,nn,-1);

    for(int i=1;i<nn;i++) tp[i]=ll(B[i].real+0.5);
    //for(int i=1;i<nn;i++) printf("%lld\n",tp[i]);
    for(int i=1;i<nn;i++) (sum+=power[(tp[i]+1)>>1]-1)%=mods;

    printf("%lld\n",(sum%mods+mods)%mods);
    return 0;
}