Day 051-Q1 | Suffix Array + BWT + FM-Index（全文検索・パターン出現位置列挙）

問題

長さ $N$ の文字列 $S$（小文字英字のみ）と $Q$ 個のクエリが与えられる。各クエリは文字列 $P_i$ であり、$S$ 中に $P_i$ が出現する全ての開始位置（0-indexed）を昇順に出力せよ。

制約

パラメータ	範囲
$N$	$1 \le N \le 2 \times 10^5$
$Q$	$1 \le Q \le 10^4$
クエリ長の総和	$\le 10^5$
文字種	小文字英字のみ
時間制限	3秒

入出力例

入力例 1

13 3
abracadabra$
3
abr
1
a
5
racad

出力例 1

0 7
0 3 5 7 10
3

"abr" は位置 0 と 7 に出現。"a" は 0,3,5,7,10 に出現。"racad" は 3 のみ。

概念図: FM-Index の LF-Mapping

ヒント（段階的開示）

ヒント1: 方向性

Suffix Array だけでパターン検索は $O(|P| \log N)$（二分探索）だが、FM-Index の LF-Mapping を使うと $O(|P|)$ で SA 上の範囲を特定できる。BWT（Burrows-Wheeler Transform）が鍵。

ヒント2: アプローチ

SA-IS: Suffix Array を $O(N)$ で構築
BWT: $BWT[i] = S[SA[i]-1]$（SA[i]=0 なら番兵 '$'）
C[c]: c より辞書順小の文字の BWT 中の総数
Occ(c, i): BWT[0..i-1] での c の出現数（チェックポイントで高速化）
LF-Mapping: $P$ を右から読み `lo = C[c] + Occ(c, lo)` で範囲を絞り込む

ヒント3: LF-Mapping の実装骨格

def fm_search(bwt, C, occ_func, SA, pattern):
    lo, hi = 0, len(bwt)
    for c in reversed(pattern):
        if c not in C:
            return []
        lo = C[c] + occ_func(c, lo)
        hi = C[c] + occ_func(c, hi)
        if lo >= hi:
            return []
    return sorted(SA[lo:hi])

チェックポイント間隔 $B$ で Occ を $O(B)$ に削減。$B = 128$ が実用的。

模範解答 (Python)

import sys
from collections import defaultdict
input = sys.stdin.readline

def build_sa_is(s_int, sigma):
    n = len(s_int)
    if n == 1: return [0]
    if n == 2: return [0,1] if s_int[0] <= s_int[1] else [1,0]
    t = [False]*n
    for i in range(n-2,-1,-1):
        t[i] = (s_int[i] > s_int[i+1]) or (s_int[i]==s_int[i+1] and t[i+1])
    def is_lms(i): return i>0 and not t[i] and t[i-1]
    lms_pos = [i for i in range(1,n) if is_lms(i)]
    def induced_sort(lms_sorted):
        sa = [-1]*n
        cnt = [0]*(sigma+1)
        for c in s_int: cnt[c] += 1
        end = [0]*(sigma+1)
        acc = 0
        for i in range(sigma+1):
            acc += cnt[i]; end[i] = acc
        start = [end[i]-cnt[i] for i in range(sigma+1)]
        ends = list(end)
        for p in reversed(lms_sorted):
            c = s_int[p]; ends[c] -= 1; sa[ends[c]] = p
        starts = list(start)
        for i in range(n):
            if sa[i]>0 and t[sa[i]-1]:
                c = s_int[sa[i]-1]; sa[starts[c]] = sa[i]-1; starts[c] += 1
        ends = list(end)
        for i in range(n-1,-1,-1):
            if sa[i]>0 and not t[sa[i]-1]:
                c = s_int[sa[i]-1]; ends[c] -= 1; sa[ends[c]] = sa[i]-1
        return sa
    sa = induced_sort(lms_pos)
    lms_order = [p for p in sa if is_lms(p)]
    rank = [-1]*n; cur = -1; prev = -1
    for p in lms_order:
        cur += 1; rank[p] = cur
    lms_ranks = [rank[p] for p in lms_pos]
    if len(set(lms_ranks)) < len(lms_pos):
        rec_sa = build_sa_is(lms_ranks, max(lms_ranks)+1)
        lms_sorted = [lms_pos[i] for i in rec_sa]
    else:
        lms_sorted = [0]*len(lms_pos)
        for i,p in enumerate(lms_pos): lms_sorted[rank[p]] = p
    return induced_sort(lms_sorted)

