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