万径人踪灭

Posted 2021-02-22 troverld

VII.万径人踪灭

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

跑题了跑题了

我们可以知道，

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

关于回文子串数，我们可以使用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^{leftlfloor{frac{h(x)+1}{2}} ight floor}-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]=‘‘;
}
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
",res);
	return 0;
}