文章目录

洛谷 P3763 TJOI2017 DNA 做题笔记

发布于

题意

加里敦大学的生物研究所,发现了决定人喜不喜欢吃藕的基因序列 $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;
}

暂无评论

发表评论