E - DNA配列のパターン検索 / Pattern Search in DNA Sequences 解説 by admin
gemini-3.5-flash-thinking概要
本問題は、与えられた基準パターン \(P\) の一部を反転させた変形パターン \(P'_j\) が、指定された文字列 \(S_{i_j}\) の中に何回出現するかを、複数のクエリに対して高速に求める問題です。
文字列の長さやクエリ数が非常に大きいため、各クエリに対して愚直に変形パターンを作成して検索を行うのではなく、ローリングハッシュ（Rolling Hash）と二分探索を組み合わせることで、クエリあたり \(O(\log |S_{i_j}|)\) で高速に処理します。
考察
素朴なアプローチとその限界
各クエリ \((i_j, l_j, r_j)\) に対して、実際に \(P\) の区間 \([l_j, r_j]\) を反転して \(P'_j\) を作成し、KMP法やローリングハッシュ等で \(S_{i_j}\) から \(P'_j\) を検索すると、1クエリあたり \(O(|P| + |S_{i_j}|)\) の時間がかかります。 クエリ数 \(Q \le 2 \times 10^5\)、文字列の長さが最大 \(5 \times 10^5\) であるため、最悪の場合全体の計算量が \(O(Q (|P| + |S_i|))\) となり、実行時間制限（TLE）になってしまいます。
高速化のアイデア
クエリごとに検索対象の文字列 \(S_{i_j}\) と変形パターン \(P'_j\) は変化しますが、以下の2つの点に注目します。
\(S_i\) の部分文字列ハッシュの事前計算 各 \(S_i\) について、長さ \(|P|\) のすべての連続部分文字列のハッシュ値を事前に計算し、ソートして保持しておきます。もし変形パターン \(P'_j\) のハッシュ値 \(HP'_j\) を高速に計算できれば、ソート済みのハッシュ値配列に対して二分探索（
std::equal_rangeなど）を行うことで、出現回数を高速にカウントできます。変形パターン \(P'_j\) のハッシュ値の高速合成 変形パターン \(P'_j\) は、元の基準パターン \(P\) を3つの部分に分割したものとして表せます（0-indexed で考えます。\(l' = l_j - 1, r' = r_j - 1\) とします）。
- 左側（反転しない）: \(P[0 \dots l'-1]\)
- 中央（反転する）: \(P[l' \dots r']\) の反転
- 右側（反転しない）: \(P[r'+1 \dots |P|-1]\)
これら3つの部分のハッシュ値がそれぞれ \(O(1)\) で求まれば、それらを結合した全体のハッシュ値 \(HP'_j\) も \(O(1)\) で計算可能です。
アルゴリズム
1. ローリングハッシュの準備
ハッシュの衝突（ハッシュ値が同じなのに異なる文字列である現象）を防ぐため、非常に大きな素数 \(MOD = 2^{61}-1\) を法とし、ランダムに決定した基数 \(B\) を用いたローリングハッシュを使用します。
2. 事前テーブルの構築
- 基準パターン \(P\) およびそれを反転した文字列 \(P_{rev}\) の累積ハッシュ配列を構築します。これにより、任意の区間のハッシュ値を \(O(1)\) で取得できるようになります。
- 各 \(S_i\) について、長さ \(|P|\) のすべての連続部分文字列のハッシュ値をスライドしながら計算し、配列
V[i]に格納して昇順にソートします。
3. クエリの処理（ハッシュ値の合成）
各クエリについて、変形パターン \(P'_j\) のハッシュ値 \(HP'_j\) を以下のように合成します。
- 左側のハッシュ \(H_A\): \(P\) の区間 \([0, l'-1]\) のハッシュ値。元の \(P\) の累積ハッシュから \(O(1)\) で取得。
- 中央のハッシュ \(H_B\): \(P\) の区間 \([l', r']\) を反転したもののハッシュ値。これは反転文字列 \(P_{rev}\) における区間 \([|P|-1-r', |P|-1-l']\) のハッシュ値に等しいため、\(P_{rev}\) の累積ハッシュから \(O(1)\) で取得。
- 右側のハッシュ \(H_C\): \(P\) の区間 \([r'+1, |P|-1]\) のハッシュ値。元の \(P\) の累積ハッシュから \(O(1)\) で取得。
これらを適切な基数 \(B\) の冪乗を掛けて足し合わせることで、全体のハッシュ値を合成します。 $\(HP'_j = H_A \times B^{|P| - l'} + H_B \times B^{|P| - 1 - r'} + H_C \pmod{MOD}\)$
4. 二分探索によるカウント
合成したハッシュ値 \(HP'_j\) を用いて、ソート済みの配列 V[i] に対して二分探索を行い、ハッシュ値が一致する要素の個数を \(O(\log |S_i|)\) で求めます。
具体例でのハッシュ合成
\(P = \) HMLLM (\(|P| = 5\)), \(l_j = 2, r_j = 4\) の場合（0-indexedで \(l' = 1, r' = 3\)）
- 元パターン:
H M L L M - 反転パターン \(P_{rev}\):
M L L M H - 変形パターン \(P'_j\):
H L L M M
- 左側 \(A\) (
H): \(P[0 \dots 0]\) のハッシュ \(\rightarrow H_A\) - 中央 \(B\) (
LLM): \(P[1 \dots 3]\) の反転。これは \(P_{rev}[5-1-3 \dots 5-1-1] = P_{rev}[1 \dots 3]\) （LLM）のハッシュ \(\rightarrow H_B\) - 右側 \(C\) (
M): \(P[4 \dots 4]\) のハッシュ \(\rightarrow H_C\)
これらを結合して \(P'_j\) のハッシュ値を得ます。
計算量
時間計算量
事前準備:
- \(P\) および \(P_{rev}\) のハッシュ構築: \(O(|P|)\)
- 全ての \(S_i\) のハッシュ構築および部分文字列ハッシュのソート: 各 \(S_i\) のソートに \(O(|S_i| \log |S_i|)\) かかるため、全体で \(O(\sum_{i=1}^{N} |S_i| \log |S_i|)\)
クエリ処理:
- 1クエリあたり、ハッシュの合成に \(O(1)\)、二分探索に \(O(\log |S_{i_j}|)\)
- 全 \(Q\) クエリで \(O(Q \log (\max |S_i|))\)
総合時間計算量: \(O(|P| + \sum_{i=1}^{N} |S_i| \log |S_i| + Q \log (\max |S_i|))\) 制約の \(\sum |S_i| \le 5 \times 10^5\)、 \(Q \le 2 \times 10^5\) に対し、実行時間制限内に余裕で間に合います。
空間計算量
- \(P, P_{rev}\) のハッシュテーブル、および各 \(S_i\) の部分文字列ハッシュを保持する配列
V - 総合空間計算量: \(O(|P| + \sum_{i=1}^{N} |S_i|)\) メモリ制限に対しても非常に省メモリで動作します。
実装のポイント
61-bit Mersenne素数ハッシュ: ハッシュの法（MOD）として \(2^{61}-1\) を使用することで、ハッシュ衝突の確率を実用上ゼロに抑えられます。また、この法に対する剰余算はビット演算を用いて高速化できるため、定数倍が非常に軽くなります。
基数 \(B\) のランダム化: ハッシュ衝突を意図的に引き起こすような「ハックケース」を防ぐため、プログラムの実行時に乱数で基数 \(B\) を決定しています。
ソースコード
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <random>
#include <chrono>
using namespace std;
using ull = unsigned long long;
const ull MOD = (1ULL << 61) - 1;
ull mul(ull a, ull b) {
__uint128_t t = (__uint128_t)a * b;
ull lo = (ull)t & MOD;
ull hi = (ull)(t >> 61);
ull res = lo + hi;
if (res >= MOD) res -= MOD;
return res;
}
ull add(ull a, ull b) {
ull res = a + b;
if (res >= MOD) res -= MOD;
return res;
}
ull sub(ull a, ull b) {
ull res = a + MOD - b;
if (res >= MOD) res -= MOD;
return res;
}
ull get_base() {
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
uniform_int_distribution<ull> dist(300, MOD - 2);
ull b = dist(rng);
if (b % 2 == 0) b++;
return b;
}
vector<ull> powerB;
void init_power(ull B, int max_len) {
powerB.resize(max_len + 1);
powerB[0] = 1;
for (int i = 1; i <= max_len; i++) {
powerB[i] = mul(powerB[i - 1], B);
}
}
vector<ull> build_hash(const string& S, ull B) {
int n = S.length();
vector<ull> h(n + 1, 0);
for (int i = 0; i < n; i++) {
h[i + 1] = add(mul(h[i], B), S[i]);
}
return h;
}
ull get_hash(const vector<ull>& h, int l, int r) {
if (l > r) return 0;
int len = r - l + 1;
return sub(h[r + 1], mul(h[l], powerB[len]));
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int N, Q;
if (!(cin >> N >> Q)) return 0;
string P;
cin >> P;
int M = P.length();
vector<string> S(N);
int max_len = M;
for (int i = 0; i < N; i++) {
cin >> S[i];
max_len = max(max_len, (int)S[i].length());
}
ull B = get_base();
init_power(B, max_len);
vector<ull> hp = build_hash(P, B);
string P_rev = P;
reverse(P_rev.begin(), P_rev.end());
vector<ull> hprev = build_hash(P_rev, B);
vector<vector<ull>> V(N);
for (int i = 0; i < N; i++) {
int len = S[i].length();
if (len >= M) {
vector<ull> hs = build_hash(S[i], B);
V[i].resize(len - M + 1);
for (int a = 0; a <= len - M; a++) {
V[i][a] = get_hash(hs, a, a + M - 1);
}
sort(V[i].begin(), V[i].end());
}
}
for (int q = 0; q < Q; q++) {
int i, l, r;
cin >> i >> l >> r;
i--;
if (V[i].empty()) {
cout << 0 << "\n";
continue;
}
int l_prime = l - 1;
int r_prime = r - 1;
ull HA = get_hash(hp, 0, l_prime - 1);
ull HB = get_hash(hprev, M - 1 - r_prime, M - 1 - l_prime);
ull HC = get_hash(hp, r_prime + 1, M - 1);
ull HP_prime = add(mul(HA, powerB[M - l_prime]), mul(HB, powerB[M - 1 - r_prime]));
HP_prime = add(HP_prime, HC);
auto range = equal_range(V[i].begin(), V[i].end(), HP_prime);
cout << distance(range.first, range.second) << "\n";
}
return 0;
}
この解説は gemini-3.5-flash-thinking によって生成されました。
投稿日時:
最終更新: