BZOJ 3160: 万径人踪灭 [fft manacher]

3160: 万径人踪灭

题意：求一个序列有多少不连续的回文子序列

---

一开始zz了直接用(2^{r_i}-1)

总-回文子串

后者用manacher处理

前者，考虑回文有两种对称形式(以元素/缝隙作为对称轴)

f[i]，i为奇数表示以缝隙对称，偶数表示以元素i>>1对称，对答案的贡献就是(2^{f[i]}-1)

[f[i] = sum_{j=1}^{i-1} [s_j = s_{i-j}] ]

就是裸卷积

因为只有a，b两个字符，可以先后令a或b=1分别求

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define fir first
#define sec second
typedef long long ll;
const int N=(1<<19)+5, mo=1e9+7;
const double PI = acos(-1);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}
    return x*f;
}

struct meow{
	double x, y;
	meow(double a=0, double b=0):x(a), y(b){}
};
meow operator +(meow a, meow b) {return meow(a.x+b.x, a.y+b.y);}
meow operator -(meow a, meow b) {return meow(a.x-b.x, a.y-b.y);}
meow operator *(meow a, meow b) {return meow(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);}
meow conj(meow a) {return meow(a.x, -a.y);}
typedef meow cd;

namespace fft {
	int n, rev[N];
	void ini(int lim) {
		n=1; int k=0; while(n<lim) n<<=1, k++;
		for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));
	}
	void dft(cd *a, int flag) {
		for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
		for(int l=2; l<=n; l<<=1) {
			int m = l>>1;
			cd wn = cd(cos(2*PI/l), flag * sin(2*PI/l));
			for(cd *p=a; p != a+n; p+=l) {
				cd w(1, 0);
				for(int k=0; k<m; k++) {
					cd t = w * p[k+m];
					p[k+m] = p[k] - t;
					p[k] = p[k] + t;
					w = w*wn;
				}
			}
		}
		if(flag == -1) for(int i=0; i<n; i++) a[i].x /= n;
	}
	void mul(cd *a) {
		dft(a, 1);
		for(int i=0; i<n; i++) a[i] = a[i] * a[i];
		dft(a, -1);
	}
} 

cd a[N];
int n, f[N]; ll ans;
char s[N];
ll Pow(ll a, int b) {
	ll ans=1;
	for(; b; b>>=1, a=a*a%mo)
		if(b&1) ans=ans*a%mo;
	return ans;
}
inline void mod(ll &x) {if(x<0) x+=mo; else if(x>=mo) x-=mo;}

namespace ma {
	int r[N]; char a[N];
	ll manacher(char *s, int n) {
		int p=0, a; ll ans=0;
		for(int i=1; i<=n; i++) {
			r[i] = i<p ? min(p-i+1, r[2*a-i]) : 1;
			while(s[ i-r[i] ] == s[ i+r[i] ]) r[i]++;
			if(i+r[i]-1 > p) p = i+r[i]-1, a=i;
			mod(ans += r[i]>>1);
		}
		return ans;
	}
	ll cal(char *s) {
		for(int i=1; i<=n; i++) a[(i<<1)-1] = '#', a[i<<1] = s[i];
		a[(n<<1)+1] = '#';
		a[0] = '@'; a[(n<<1)+2] = '$';
		return manacher(a, (n<<1)+1);
	}
}

int main() {
	freopen("in","r",stdin);
	scanf("%s", s+1); n = strlen(s+1);
	fft::ini(n+n+1);
	for(int i=1; i<=n; i++) a[i].x = s[i]=='a';
	fft::mul(a);
	for(int i=1; i<=n+n; i++) f[i] = (int)floor(a[i].x + 0.5);

	for(int i=0; i<fft::n; i++) a[i] = cd(0, 0);
	for(int i=1; i<=n; i++) a[i].x = s[i]=='b';
	fft::mul(a);
	for(int i=1; i<=n+n; i++) f[i] += (int)floor(a[i].x + 0.5);
	for(int i=1; i<=n+n; i++) f[i] = (f[i]+1)>>1;
	//for(int i=1; i<=n+n; i++) printf("f %d  %d
", i, f[i]);

	for(int i=1; i<=n+n; i++) mod(ans += Pow(2, f[i]) - 1);
	
	//printf("ans1 %lld
", ans);
	mod(ans -= ma::cal(s));
	printf("%lld", ans);
}