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

时间：2022-12-26 11:50:26

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;

}

相关文章