万径人踪灭

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;
}