题意
加里敦大学的生物研究所,发现了决定人喜不喜欢吃藕的基因序列 $T$,有这个序列的碱基序列就会表现出喜欢吃藕的性状,但是研究人员发现对碱基序列 $S$,任意修改其中不超过 $3$ 个碱基,依然能够表现出吃藕的性状。现在研究人员想知道这个基因在 DNA 链 $S$ 上的位置。所以你需要统计在一个表现出吃藕性状的人的 DNA 序列 $S$ 上,有多少个连续子串可能是该基因,即有多少个 $S$ 的连续子串修改小于等于三个字母能够变成 $T$。
$S, T$ 的长度不超过 $10^5$,$1 \leq T\leq 10$。
DNA 碱基序列只有 ATCG 四种字符。
题解
首先暴力枚举 $S$ 上的一个位置,跟 $T$ 做三次 LCP,看能不能匹配到 $T$ 的结尾即可。
还有一种 FFT 做法:
拿样例举例:
$S: \tt{ATCGCCCTA}$
$T: \tt{CTTCA}$
将 $S$ 和 $T$ 中为 $\tt C$ 的位置变为 $\tt1$,不为 $\tt C$ 的位置变为 $\tt0$ 再进行匹配。
$S: \tt{001011100}$
$T: \tt{10010}$
那么对于某个位置 $k$,$S_{k, \ldots, k + m - 1}$ 和 $T_{1, \ldots, m}$ 的匹配个数可以写为:
$$\sum\limits_{i = 1}^m S_{i + k} T_i$$
设 $A_{m - i - 1} = T_i$,
$$\sum\limits_{i + j = m - k - 1}^m S_i A_j$$
上式是卷积的形式,可以 FFT 求解。最后判断是否匹配可以有 $3$ 的误差。
代码
#include <cstdint>
#include <iostream>
#include <string>
using namespace std;
using tp = int64_t;
constexpr tp _BASE_ = 7, HASHAT_N = 2e5 + 3;
tp AL, BL;
string A, B;
uintmax_t HASH[HASHAT_N], PowTable[HASHAT_N];
uintmax_t HASHash(tp l, tp r) {
return HASH[r] - HASH[l - 1] * PowTable[r - l + 1];
}
tp LCP(tp a, tp b) {
tp l = 0, r = min(A.size() - a, B.size() - b);
while (l < r) {
tp mid = l + r >> 1;
if (HASHash(a, a + mid) == HASHash(AL + b, AL + b + mid)) {
l = mid + 1;
} else {
r = mid;
}
}
return l;
}
void Pre() {
*PowTable = 1;
for (tp i = 1; i < HASHAT_N; ++i) {
PowTable[i] = PowTable[i - 1] * _BASE_;
}
}
signed main() {
tp __TEST_CASE__;
cin >> __TEST_CASE__;
Pre();
while (__TEST_CASE__--) {
tp cnt = 0;
cin >> A >> B;
AL = A.size();
BL = B.size();
A = '\114' + A;
B = '\114' + B;
*HASH = 0;
for (tp i = 1; i < A.size(); ++i) {
HASH[i] = HASH[i - 1] * _BASE_ + A[i];
}
for (tp i = 1; i <= BL; ++i) {
HASH[AL + i] = HASH[AL + i - 1] * _BASE_ + B[i];
}
for (tp i = 1; i + BL <= A.size(); ++i) {
for (tp j = 0, a = i, b = 1; j < 4; ++j) {
tp len = LCP(a, b) + (j != 3);
a += len;
b += len;
if (b > BL) {
++cnt;
break;
}
}
}
cout << cnt << '\n';
}
return 0;
}