def solve():
    N, Q = map(int, input().split())
    S = input().strip()
    s = S + '$'; n = len(s)
    chars = sorted(set(s))
    c2i = {c:i for i,c in enumerate(chars)}
    s_int = [c2i[c] for c in s]
    SA = build_sa_is(s_int, len(chars))
    BWT = [s[SA[i]-1] if SA[i]>0 else '$' for i in range(n)]
    freq = defaultdict(int)
    for c in BWT: freq[c] += 1
    sorted_chars = sorted(freq.keys())
    C = {}; acc = 0
    for c in sorted_chars: C[c] = acc; acc += freq[c]
    B = 128; chk = (n//B)+2
    checkpoints = {c:[0]*chk for c in sorted_chars}
    for c in sorted_chars:
        cnt = 0
        for i in range(n):
            if i%B==0: checkpoints[c][i//B] = cnt
            if BWT[i]==c: cnt += 1
    def occ(c, i):
        if c not in C: return 0
        block = i//B; cnt = checkpoints[c][block]
        for j in range(block*B, i):
            if BWT[j]==c: cnt += 1
        return cnt
    def fm_search(pattern):
        lo, hi = 0, n
        for ch in reversed(pattern):
            if ch not in C: return []
            lo = C[ch]+occ(ch,lo); hi = C[ch]+occ(ch,hi)
            if lo>=hi: return []
        return sorted(SA[lo:hi])
    out = []
    for _ in range(Q):
        input()  # パターン長
        P = input().strip()
        pos = fm_search(P)
        out.append(' '.join(map(str,pos)) if pos else '')
    print('\n'.join(out))

solve()

Step-by-Step 解説

1Suffix Array の構築（SA-IS）
S に番兵 '$' を追加して SA-IS で $O(N)$ 構築。SA[i] = ソート済み i 番目の suffix の開始位置。

2BWT の計算
$BWT[i] = S[SA[i]-1]$。BWT は同じ文字が連続して出現しやすく圧縮に適した変換。

3C 配列と Occ 関数
$C[c]$ = BWT 中で $c$ より辞書順小さい文字の総数。$Occ(c, i)$ = BWT[0..i-1] での $c$ の出現数（チェックポイント + 線形スキャンで $O(B)$）。

4LF-Mapping による検索
パターン $P$ を右から読み、$(lo, hi)$ を絞り込む：
$lo = C[c] + Occ(c, lo)$、$hi = C[c] + Occ(c, hi)$
最終 $[lo, hi)$ の SA 値が出現位置。

5チェックポイントで Occ を高速化
間隔 $B$ ごとに累積カウントを保存。クエリは $O(B)$。$B = O(\sqrt{N})$ で全体 $O(|P| \cdot \sqrt{N})$。

計算量

よくあるミス

ミス	原因	正しい書き方
番兵文字を忘れる	SA の最後 suffix が正しく定まらない	`S + '$'`（'$' < 'a'）
C[c] を S ベースで計算	BWT ベースでないと LF-Mapping が破綻	`for c in BWT: freq[c] += 1`
occ(c, hi) で hi = n を超える	チェックポイント配列の境界外	`chk = (n//B) + 2` で余裕を持つ
reversed(pattern) を忘れる	検索方向が逆になる	`for ch in reversed(pattern):`

次のステップ

発展問題: 動的テキスト（挿入・削除あり）での全文検索（Dynamic BWT / Wavelet Tree BWT）
関連: SA + LCP 配列による最長共通部分文字列の $O(N)$ 計算
応用: 圧縮全文索引（CSA）— BWT + Wavelet Tree で $O(|P| \log \sigma)$ クエリ

自己評価

完全理解おおむね理解要復習