1. 程式人生 > >萬徑人蹤滅(FFT+manacher)

萬徑人蹤滅(FFT+manacher)

傳送門

這題……我覺得像我這樣的菜雞選手難以想出來……

題目要求求出一些子序列,使得其關於某個位置是對稱的,而且不能是連續一段,求這樣的子序列的個數。這個直接求很困難,但是我們可以先求出所有關於某個位置對稱的子序列,最後減去子串的個數。

子串個數可以用\(manacher\)求,至於子序列的話,我們假設以第\(i\)位為中心,那麼如果兩邊有\(x\)對相同的字元,那麼這個位置對答案的貢獻就是\(2^x-1\)或者\(2^(x+1)-1\)。(因為有可能迴文串的長度是偶數,也就是不存在中間點)

考慮怎麼求\(x_i\)\(x_i\)的形式可以寫成如下的形式:
\[\sum_{j=0}^ic[i-j] == c[i+j]\]

發現這個式子非常像卷積的形式。那麼我們先初始化兩個序列,第一個序列是原字串為‘a’,對應位置為1,第二個是原字串為'b',對應位置是1,剩下都是0。這樣結果就轉化為如下形式:

\[\sum_{j=0}^i a(i-j) * a(i+j) + b(i-j) * b(i+j) \]

然後讓他們自己和自己乘起來,結果相加一下,然後因為卷積會重複把元素計算兩遍,所以要+1再/2.

這樣得到的各項係數就是各項\(x_i\),我們就可以用快速冪計算。算完之後減去\(manacher\)求出的子串個數即可。

看一下程式碼。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('\n')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back
#define I puts("bug")

using namespace std;
typedef long long ll;
const int M = 200005;
const int INF = 1000000009;
const double eps = 1e-7;
const double pi = acos(-1);
const ll mod = 1e9+7;

int read()
{
    int ans = 0,op = 1;char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
    return ans * op;
}

struct Comp
{
   double x,y;
   Comp(){}
   Comp(double kx,double ky){x = kx,y = ky;}
   fr Comp operator + (const Comp &a,const Comp &b) {return (Comp){a.x + b.x,a.y + b.y};}
   fr Comp operator - (const Comp &a,const Comp &b) {return (Comp){a.x - b.x,a.y - b.y};}
   fr Comp operator * (const Comp &a,const Comp &b) {return (Comp){a.x * b.x - a.y * b.y,a.y * b.x + a.x * b.y};}
}a[M<<1],b[M<<1],kx,ky;

int n,len = 1,L,p[M<<1],rev[M<<1];
char s[M<<1],c[M];
ll tot,d[M<<1],ans;

int change()
{
   int l = strlen(c),j = 2;
   s[0] = '!',s[1] = '#';
   rep(i,0,l-1) s[j++] = c[i],s[j++] = '#';
   s[j] = '&';
   return j;
}

void manacher()
{
   int l = change(),mx = 1,mid = 1;
   rep(i,1,l-1)
   {
      if(i < mx) p[i] = min(mx - i,p[(mid<<1) - i]);
      else p[i] = 1;
      while(s[i-p[i]] == s[i+p[i]]) p[i]++;
      if(mx < i + p[i]) mid = i,mx = i + p[i];
      tot += (p[i] >> 1),tot %= mod;
   }
}

void fft(Comp *a,int f)
{
   rep(i,0,len-1) if(i < rev[i]) swap(a[i],a[rev[i]]);
   for(int i = 1;i < len;i <<= 1)
   {
      Comp omg1(cos(pi / i),f * sin(pi / i));
      for(int j = 0;j < len;j += (i << 1))
      {
     Comp omg(1,0);
     rep(k,0,i-1)
     {
        kx = a[k+j],ky = omg * a[k+j+i];
        a[k+j] = kx + ky,a[k+j+i] = kx - ky,omg = omg * omg1;
     }
      }
   }
}

ll qpow(ll a,ll b)
{
   ll p = 1;
   while(b)
   {
      if(b & 1) p *= a,p %= mod;
      a *= a,a %= mod;
      b >>= 1;
   }
   return p;
}

int main()
{
   scanf("%s",c);
   n = strlen(c);
   rep(i,0,n-1)
   {
      if(c[i] == 'a') a[i].x = 1;
      else b[i].x = 1;
   }
   while(len <= n << 1) len <<= 1,L++;
   //I;
   rep(i,0,len-1) rev[i] = (rev[i>>1] >> 1) | ((i&1) << (L-1));
   fft(a,1),fft(b,1);
   rep(i,0,len-1) a[i] = a[i] * a[i] + b[i] * b[i];
   fft(a,-1);
   rep(i,0,len-1) d[i] = ((ll)floor(a[i].x / len + 0.5) + 1) >> 1;
   rep(i,0,len-1) ans += (qpow(2,d[i]) - 1),ans %= mod;
   manacher(),ans -= tot,ans %= mod;
   while(ans < 0) ans += mod;
   printf("%lld\n",ans);
   return 0;
}