万径人踪灭

VII.万径人踪灭

千山鸟飞绝，万径人踪灭。
孤舟蓑笠翁，独钓寒江雪。
  ——【唐】柳宗元《江雪》

跑题了跑题了

我们可以知道，

\(\text{答案=位置对称且字符对称的子序列的数量-回文子串数}\)

关于回文子串数，我们可以使用Manacher算法在\(O(n)\)时间内实现。如果不会的话，可以参加鄙人的拙作，这里不再赘述。

那么如何求出位置对称且字符对称的子序列的数量呢？

考虑如果对于某个对称中心来说，有\(k\)对字符关于它位置对称并且字符对称，这个对称中心的贡献是多少。

当这个对称中心是一个具体的字符时，贡献为\(2^{k+1}-1\)，这\(k\)对字符，连带它自己，一共\(k+1\)个物品，都有选与不选两种情况，共\(2^{k+1}\)中方案。当然，不能一个物品都不选，因此要减去\(\varnothing\)时的方案，即减一。

当这个对称中心是两个字符中间的夹缝时，贡献为\(2^k-1\)。理由同上。

考虑将\(a\)和\(b\)分开考虑。以考虑\(a\)为例。则我们令\(f(x)=[s_x='a']\)，用\(f\)自乘。我们看看结果有何含义。

当\(f(x)=f(y)=1\)时，\(f(x+y)\)会加一。而它的实际含义是，\(\dfrac{x+y}{2}\)处的\(k\)增加了\(1\)。也就是说，\(f(i)\)的实际含义是\(\dfrac{i}{2}\)处的\(k\)值。

考虑\(b\)时同理，我们令\(g(x)=[s_x='b']\)，再用\(g\)自乘。之后，将\(f\)和\(g\)加在一起，得到每个位置的\(k\)。

但这个\(k\)并不是真实的\(k\)。在卷积时，\(f(x)f(y)\)与\(f(y)f(x)\)，除非\(x=y\)，否则都被考虑了两次。因此，如果\(x\)是一个字符，\(h(x)\)必定为奇，因为除了它自己被考虑一次以外，其他的对都被考虑两次。如果\(x\)是一个夹缝，\(h(x)\)必定为偶，因为所有对都被考虑两次。

这样的话，我们可以令答案增加\(\large2^{\left\lfloor{\frac{h(x)+1}{2}}\right\rfloor}-1\)，当\(h(x)\)为奇时，实际上增加了\(1\)；当\(h(x)\)为偶数时，实际上什么也没有增加。

代码：

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1);
const int mod=1e9+7;
int S,lim=1,lg,rev[1<<20],h[1<<20],pov[1<<20],res,rad[1<<20];
char s[1<<20];
struct cp{
	double x,y;
	cp(double u=0,double v=0){x=u,y=v;}
	friend cp operator +(const cp &u,const cp &v){return cp(u.x+v.x,u.y+v.y);}
	friend cp operator -(const cp &u,const cp &v){return cp(u.x-v.x,u.y-v.y);}
	friend cp operator *(const cp &u,const cp &v){return cp(u.x*v.x-u.y*v.y,u.x*v.y+u.y*v.x);}
}f[1<<20],g[1<<20];
void FFT(cp *a,int tp){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int md=1;md<lim;md<<=1){
		cp rt=cp(cos(pi/md),tp*sin(pi/md));
		for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
			cp w=cp(1,0);
			for(int i=0;i<md;i++,w=w*rt){
				cp x=a[pos+i],y=w*a[pos+md+i];
				a[pos+i]=x+y;
				a[pos+md+i]=x-y;
			}
		}
	}
}
void Prepare(){
	for(int i=S-1;i>=0;i--)s[2*i+1]=s[i];
	S=2*S+1;
	for(int i=0;i<S;i+=2)s[i]='#';
	s[S]='\0';
}
void Manacher(){
	Prepare();
	int Rpos=-1,Centre=-1;
	for(int i=0;i<S;i++){
		rad[i]=(i<Rpos)?min(Rpos-i,rad[2*Centre-i]):1;
		while(i-rad[i]>=0&&i+rad[i]<S)if(s[i-rad[i]]==s[i+rad[i]])rad[i]++;else break;
		if(i+rad[i]>Rpos)Rpos=i+rad[i],Centre=i;
		res-=rad[i]>>1;
		if(res<0)res+=mod;
	}
}
int main(){
	scanf("%s",s),S=strlen(s);
	while(lim<=S*2)lim<<=1,lg++;
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	pov[0]=1;
	for(int i=1;i<lim;i++)pov[i]=(pov[i-1]<<1)%mod;
	for(int i=0;i<S;i++)f[i]=cp(s[i]=='a',0),g[i]=cp(s[i]=='b',0);
	FFT(f,1),FFT(g,1);
	for(int i=0;i<lim;i++)f[i]=f[i]*f[i],g[i]=g[i]*g[i];
	FFT(f,-1),FFT(g,-1);
	for(int i=0;i<lim;i++)h[i]=(int)(f[i].x/lim+0.5)+(int)(g[i].x/lim+0.5),res=(res+pov[(h[i]+1)>>1]-1)%mod;
	Manacher(); 
	printf("%d\n",res);
	return 0;
}