BZOJ 3160: 万径人踪灭(FFT+快速幂+manacher)

题目传送门

---

题解

题目大意就是求一个只有’a’和’b’的字符串中不连续的回文子序列的数量。所谓的回文子序列就是说这个子序列关于某个字符或某个间隙左右对称。

举个栗子：在babbab中，b_bb_b就是一个合法的子序列。

字符串只有’a’和’b’，直接设为0或1卷积一下两下。
于是位置和相同的字符对们就都被我们求出来了，回文子序列就是它们组成的。

怎么来的还是要引用tututu大神的证明（括号内是我的补充）：

先不考虑对称轴为间隙的情况

假设对于位置$x$，设$f[x]$为其左右各有$f[x]$个字符关于$x$对称（包括$x$），那么根据二项式定理(取多少个的情况加起来就是2的幂)，这个对答案的贡献为 $2^{f[x]}-1$（集合非空子集数）。

$f[x]$显然可以暴力求，设原串为$s$,下标从$0$开始，大概是这样的

$$f[x]=\sum_{i=0}^{x}[s[x−i]==s[x+i]]$$ (不考虑越界)

由于$x-i+x+i=2*x$，这就是一个果的卷积，跑一遍NTT（改为FFT）即可，对称轴为间隙也类似。

然后求出所有的回文子序列后，再用manacher算法(好久没打了)求出所有回文子串(即连续的回文子序列)，然后减去即可。中间求$2^{f[x]}-1$时用快速幂加速一下。其他的细节就是如何把求出来的东西搞一下之类的细节了，自行yy就行了(或者看我代码)。

这样问题就愉快地解决了，不过这题目叫“万径人踪灭”，差点把我给灭了，一开始FFT写错，2*PI写成2/PI，我居然没发现？？后来算答案上界搞成n，然而我搞得LJ方法超上界会算多！！从n改到(len<<1)后，样例终于过了。一交就是一个RE。又查啊查，发现“马拉车”拉错了，2*p-i打成了2*p+i，导致越界了，QAQ。。。

改完终于AC了，菜不成声啊。。。

总结一下：正着做卷积位置和相同的卷到一起，反着做卷积“位置”相同的卷到一起(这里的位置可以理解为相对位置、位置差之类的)。

---

代码

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <cstdlib>
#define MAXN 100010
#define Mod 1000000007

using namespace std;


typedef long long LL;
int len, n;
char str[MAXN], s[MAXN<<1];
const double PI = acos(-1.0);

struct Complex{
    double real, image;
    Complex() {}
    Complex(double _real, double _image){
      real = _real;
      image = _image;
    }
    friend Complex operator + (Complex A, Complex B){return Complex(A.real + B.real, A.image + B.image);}
    friend Complex operator - (Complex A, Complex B){return Complex(A.real - B.real, A.image - B.image);}
    friend Complex operator * (Complex A, Complex B){return Complex(A.real * B.real - A.image * B.image, A.real * B.image + A.image * B.real);}
}a[MAXN<<2], b[MAXN<<2], c1[MAXN<<2], c2[MAXN<<2];

int P[MAXN], pali[MAXN<<1];
LL ans;


void Get_P(){
    for(int i = 1, t = 1; i <= len; i++){
      if(t < i)  t <<= 1;
      P[i] = t;
    }
}

void Reverse(Complex *A){
    for(int i = 0; i < n-1; i++){
      int j = 0;
      for(int k = 1, tmp = i; k < n; k <<= 1, tmp >>= 1)
        j = ((j << 1) | (tmp & 1));
      if(j > i)  swap(A[i], A[j]);
    }
}

void FFT(Complex *A, int DFT){
    Reverse(A);
    for(int s = 1; (1<<s) <= n; s++){
      int m = (1 << s);
      Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m));
      for(int k = 0; k < n; k += m){
        Complex w = Complex(1, 0);
        for(int j = 0; j < (m>>1); j++){
          Complex u = A[k+j], t = w * A[k+j+(m>>1)];
          A[k+j] = u + t;
          A[k+j+(m>>1)] = u - t;
          w = w * wm;
        }
      }
    }
    if(DFT == -1)  for(int i = 0; i < n; i++)  A[i].real /= n, A[i].image /= n;
}

LL Pow(LL x, int cnt){
    LL res = 1LL;
    while(cnt){
      if(cnt & 1)  res = res * x % Mod;
      x = x * x % Mod;
      cnt >>= 1;
    }
    return res - 1LL;
}

void Manacher(){

    int l = 0;
    s[0] = '@';
    for(int i = 0; i < len; i++){
      s[++l] = str[i];
      if(i != len-1)  s[++l] = '#';
    }
    s[++l] = '$';

    pali[0] = pali[1] = 0;
    int p = 1;

    for(int i = 2; i <= l; i++){
      pali[i] = max(0, min(pali[2*p-i], p+pali[p]-i));
      while(s[i-pali[i]-1] == s[i+pali[i]+1])  pali[i] ++;
      if(i + pali[i] > p + pali[p])  p = i;
    }
    for(int i = 1; i < l; i++){
      if(s[i] != '#')  ans = (ans - LL(pali[i] / 2 + 1) + Mod) % Mod;
      else  ans = (ans - LL(pali[i] + 1) / 2 + Mod) % Mod;
    }
}


int main(){

    freopen("bzoj3160.in", "r", stdin);
    freopen("bzoj3160.out", "w", stdout);

    scanf("%s", &str);
    len = strlen(str);

    Get_P();
    n = P[len] << 1;

    for(int i = 0; i < len; i++)  a[i] = b[i] = Complex(str[i]=='a', 0);
    for(int i = len; i < n; i++)  a[i] = b[i] = Complex(0, 0);

    FFT(a, 1);
    FFT(b, 1);
    for(int i = 0; i < n; i++)  c1[i] = a[i] * b[i];
    FFT(c1, -1);


    for(int i = 0; i < len; i++)  a[i] = b[i] = Complex(str[i]=='b', 0);
    for(int i = len; i < n; i++)  a[i] = b[i] = Complex(0, 0);

    FFT(a, 1);
    FFT(b, 1);
    for(int i = 0; i < n; i++)  c2[i] = a[i] * b[i];
    FFT(c2, -1);


    for(int i = 0; i < (len<<1); i++){
      int v = (int)(c1[i].real+0.5) + (int)(c2[i].real+0.5);
      v = (v >> 1) + (!(i & 1));
      ans = (ans + Pow(2LL, v)) % Mod;
    }

    Manacher();

    printf("%lld\n", ans);

    return 0;
}

---

我只想站在比你高的地方，用人类最纯粹的痛苦与烦恼给你一记响亮的耳光。

——《阴火》