![]()

公共回文子串:求两个串的公共回文子串个数。
4.7华丽地处理字符串
后缀数组
拼接后先预处理出后缀数组、高度数组和回文半径数组。其中,每个位置的回文半径用Manacher算法线性时间求出。
定义a[plen]为半径=plen的公共回文子串的个数,下面看怎么逐步求a[plen]。按照后缀字典序遍历公共子串,维护目前遇到过的包含该公共子串的回文子串。
假设回文串长度为奇数,比如例题的遍历顺序如下(后缀|公共子串|回文串):
$CPCPC | | CPC$CPC C | | C C$CPCPC | C | C CPC | C | CPCPC CPC$CPCPC | CPC | C CPCPC | CPC | C ICPC$CPCPC | | I PC | | CPC PC$CPCPC | PC | CPC PCPC | PC | CPC
在每遇到一个长clen公共子串c的过程中,都遍历长度plen大于clen的公共回文子串p,将p的个数加到区间[clen, plen]中的每一个元素上去。因为每出现一次长plen的公共回文子串p,都代表出现了一次plen-1,plen-2, … clen的公共回文子串,示意图如下:

涉及到区域批量修改,所以要用BIT。
对每个位置k,都将区域和a[1,plen]累积到答案上即可。
#include <iostream>
#include <string>
#include <cstring>
#include <map>
#include <algorithm>
using namespace std;
typedef long long LL;
struct SA
{
static const int N = 200000; // len+1<=N
int a[N], lcp[N];
void build(const string &s)
{ build(s.c_str(), s.size()); }
void build(const char *t, int n)
{
// suffix array
static int g[N], b[N];
for (int i = 0; i < n + 1; ++i)
{
a[i] = i, g[i] = t[i];
}
b[0] = b[n] = 0;
sort(a, a + n + 1, SAComp(0, g));
for (int h = 1; b[n] != n; h *= 2)
{
SAComp comp(h, g);
sort(a, a + n + 1, comp);
for (int i = 0; i < n; ++i)
{
b[i + 1] = b[i] + comp(a[i], a[i + 1]);
}
for (int i = 0; i < n + 1; ++i)
{
g[a[i]] = b[i];
}
}
// longest common prefix (you can omit here)
for (int i = 0; i < n + 1; ++i)
{
b[a[i]] = i;
}
int h = 0;
for (int i = 0; i < n; ++i)
{
const int j = a[b[i] - 1];
while (j + h < n && i + h < n && t[j + h] == t[i + h])
h++;
lcp[b[i]] = h;
if (h)
h--;
}
lcp[0] = -1;
}
struct SAComp
{
const int h, *g;
SAComp(int h, int *g) : h(h), g(g)
{}
bool operator()(int l, int r)
{
return g[l] != g[r] ? g[l] < g[r] : g[l + h] < g[r + h];
}
};
} sa;
struct Manacher
{
static const int N = 200000 * 2; // len*2
int rad[N];
void build(const string &s)
{ build(s.c_str(), s.size()); }
void build(const char *t, int n)
{
int k;
for (int i = 0, j = 0; i < 2 * n; i += k, j = max(j - k, 0))
{
while (i - j >= 0 && i + j + 1 < 2 * n && t[(i - j) / 2] == t[(i + j + 1) / 2])
j++;
rad[i] = j;
for (k = 1; i - k >= 0 && rad[i] - k >= 0 && rad[i - k] != rad[i] - k; k++)
{
rad[i + k] = min(rad[i - k], rad[i] - k);
}
}
}
} mana;
struct BIT
{
static const int N = 80000;
LL x[N];
void init()
{ memset(x, 0, sizeof(x)); }
void add(int k, LL a)
{
for (; k < N; k |= k + 1)
x[k] += a;
}
LL sum(int k)
{
LL s = 0;
for (; k >= 0; k = (k & (k + 1)) - 1)
s += x[k];
return s;
}
LL rsum(int i, int j)
{ return i == 0 ? sum(j) : sum(j) - sum(i - 1); }
LL esum(int i)
{ return rsum(i, N - 1); }
} sum[2], cnt[2];
int of[200000];
int main()
{
string S, T;
cin >> S >> T;
const string F(S + '$' + T);
sa.build(F);
mana.build(F);
for (int i = 0; i < F.size() + 1; ++i)
{
of[i] = sa.a[i] > (int) S.size();
}
LL ans = 0;
for (int odd = 0; odd < 2; ++odd)// 回文长度是否为奇数
{
for (int i = 0; i < 2; ++i)
sum[i].init(), cnt[i].init();
map<int, int> palindrome[2];// AB各自的回文子串半径plen的个数
for (int k = 1; k < (int) F.size() + 1; k++)
{
const int clen = sa.lcp[k];// 公共子串的长度
for (int i = 0; i < 2; ++i)
{
map<int, int>::iterator it = palindrome[i].lower_bound(clen + 1);
while (it != palindrome[i].end())// 遍历半径大于公共子串clen的回文子串
{//维护Binary Indexed Tree
const int plen = it->first, p = it->second;// 回文子串长度及个数
sum[i].add(plen, -(LL) plen * p);//在[clen,plen]区间同时加上p
cnt[i].add(plen, -p);
sum[i].add(clen, (LL) clen * p);
cnt[i].add(clen, p);
palindrome[i][clen] += p;
palindrome[i].erase(it++);
}
}// 在sa[k]位置的回文半径
const int plen = odd ? mana.rad[sa.a[k] * 2] / 2 + 1 : sa.a[k] ? mana.rad[sa.a[k] * 2 - 1] / 2 : 0;
ans += sum[1 - of[k]].sum(plen) + (LL) plen * cnt[1 - of[k]].esum(plen + 1);//前plen项之和
sum[of[k]].add(plen, plen); // a[plen] += 1
cnt[of[k]].add(plen, 1);
palindrome[of[k]][plen]++;
}
}
cout << ans << endl;
return 0;
}

Reference
http://d.hatena.ne.jp/kohyatoh/
知识共享署名-非商业性使用-相同方式共享:码农场 » AOJ 2292 Common Palindromes 题解《挑战程序设计竞赛》
码农场