RRR

ウェーブレット行列というデータ構造を使えば、文字列に対するrankクエリ（i.e.文字列の先頭から位置iまでの間にアルファベットcが何回出現したか）をO(log C)で処理できます。ただしCはアルファベットの数です。このテクニックを使うと、長さMのクエリ文字列に対して完全一致検索クエリをO(MlogC)で処理できます。このアルゴリズムをFM-indexと呼びます。FM-indexは時間計算量が小さいだけではなく、そのために事前に用意すべき追加データがO(N)で済む（ただしNはリファレンス文字列の長さ。正確にはO(NlogC)でしょうが、慣習としてここのlogCは省略します）という長所があります。ちなみにこの「検索を高速化するために事前に用意する追加データ」は『索引』という単語の意味そのものです。全文検索用の索引で、かつリファレンス文字列自体よりも小さな索引を『圧縮全文索引』と呼びます。

以前の記事で、FM-indexによる全文検索の具体的なコードを示しましたが、そのコードの索引は圧縮全文索引にはなっていませんでした。まず、示したコードの手順は以下の通りでした。（Compressed Suffix Arrayは作っていませんでしたが、作るとした場合の手順です）
・構築ステップ
(1)リファレンス文字列を使ってSuffix Arrayを構築する。
(2)Suffix Arrayを使ってCompressed Suffix ArrayとBWT文字列を構築する。
(3)BWT文字列を使ってWavelet Matrixを構築する。
・検索ステップ
(4)Wavelet Matrixを使って、完全一致箇所を意味するSuffix Array上の区間を求める。
(5)Compressed Suffix Arrayを読み、その区間内の値を列挙する。それらがリファレンス文字列上の完全一致箇所そのものである。
構築ステップを終えて検索ステップに移るとき、生のSuffix ArrayとBWT文字列は捨ててよいです。Wavelet MatrixとCompressed Suffix Arrayだけ保持する必要があります。Compressed Suffix Arrayは好きなだけ圧縮できますが、問題はWavelet Matrixです。Wavelet MatrixはLogC個のSuccinct Bit Vectorから成ります。そしてあのコードにおいて各Succinct Bit Vectorは長さNのビット列と追加データから成っていました。結局、あのコードにおけるWavelet Matrixは元文字列のコピーと同じ大きさのビット列を含んでおり、ゆえに元文字列よりも大きいものでした。『圧縮全文索引』であるためにはその『長さNのビット列』を何らかの方法で可逆圧縮する必要があります。そのビット列は具体的にはBWT文字列をある意味でソートしてからBit Slicingして得られたビット列でした。一般にビット列の1と0の出現確率が大きく異なれば、少ないほうの出現位置を可変長配列で持つことで自明に可逆圧縮できますが、それは期待できません。しかし幸いなことに、実データに対するBWTの性質からして同じ文字が連続しがちなはずであり、ゆえに『ビット列』では同じビットが連続しがちなはずです。この性質を利用してビット列を可逆圧縮する方法として、RRRというアルゴリズムが岡野原本のp52以降で紹介されています。

RRRは元論文の筆者たちのイニシャルです。具体的には岡野原本を読んでください。

LZ77 + Hybrid Indexing

Hybrid Indexing Revisited (ALENEX2018) で用いられているアルゴリズムについて説明します。

この論文の筆者らは、先行研究においてHybrid Indexingという圧縮全文索引を提案していました。
[1306.4037] Hybrid Indexes for Repetitive Datasets
https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2013.0137
ALENEX2018の論文では、その提案手法が他のいくつかの索引手法と比べて優れていることを様々な実データで定量的に実証しています。実装は以下のリポジトリにあります。
github.com

この記事では準備としてLZ圧縮について述べたあと、Hybrid Indexingそのものについて説明します。

LZ圧縮(or LZ分解)

LZ圧縮は辞書式圧縮という文字列圧縮手法の一つです。（LZ圧縮とLZ分解は異なる概念ですが、以降の文では気にせず全てLZ圧縮と書いています）具体的には、文字列を頭から順番に見ていって、以前に登場した部分文字列を見つけたら、その部分文字列の開始位置と長さで置き換えます。詳しく言うと、ある文字iから始まるsuffixが、j<iなるjから始まるsuffixと長さkのcommon prefixを持つとき、部分文字列[i,i+k)をタプル(j,k)で置き換えます。jが存在しない場合は(文字、0)としたり置き換えなかったりします。jが存在しない場合の具体的な扱いによってLZ77の他にLZSSなどいろいろなバリエーションがありますが、今回は特に気にせずに説明します。以下の説明では文字の添字は0-originとします。

例えば"abcabcd"をLZ圧縮すると、"(a,0),(b,0),(c,0),(0,3),(d,0)"となります。最初の3文字は初出なのでそのまま出力しています。3番目からのsuffixは0番目からのsuffixと長さ3のcommon prefixを持つので、第4項でそれを表しています。6文字目は初出なのでそのまま出力しています。

（論文中の例ですが）例えば"zzzzzapzap"をLZ圧縮すると、"(z,0),(0,4),(a,0),(p,0),(4,3)"となります。この例のポイントは第2項で、置き換えが発生する際に自分自身の一部を参照するような置き換えが許されています。もしそのような置き換えが許されないならば、"(z,0),(0,1),(0,2),(0,1),(a,0),(p,0),(4,3)"などとせざるを得ず、圧縮効率が悪くなります。（ちなみにLZ系の他の亜種ではこれの扱いが色々異なります）とはいえ、Hybrid IndexingにおいてはLZ77圧縮がどんなに非効率な形であっても全文検索の出力は等価で計算時間と索引の大きさが変わるだけです。

Observation

以降、リファレンス文字列Rにおけるクエリ文字列Qの完全一致検索クエリを考えます。複数の完全一致があるとき、一致箇所がオーバーラップしているものも含めて全て報告するとします。（例えばR="zzzzzapzap"、Q="zzz"のとき、[0,3),[1,4),[2,5)の3箇所で一致するとします。）クエリ文字列の長さは一定値とし、Mとします。M>1とします。ここでObservationとして、R上でEから始まる部分文字列[E,E+M)がQと完全一致するとき

[E,E+M)はRをLZ圧縮したときの項と項の境界をまたがない ⇒ Eより前に他の完全一致箇所が存在する

が成り立ちます。これは自明でして、LZ圧縮の2文字以上の項は、それより前の文字列との完全一致箇所だからです。

そして、このObservationの対偶は

[E,E+M)はR上で最初に出現する完全一致である ⇒ [E,E+M)はRをLZ圧縮した時の項と項の境界をまたぐ

となります。

例えばリファレンス文字列"abcabcd"からクエリ文字列"ab"を探したいとします。リファレンス文字列をLZ圧縮すると"(a,0),(b,0),(c,0),(0,3),(d,0)"となります。[3,5)は完全一致であり、LZ圧縮の第4項に完全に含まれますが、実は[3,5)より前にも完全一致箇所[0,2)が存在します。[0,2)は最初の完全一致ですが、LZ圧縮の第1項と第2項にまたがっています。

ちなみに逆は成り立ちません。例えばR="faabcdefcdefabcd"とすると、LZ圧縮は(f,0),(a,0),(1,1),(b,0),(c,0),(d,0),(e,0),(0,1),(4,4),(2,4)となります。ここでQ="fa"とすると完全一致箇所は[0,2),[11,13)の2つですが、両方とも境界をまたいでいます。

Kernel Sequence

LZ圧縮の項のうち、(M*2-1)文字以上の項に関して、両端(M-1)文字だけを残して真ん中はダミー文字で置き換えた文字列をKernel Sequenceと呼びます。

例えばR="abcdefghijabcdefghij"でM=3のとき、ダミー文字を#とするとK="abcdefghijab#ij"となります。リファレンス文字列の中に完全な長い反復が大量に存在するとき、Kernel Sequenceは非常に短くなり得ます。

Kernel Sequence上で既存の（FM-indexなどを使った）完全一致検索を行うと、RとQの完全一致箇所が1つ以上存在するならば、少なくとも最初の完全一致は必ず見つかることが前述のObservationにより言えます。ここで見つかった一致箇所をPrimary Matchと呼ぶことにします。

Secondary Match

少なくとも最初の完全一致は必ず見つかりますが、2つ目以降の完全一致(Secondary Match)は見つからない場合があります。それは完全一致箇所がダミー文字で潰された領域を含む場合です。ここで追加のデータ構造としてRMQを使えばSecondary Matchを高速に全て見つけられます。

(M*2-1)文字以上の各項[i,j)について、開始位置と終端位置の組(i,j)を集めたとします。（ここでi,jと言っているのは参照している先の番号です。例えばR="abcdefghijabcdefghij"でM=3のとき、(i,j)=(0,10)のことです。(10,20)ではありません）ここで高速に実行したいクエリは、既に見つかっている完全一致箇所[s,e)に対して、なる項(i,j)を全て見つけることです。それができれば、項[i,j)には完全一致箇所が含まれることがわかります。見つけた(i,j)で再帰的に同じクエリを投げて深さ優先探索的なことをやれば、全ての完全一致箇所を得ることができます。

方法ですが、まずすべての(i,j)をiで昇順にソートします。iが同じ場合のjの順番は気にしなくても良いです。RMQを使えば、の範囲内におけるjの最大値とそのときの(i,j)のインデックスを高速に得られます。得た(i,j)のjがe以上ならば、インデックスの前後で区切って再帰的にRMQを投げることで、なる項(i,j)を全て見つけられます。

まとめ

Kernel Sequence上のFM-index（あるいは他の既存の全文索引）とRMQ用のデータ構造を組み合わせることで、完全一致検索を行えます。LZ圧縮がうまく働くような文字列に対して、もともとの文字列でFM-indexを作るよりも索引の容量が小さくなる利点があります。また計算時間も（ALENEX2018の論文によると）低く抑えられます。

Wavelet Matrix + BWT + FM-index

以前の記事で何度か触れましたが、FM-indexという方法を使えば完全一致検索を高速に処理できます。他の良い解説の紹介とコード例をあげておきます。

		構築時間	検索時間(平均)	検索時間(最悪)
(1)	Suffix Array	O(N) (SAIS)	O(M+logN)	O(MlogN)
(2)	(1)+LCP Array+RMQ	O(N) (Kasaiの方法)	O(M+logN)	O(M+logN)
(3)	(1)+Compression	O(N)†	O(M+logN)	O(MlogN)	(†時間は定数倍大きく、空間は定数倍小さい)
(4)	(2)+(3)	O(N)†	O(M+logN)	O(M+logN)
(5)	(4)+BWT+Wavelet Matrix+FM-index	O(N)†	O(MlogC)	O(MlogC)

Nはリファレンス文字列の長さ、Mはクエリ文字列の長さ、Cは文字の種類

Wavelet Matrix(ウェーブレット行列)

詳しくは説明しませんが、前回の記事で紹介したSuccinct Bit Vectorを応用して、文字列に対してrank操作を可能にするものです。 ウェーブレット木の効率的で簡単な実装 "The Wavelet Matrix" - EchizenBlog-Zwei によると、元論文は The Wavelet Matrix です。私が知る限り最もわかりやすい解説ブログ記事は　中学生にもわかるウェーブレット行列 - アスペ日記です。岡野原本(高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学))にはp74で説明されていますが、ウェーブレット行列の構築とrank操作の具体的な実装方法は書かれていません。

既存の高速なライブラリとしては GitHub - simongog/sdsl-lite: Succinct Data Structure Library 2.0 GitHub - herumi/cybozulib: a tiny library for C++などがあります。

以下はコード例です。

class BruteForceWaveletMatrix {
private:

	std::vector<uint8_t>v;

public:

	BruteForceWaveletMatrix(const std::vector<uint8_t>& input) {
		v = input;
	}

	int64_t rank(const int64_t index, const uint8_t character) {
		//リファレンス文字列の[0,index)のうち、characterが何回出現したか返す。
		int64_t answer = 0;
		for (int i = 0; i < index; ++i)if (v[i] == character)++answer;
		return answer;
	}
};
class WaveletMatrix {
private:

	std::vector<BitVector>bv;
	std::vector<int64_t>acc;
	int64_t size;

	//x以上の2べき乗数のうち最小の数が2^kとしてkを返す。
	static int64_t log2_ceiling(uint64_t x) {
		if (x <= 0)return 0;
		x--;
		for (int i = 0; i < 6; ++i)x |= x >> (uint64_t(1) << i);
		x = (x & 0x5555555555555555ULL) + ((x & 0xAAAAAAAAAAAAAAAAULL) >> 1);
		x = (x & 0x3333333333333333ULL) + ((x & 0xCCCCCCCCCCCCCCCCULL) >> 2);
		x = (x & 0x0F0F0F0F0F0F0F0FULL) + ((x & 0xF0F0F0F0F0F0F0F0ULL) >> 4);
		x *= 0x0101010101010101ULL;
		return x >> 56;
	}

public:

	WaveletMatrix(const std::vector<uint8_t>& input, const int max_char) {
		//inputは入力文字列で、max_charは入力文字列が含みうる文字の最大値。
		//実際に含んでいる文字の最大値ではない。これを用意する理由は高速化で、
		//具体的にはrank操作においてmax_charを超える引数を受け取らないと仮定する。

		assert(1 <= max_char && max_char <= 255);
		for (const auto x : input)assert(x <= max_char);

		size = input.size();

		//max_charのビット数ぶんBitVectorを用意する。
		const int64_t bits = log2_ceiling(max_char + 1);
		bv = std::vector<BitVector>(bits, BitVector());

		//BitVectorの最上位は、inputの最上位ビットの列を格納する。
		std::vector<uint64_t>bv_tmp((input.size() + 63) / 64, 0);
		const auto SetBit1 = [&](const int i) {bv_tmp[i / 64] |= uint64_t(1) << (i % 64); };
		for (int i = 0; i < input.size(); ++i) {
			if (input[i] & (1 << (bits - 1)))SetBit1(i);
		}
		bv[bits - 1]= BitVector(bv_tmp);

		auto input_tmp = input;
		//BitVectorの最上位以外は、input_tmp配列を(bit_pos+1)番目のビットで安定ソートした後における、
		//bit_pos番目のビットの列を格納する。cf.高速文字列解析の世界p74
		for (int bit_pos = bits - 2; bit_pos >= 0; --bit_pos) {
			int count = 0;
			std::vector<uint8_t>input_tmp2(input.size(), 0);
			for (int i = 0; i < input.size(); ++i) {
				if ((input_tmp[i] & (1 << (bit_pos + 1))) == 0)input_tmp2[count++] = input_tmp[i];
			}
			for (int i = 0; i < input.size(); ++i) {
				if ((input_tmp[i] & (1 << (bit_pos + 1))) != 0)input_tmp2[count++] = input_tmp[i];
			}
			for (int i = 0; i < bv_tmp.size(); ++i)bv_tmp[i] = 0;
			for (int i = 0; i < input.size(); ++i) {
				if (input_tmp2[i] & (1 << bit_pos))SetBit1(i);
			}
			bv[bit_pos] = BitVector(bv_tmp);
			input_tmp = input_tmp2;
		}

		{
			int count = 0;
			std::vector<uint8_t>input_tmp2(input.size(), 0);
			for (int i = 0; i < input.size(); ++i) {
				if ((input_tmp[i] & 1) == 0)input_tmp2[count++] = input_tmp[i];
			}
			for (int i = 0; i < input.size(); ++i) {
				if ((input_tmp[i] & 1) != 0)input_tmp2[count++] = input_tmp[i];
			}
			//acc[i]にはinput_tmp2中でiが最初に出現した位置を格納する。出現しない場合は-1とする。
			acc = std::vector<int64_t>(max_char + 1, -1);
			for (int i = 0; i < input_tmp2.size(); ++i)if (acc[input_tmp2[i]] == -1)acc[input_tmp2[i]] = i;
		}
	}

	int64_t rank(const int64_t index, const uint8_t character) {
		//リファレンス文字列の[0,index)のうち、characterが何回出現したか返す。

		if (acc[character] == -1)return 0;

		//リファレンス文字列の[0,index)のうち、最上位ビットがcharacterに等しいものを数える。
		uint64_t tmp_rank =
			(character & (1 << (bv.size() - 1))) ?
			bv[bv.size() - 1].rank1(index) :
			bv[bv.size() - 1].rank0(index);

		//最上位-1ビットから下位ビットに向けて処理していく。
		for (int bit_pos = bv.size() - 2; bit_pos >= 0; --bit_pos) {
			if (character & (1 << (bit_pos + 1))) {
				const uint64_t start_pos = size - bv[bit_pos + 1].sum1();
				tmp_rank = (character & (1 << bit_pos)) ?
					bv[bit_pos].rank1(start_pos + tmp_rank) :
					bv[bit_pos].rank0(start_pos + tmp_rank);
			}
			else {
				tmp_rank = (character & (1 << bit_pos)) ?
					bv[bit_pos].rank1(tmp_rank) :
					bv[bit_pos].rank0(tmp_rank);
			}
		}
		if (character & 1)tmp_rank += size - bv[0].sum1();
		return tmp_rank - acc[character];
	}
};

BWT

Suffix Arrayのsuffixたちの1個前の文字を集めたものです。

std::vector<uint8_t>BWT(const std::vector<uint8_t>& input) {

	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] != 0);
	assert(input.back() == 0);

	const auto suffix_array = SAISSuffixArrayConstruction(input);

	std::vector<uint8_t>bwt(suffix_array.size(), 0);
	for (int i = 0; i < suffix_array.size(); ++i) {
		if (suffix_array[i] == 0)bwt[i] = input.back();
		else bwt[i] = input[suffix_array[i] - 1];
	}
	return bwt;
}

FM-index

FM-indexの解説は岡野原本の7.5.2節が鉄板なので読みましょう。ライブラリは先ほどと同じくSDSLなどが鉄板として知られています。

以下のコードは簡単のため、queryが来るたびにFM-indexを都度構築するように書かれています。実際には最初に一度だけ構築すれば十分です。

std::vector<int> FMIndexMatch(
	const std::vector<uint8_t>& ref,
	const std::vector<uint8_t>& query) {

	for (const auto x : ref)assert(x != 0);
	for (const auto x : query)assert(x != 0);

	auto ref_ = ref;
	ref_.push_back(0);
	const auto suffix_array = SAISSuffixArrayConstruction(ref_);
	const auto bwt = BWT(ref_);
	//BruteForceWaveletMatrix wavelet_matrix(bwt);
	WaveletMatrix wavelet_matrix(bwt, 255);

	//bwt中でのiより小さい文字の出現回数を求めてacc_num[i]に格納する。
	std::vector<int64_t>num(256, 0), acc_num(256, 0);
	for (const auto x : bwt)++num[x];
	for (int i = 1; i < 256; ++i)acc_num[i] = num[i - 1] + acc_num[i - 1];

	//ここまで構築で、ここから検索

	int64_t start_pos = 0, end_pos = ref_.size();
	for (int i = query.size() - 1; i >= 0; --i) {
		const auto c = query[i];
		start_pos = acc_num[c] + wavelet_matrix.rank(start_pos, c);
		end_pos = acc_num[c] + wavelet_matrix.rank(end_pos, c);
		if (start_pos >= end_pos)return std::vector<int>{};
	}

	std::vector<int>match_start_pos;
	for (int64_t i = start_pos; i < end_pos; ++i)match_start_pos.push_back(suffix_array[i]);
	sort(match_start_pos.begin(), match_start_pos.end());
	return match_start_pos;
}

Vertical Code + Succinct Bit Vector + Compressed Suffix Array

前回までの記事では、完全一致検索を高速化するためにリファレンス文字列に加えてSuffix Arrayを保持していました。しかし、文字列自体よりもSuffix Arrayのほうが大きいのが問題でした。この記事ではSuffix Arrayを圧縮して保持する方法の一つを説明します。文字列そのものよりも小さな追加データ(Compressed Suffix Array, CSA)を用いて、引数iを受け取ってsuffix_array[i]を返す機能を作れます。そのための基礎としてVertical CodeとSuccinct Bit Vectorが必要になりますが、この記事ではそれらについても紹介します。

d.hatena.ne.jp

岡野原大輔, 現実的な圧縮付全文索引

これらの記事を参考にしました。後者の論文の3.1節の最初において、CSAの方法としてGrossiらの論文とSadakaneらの論文がreferされています。後者の論文自身はSadakaneらの論文に基づいています。一方前者のブログ記事はGrossiらの論文に基づいています。今回の記事ではSadakaneらの論文の方法に基づいています。

Vertical Code

現実的な圧縮付全文索引、Vertical Codeを調べたよ - EchizenBlog-Zweiで説明されています。

以下のコード例を示します。popcount関数は引数の立っているビットを数えるもので、SSE4.2以降にpopcntという等価な命令がありますが今回は基本演算で実装しました。トリッキーな実装になっていますが最速ではありません。cf. 2009-07-04 - 当面C#と.NETな記録

ブロックサイズが64で決め打ちされていますが、これは簡単のためであり、ブロックサイズ64で性能が最善になるとは限りません。

class VerticalCode {
private:

	static constexpr int block = 64;
	uint64_t size;
	std::vector<uint64_t>psi;
	std::vector<uint64_t>code;
	std::vector<uint64_t>code_index;

	static uint64_t popcount(uint64_t x) {
		x = (x & 0x5555555555555555ULL) + ((x & 0xAAAAAAAAAAAAAAAAULL) >> 1);
		x = (x & 0x3333333333333333ULL) + ((x & 0xCCCCCCCCCCCCCCCCULL) >> 2);
		x = (x & 0x0F0F0F0F0F0F0F0FULL) + ((x & 0xF0F0F0F0F0F0F0F0ULL) >> 4);
		x *= 0x0101010101010101ULL;
		return x >> 56;
	}

public:

	VerticalCode(const std::vector<uint64_t>& input) {
		//inputが元の数列だとして、差分化してvertical codeで保持する。
		//inputの各要素は[0,N)の並べ替えであることを仮定する。

		{
			auto testarray = input;
			sort(testarray.begin(), testarray.end());
			for (int i = 0; i < testarray.size(); ++i)assert(testarray[i] == i);
		}

		size = input.size();
		const int64_t code_length = (input.size() + block - 1) / block;

		//psi[i]には、0～(i-1)ブロック内の差分列の総和をsizeで割った余りを格納する。
		//「差分列の総和」は元の数列の値とは異なる。なぜなら、
		//ここで言う「差分列」は、圧縮しやすい単調増加列にするため、差分に-1や+Nしているからである。
		psi = std::vector<uint64_t>(code_length, 0);
		//for (int64_t i = 1; i < code_length; ++i)psi[i] = input[i * block - 1];

		code.clear();
		code_index = std::vector<uint64_t>(code_length + 1, 0);

		for (int64_t b = 0; b < code_length; ++b) {
			const int64_t offset = b * block;

			//offset番目からblock個について、単調増加な差分列を作る。
			std::vector<uint64_t>diff(std::min(int64_t(64), int64_t(input.size() - offset)), 0);
			uint64_t diffmax = 0;
			for (int64_t i = 0; i < diff.size(); ++i) {
				const int64_t j = offset + i;
				if (j == 0)diff[i] = input[j];
				else if (input[j - 1] < input[j])diff[i] = input[j] - input[j - 1] - 1;
				else {
					//↓はinputが[0,N)の並べ替えであれば必ず成り立つ。
					//assert(input[j - 1] - input[j] < input.size());
					diff[i] = input[j] - input[j - 1] - 1 + input.size();
				}
				diffmax = std::max(diffmax, diff[i]);
			}
			const int64_t bit_demand = SparseTableMin<uint64_t>::log2_ceiling(diffmax + 1);
			
			code.resize(code.size() + bit_demand);
			code_index[b + 1] = code.size();
			for (int i = 0; i < diff.size(); ++i) {
				for (uint64_t num = diff[i], j = code_index[b]; num; num /= 2, ++j) {
					code[j] |= (num & uint64_t(1)) << i;
				}
			}

			if (b < code_length - 1) {
				psi[b + 1] = psi[b];
				for (const int x : diff)psi[b + 1] += x;
				psi[b + 1] %= size;
			}
		}
	}

	VerticalCode() {
		VerticalCode(std::vector<uint64_t>{0});
	}

	uint64_t Get(const int64_t index) {
		const int64_t b = index / block;
		const auto bits = index - (b * block);

		const uint64_t bitmask = 0xFFFF'FFFF'FFFF'FFFFULL >> (63 - bits);

		uint64_t answer = psi[b];
		for (uint64_t i = code_index[b]; i < code_index[b + 1]; ++i) {
			answer += popcount(code[i] & bitmask) << (i - code_index[b]);
		}
		return (answer + index) % size;
	}
};

Succinct Bit Vector

長さNのビット列に対して、少ない追加データを用意することで、rank操作とselect操作という2種類の操作を定数時間で処理できるようにしたものをSuccinct Bit Vectorと呼びます。曖昧な言い方になっています。実は、長さNのビット列に対してo(N)の追加データによりrankとselectを定数時間で処理できるような方法が存在しますが、以下のコード例はそうなっていません。具体的には、以下のコードでは256bitごとと32bitごとのテーブルが用意されますが、o(N)を実現するためにはこの値を(256と32の決め打ちではなく)Nに応じた適切な値とする必要があります。Bit Vectorを実際に使いたい人は有名なライブラリがあるのでそれを使いましょう。あとCompressed Suffix Arrayではselect操作は構築のときにのみ使われて、クエリ処理のときには使われません。そのため以下のコード例ではselectが単純で低速な実装になっています。

class BruteForceBitVector {
private:

	std::vector<uint64_t>B;

public:

	BruteForceBitVector(const std::vector<uint64_t>& input) {
		B = input;
	}
	BruteForceBitVector() {
		BruteForceBitVector(std::vector<uint64_t>{0});
	}

	bool access(const uint64_t index) {
		return (B[index / 64] & (uint64_t(1) << (index % 64))) != 0;
	}

	uint64_t rank(const uint64_t index) {
		//[0,index)にいくつ立っているビットがあるか求めて返す。
		uint64_t answer = 0;
		for (int i = 0; i < index; ++i)if (access(i))++answer;
		return answer;
	}

	uint64_t select1(const uint64_t count) {
		//count番目(0-origin)の1のビットの位置を求めて返す。低速。
		uint64_t acc = 0;
		for (int i = 0; i < B.size() * 64; ++i)if (access(i))if (acc++ == count)return i;
		return std::numeric_limits<uint64_t>::max();
	}
};
class BitVector {
private:

	std::vector<uint64_t>B;
	std::vector<uint64_t>rank_global_256;
	std::vector<uint8_t>rank_local_64;

public:

	BitVector(const std::vector<uint64_t>& input) {

		B = input;

		//rank_global_256[i]の末尾以外には、[0,256*i)の範囲に立っているビットの数を格納する。
		//末尾にはデータ全体で立っているビットの数を格納する。
		rank_global_256.resize(((input.size() + 3) / 4) + 1);

		//rank_local_64[i]には、[256*(i/4),64*i)の範囲に立っているビットの数を格納する。
		rank_local_64.resize(input.size());

		uint64_t global_sum = 0;
		uint64_t local_sum = 0;
		for (int64_t i = 0; i < input.size(); ++i) {
			if (i % 4 == 0) {
				local_sum = 0;
			}
			rank_local_64[i] = local_sum;
			const uint64_t pop = popcount(input[i]);
			global_sum += pop;
			local_sum += pop;
			if (i % 4 == 3) {
				rank_global_256[i / 4 + 1] = global_sum;
			}
		}
		rank_global_256.back() = global_sum;
	}

	BitVector() {
		BitVector(std::vector<uint64_t>{0});
	}

	BitVector(const BitVector& obj) {
		this->B = obj.B;
		this->rank_global_256 = obj.rank_global_256;
		this->rank_local_64 = obj.rank_local_64;
	}

	bool access(const uint64_t index) {
		return (B[index / 64] & (uint64_t(1) << (index % 64))) != 0;
	}

	uint64_t rank(const uint64_t index) {
		//[0,index)にいくつ立っているビットがあるか求めて返す。
		uint64_t answer = rank_global_256[index / 256] + rank_local_64[index / 64];
		answer += popcount(B[index / 64] & ((uint64_t(1) << (index % 64)) - 1));
		return answer;
	}

	uint64_t select1(const uint64_t count) {
		//count番目(0-origin)の1のビットの位置を求めて返す。低速。

		int64_t L = 0, R = B.size() * 64;
		if (rank_global_256.back() <= count)return std::numeric_limits<uint64_t>::max();
		if (rank_global_256.back() == count - 1 && access(R - 1))return R - 1;

		//欲しい位置は[L,R)のどこかだとする。
		while (L + 1 < R) {
			const int64_t mid = (L + R) / 2;
			if (rank(mid) <= count)L = mid;
			else R = mid;
		}
		return L;
	}

	uint64_t sum1() { return rank_global_256.back(); }

	static uint64_t popcount(uint64_t x) {
		x = (x & 0x5555555555555555ULL) + ((x & 0xAAAAAAAAAAAAAAAAULL) >> 1);
		x = (x & 0x3333333333333333ULL) + ((x & 0xCCCCCCCCCCCCCCCCULL) >> 2);
		x = (x & 0x0F0F0F0F0F0F0F0FULL) + ((x & 0xF0F0F0F0F0F0F0F0ULL) >> 4);
		x *= 0x0101010101010101ULL;
		return x >> 56;
	}
};

Compressed Suffix Array

下記のコードでは、現実的な圧縮付全文索引で説明されているアルゴリズムが実装されています。SAISSuffixArrayConstructionは前回の記事のコードを想定しています。

template<uint64_t compress_order>class CompressedSuffixArray {

	//2^compress_order個ごとにサンプルする。

private:

	VerticalCode psi_vertical_code;
	std::vector<uint64_t>suffix_array_sumple;
	BitVector B;
	int64_t suffix_array_size;

public:

	CompressedSuffixArray(const std::vector<uint8_t>& input) {

		assert(1 <= compress_order && compress_order <= 10);
		for (int64_t i = 0; i < input.size() - 1; ++i)assert(input[i] != 0);
		assert(input.back() == 0);

		const auto suffix_array = SAISSuffixArrayConstruction(input);
		suffix_array_size = suffix_array.size();

		//Inverse Suffix Arrayを作る。
		std::vector<int>inverse_suffix_array(suffix_array.size(), 0);
		for (int64_t i = 0; i < inverse_suffix_array.size(); ++i) {
			inverse_suffix_array[suffix_array[i]] = i;
		}

		//Inverse Suffix Arrayを使ってPsi Arrayを作り、Vertical Codeにする。
		std::vector<uint64_t>psi(suffix_array.size(), 0);
		for (int64_t i = 0; i < psi.size(); ++i) {
			psi[i] =
				(suffix_array[i] == suffix_array.size() - 1) ?
				inverse_suffix_array[0] :
				inverse_suffix_array[suffix_array[i] + 1];
		}
		psi_vertical_code = VerticalCode(psi);

		//ビット列Bを作る。定義は
		//Bのi番目のビットがtrue⇔SA[i]が2^compress_orderで割り切れる
		//とする。(元の岡野原論文の定義)
		//加えてSA[i]==SA.size()-1のときもサンプルを保存しておく。
		//そうすることでGetのときに剰余を取る必要がなくなる。
		const auto interval = uint64_t(1) << compress_order;
		std::vector<uint64_t>b_tmp((suffix_array.size() + 63) / 64);
		for (int i = 0; i < suffix_array.size(); ++i) {
			if ((suffix_array[i] % interval) == 0 || suffix_array[i] == suffix_array.size() - 1) {
				b_tmp[i / 64] |= uint64_t(1) << (i % 64);
			}
		}
		B = BitVector(b_tmp);

		//Bのビットが立っている位置のsuffix arrayだけをサンプリングして保存しておく。
		const int64_t sample_size = B.sum1();
		suffix_array_sumple = std::vector<uint64_t>(sample_size, 0);
		for (int64_t i = 0; i < sample_size; ++i) {
			suffix_array_sumple[i] = suffix_array[B.select1(i)];
		}
	}

	uint64_t Get(int64_t index) {

		//suffix_array[index]を返す。

		for(int64_t move_count = 0;;++move_count) {

			if (B.access(index)) {

				//剰余について
				//この記述の中ではk=2^compress_order, t=move_countとする。
				//元の岡野原論文の定義だと、B.access(index)がtrueになる⇔SA[index]がkで割り切れる　であった。
				//B.access(index)がfalseの場合、index:=psi[index]とすることで、"返り値が1増える"効果があるので、
				//たかだかk-1回繰り返しpsiを作用させることで、B.access(index)がtrueになる。
				//そのとき、psiを作用させた回数をtとすると、SA_sample[B.rank(index)]-tを返せば良いと書かれていた。
				//しかし実際には、SA.size()がkで割り切れないとき、SA_sample[B.rank(index)]==0になることがある。
				//そのため、元論文の通りSA[index]がkで割り切れるもののみサンプリングしたなら、
				//元論文と違いSA.size()で割った余りを返す必要がある。
				//今回の実装では元論文と違いSA[index]==SA.size()-1のときもサンプリングしているので、
				//tが非ゼロのときSA_sample[B.rank(index)]==0にはならず、剰余の処理は必要ない。

				return (int64_t(suffix_array_sumple[B.rank(index)]) - move_count);
				//return (int64_t(suffix_array_sumple[B.rank(index)]) - move_count + suffix_array_size) % suffix_array_size;
			}
			index = psi_vertical_code.Get(index);
		}

		assert(0);
		return 0;
	}
};

Suffix Array + LCP Array + RMQ

前回の記事で、リファレンス文字列に加えてSuffix Arrayが与えられれば完全一致検索クエリを高速に処理できることを紹介しました。そのコードの最悪時間計算量は(ただしNはリファレンス文字列の長さ、Mはクエリ文字列の長さ)でした。実はSuffix Arrayに加えてLongest Common Prefix(LCP) Arrayというデータを追加で持っておくことで、完全一致検索クエリを最悪時間計算量で処理できます。この記事ではLCP Arrayを使って完全一致検索する具体的な方法を紹介します。ただしちなみに、前回の記事の方法でもリファレンス・クエリ文字列がともにランダムならば平均時間計算量はです。実際のデータの性質によっては計算時間も前回の記事の方法のほうが速い場合があります。また、実際にはFM-indexなどのより新しい方法がよく使われます。

LCP arrayとは

短く言うと、Suffix Arrayで隣接するSuffix間の共通するPrefixの長さです。詳しい説明はネット上にいっぱいあるので省きます。"banana$"を例に取ると、下記の通り{0,1,3,0,0,2}となります。

開始位置	suffix	次のsuffixと共通するprefixの最大長
6	$	0
5	a$	1
3	ana$	3
1	anana$	0
0	banana$	0
4	na$	2
2	nana$

LCP Arrayは愚直に求めようとすると最悪時間計算量かかります。具体的には例えば同じ文字がひたすら並んでいるケースでそうなります。Kasaiの方法を使うと最悪時間計算量で構築できます。詳しい説明は省きますがコード例は下記の通りです。

std::vector<int>BruteForceLCPArrayConstruction(
	const std::vector<uint8_t>& input,
	const std::vector<int>& suffix_array) {

	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] != 0);
	assert(input.back() == 0);
	assert(input.size() == suffix_array.size());

	std::vector<int>lcp_array(input.size() - 1, 0);
	for (int i = 0; i < input.size() - 1; ++i) {
		while (
			input[suffix_array[i] + lcp_array[i]] ==
			input[suffix_array[i + 1] + lcp_array[i]]) {
			lcp_array[i]++;
		}
	}

	return lcp_array;
}

std::vector<int>KasaiLCPArrayConstruction(
	const std::vector<uint8_t>& input,
	const std::vector<int>& suffix_array) {

	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] != 0);
	assert(input.back() == 0);
	assert(input.size() == suffix_array.size());

	std::vector<int>inverse_suffix_array(suffix_array.size(), 0);
	for (int i = 0; i < inverse_suffix_array.size(); ++i) {
		inverse_suffix_array[suffix_array[i]] = i;
	}

	std::vector<int>lcp_array(input.size() - 1, 0);
	int lcp = 0;
	for (int i = 0; i < input.size() - 1; ++i) {
		if (inverse_suffix_array[i] == input.size() - 1)continue;
		const int pos1 = suffix_array[inverse_suffix_array[i]];
		const int pos2 = suffix_array[inverse_suffix_array[i] + 1];
		while (input[pos1 + lcp] == input[pos2 + lcp])lcp++;
		lcp_array[inverse_suffix_array[i]] = lcp;
		if (lcp > 0)--lcp;
	}
	return lcp_array;
}

RMQについて

LCPを利用して最悪時間計算量で完全一致検索を行うためには、任意のSuffix間のLCPをで求める必要があります。ここで任意の(i,j)について、辞書順でi番目とj番目のSuffix間のLCPはであることを利用します。するともし数列の任意の区間内の最小値をで求められるのならば、LCP Arrayの区間[i,j)内の最小値をその方法で求めれば、その値が辞書順でi番目とj番目のSuffix間のLCPとなります。実は、数列の任意の区間内の最小値を求めるクエリはRange Minimum Query(RMQ)という名前で知られており、で処理できるデータ構造が存在します。

	構築	更新	RMQ
ナイーブ	O(N)	O(1)	O(N)
Segment Tree	O(N)	O(logN)	O(logN)
Sparse Table	O(NlogN)	O(N)	O(1)

Segment TreeとSparse Tableは蟻本に載っていて、ネット上にも色々情報があります。ちなみに、で構築できるSparse Tableもあります。
qiita.com

コード例は以下の通りです。

template<typename T>class BruteForceMin {
private:
	std::vector<T>v;

public:

	BruteForceMin(const std::vector<T>& input) {
		v = input;
	}
	BruteForceMin(const int64_t N) {
		v = std::vector<T>(N, std::numeric_limits<T>::max());
	}
	BruteForceMin() {
		BruteForceMin(1);
	}

	//v[index]=numberとする。
	void update(const int64_t index, const T number) {
		v[index] = number;
	}

	//[index]の値を返す。
	T getnum(const int64_t index) {
		return v[index];
	}

	//[L,R)の最小値を求める。
	int getmin(const int64_t L, const int64_t R) {
		T answer = std::numeric_limits<T>::max();
		for (int64_t i = L; i < R; ++i)if (v[i] < answer)answer = v[i];
		return answer;
	}
};
template<typename T>class SegTreeMin {
private:

	int64_t SIZE;
	std::vector<T>v;

	T getmin_inner(
		const int64_t queryL, const int64_t queryR, const int64_t index,
		const int64_t segL, const int64_t segR) {
		if (queryR <= segL || segR <= queryL)return std::numeric_limits<T>::max();
		if (queryL <= segL && segR <= queryR)return v[index];
		return std::min(
			getmin_inner(queryL, queryR, index * 2,
				segL, (segL + segR) / 2),
			getmin_inner(queryL, queryR, index * 2 + 1,
				(segL + segR) / 2, segR));
	}

public:

	SegTreeMin(const std::vector<T>& input) {
		SIZE = roundup_pow2(input.size());
		v = std::vector<T>(SIZE * 2, std::numeric_limits<T>::max());
		for (int64_t i = 0; i < input.size(); ++i)v[i + SIZE] = input[i];
		for (int64_t i = SIZE - 1; i; --i)v[i] = std::min(v[i * 2], v[i * 2 + 1]);
	}

	SegTreeMin(const int64_t N) {
		SegTreeMin(std::vector<T>(N, std::numeric_limits<T>::max()));
	}

	SegTreeMin() {
		SegTreeMin(1);
	}

	//v[index]=numberとする。
	void update(const int64_t index, const T number) {
		v[index + SIZE] = number;
		for (int64_t i = (index + SIZE) / 2; i; i = i / 2) {
			v[i] = std::min(v[i * 2], v[i * 2 + 1]);
		}
	}

	//[index]の値を返す。
	T getnum(const int64_t index) {
		return v[index + SIZE];
	}

	//[L,R)の最小値を返す。
	T getmin(const int64_t L, const int64_t R) {
		return getmin_inner(L, R, 1, 0, SIZE);
	}

	//x以上であるような2のべき乗数のうち最小のものを返す。
	static uint64_t roundup_pow2(uint64_t x) {
		x--;
		for(int i = 0; i < 6; ++i)x |= x >> (1 << i);
		return x + 1;
	}
};
template<typename T>class SparseTableMin{
private:

	int64_t rank;
	std::vector<std::vector<T>>v;

public:

	SparseTableMin(const std::vector<T>& input) {
		rank = std::max(int64_t(1), log2_ceiling(input.size()));
		v = std::vector<std::vector<T>>(rank,
			std::vector<T>(input.size(), std::numeric_limits<T>::max()));
		for (int64_t i = 0; i < input.size(); ++i)v[0][i] = input[i];
		for (int64_t x = 1; x < rank; ++x) {
			for (int64_t i = 0; i < input.size(); ++i) {
				const int64_t pos = i + (1 << (x - 1));
				if (pos >= input.size())v[x][i] = v[x - 1][i];
				else v[x][i] = std::min(v[x - 1][i], v[x - 1][pos]);
			}
		}
	}

	SparseTableMin(const int64_t N) {
		SparseTableMin(std::vector<T>(N, std::numeric_limits<T>::max()));
	}
	SparseTableMin() {
		SparseTableMin(1);
	}

	//v[index]=numberとする。O(N)
	void update(const int64_t index, const T number) {
		v[0][index] = number;
		for (int64_t x = 1; x < rank; ++x) {
			for (int64_t i = std::max(int64_t(0), index - (1 << x) + 1);
				i <= index; ++i) {
				const int64_t pos = i + (1 << (x - 1));
				if (pos >= v[0].size())v[x][i] = v[x - 1][i];
				else v[x][i] = std::min(v[x - 1][i], v[x - 1][pos]);
			}
		}
	}

	T getnum(const int64_t index) {
		return v[0][index];
	}

	//[L,R)の最小値を返す。log2_ceilingがO(1)だとみなすとO(1)
	T getmin(const int64_t L, const int64_t R) {
		if (L + 1 == R)return v[0][L];
		if (R <= L)return std::numeric_limits<T>::max();
		const int64_t query_length = R - L;
		const int64_t index_rank = log2_ceiling(query_length) - 1;
 		return std::min(v[index_rank][L], v[index_rank][R - (1 << index_rank)]);
	}

	//x以上の2べき乗数のうち最小の数が2^kとしてkを返す。
	static int64_t log2_ceiling(uint64_t x) {

		if (x <= 0)return 0;
		x--;

		//msbの位置を求めるために、ビットを下位に伝播させてpopcntする。

		for (int i = 0; i < 6; ++i)x |= x >> (uint64_t(1) << i);
		x = (x & 0x5555555555555555ULL) + ((x & 0xAAAAAAAAAAAAAAAAULL) >> 1);
		x = (x & 0x3333333333333333ULL) + ((x & 0xCCCCCCCCCCCCCCCCULL) >> 2);
		x = (x & 0x0F0F0F0F0F0F0F0FULL) + ((x & 0xF0F0F0F0F0F0F0F0ULL) >> 4);
		x *= 0x0101010101010101ULL;
		return x >> 56;
	}
};
template<typename T>class AlstrupSparceTableMin {
private:

	int64_t block_size;
	SparseTableMin<T> macro_array;
	std::vector<T>v;
	std::vector<std::vector<uint64_t>>blocks;

	T getmin_micro(const int64_t index, const int64_t L, const int64_t R) {
		const uint64_t w = blocks[index][R] & ~((uint64_t(1) << L) - uint64_t(1));
		const int64_t lsb =
			SparseTableMin<T>::log2_ceiling(w & uint64_t(-int64_t(w)));
		return v[index * block_size + ((w == 0) ? R : lsb)];
	}

public:

	AlstrupSparceTableMin(const std::vector<T>& input) {
		v = input;
		block_size = std::max(int64_t(1),
			SparseTableMin<T>::log2_ceiling(input.size()) / 2);
		assert(block_size <= 63);
		const int64_t block_num = (input.size() + block_size - 1) / block_size;
		{
			std::vector<int64_t>M(block_num, std::numeric_limits<T>::max());
			for (int64_t i = 0; i < input.size(); ++i) {
				const int64_t index = i / block_size;
				M[index] = std::min(M[index], int64_t(input[i]));
			}
			macro_array = SparseTableMin<T>(M);
		}
		blocks = std::vector<std::vector<uint64_t>>(block_num,
			std::vector<uint64_t>(block_size, 0));
		for (int64_t x = 0; x < blocks.size(); ++x) {
			std::vector<T>B(block_size, std::numeric_limits<T>::max());
			const int64_t offset = x * block_size;
			for (int64_t i = 0;
				offset + i < input.size() && i < block_size; ++i) {
				B[i] = input[offset + i];
			}
			std::vector<int64_t>g(block_size, -1);
			std::stack<int64_t>g_stack;
			for (int64_t i = 0; i < block_size; ++i) {
				while (!g_stack.empty() && B[i] <= B[g_stack.top()]) {
					g_stack.pop();
				}
				g[i] = g_stack.empty() ? -1 : g_stack.top();
				g_stack.push(i);
			}
			for (int64_t i = 1; i < block_size; ++i) {
				blocks[x][i] = g[i] == -1 ? 0 :
					(blocks[x][g[i]] |
					(uint64_t(1) << (uint64_t(g[i]))));
			}
		}
	}

	AlstrupSparceTableMin(const int64_t N) {
		AlstrupSparceTableMin(std::vector<T>(N, std::numeric_limits<T>::max()));
	}

	AlstrupSparceTableMin() {
		AlstrupSparceTableMin(1);
	}

	//v[index]=numberとする。O(N)
	void update(const int64_t index, const T number) {

		v[index] = number;

		const int64_t block_index = index / block_size;
		int64_t macro_min = std::numeric_limits<T>::max();
		const int64_t offset = block_index * block_size;
		for (int64_t i = 0; i < block_size && offset + i < v.size(); ++i) {
			macro_min = std::min(macro_min, v[offset + i]);
		}
		macro_array.update(block_index, macro_min);

		std::vector<T>B(block_size, std::numeric_limits<T>::max());
		for (int64_t i = 0; offset + i < v.size() && i < block_size; ++i) {
			B[i] = v[offset + i];
		}
		std::vector<int64_t>g(block_size, -1);
		std::stack<int64_t>g_stack;
		for (int64_t i = 0; i < block_size; ++i) {
			while (!g_stack.empty() && B[i] <= B[g_stack.top()])g_stack.pop();
			g[i] = g_stack.empty() ? -1 : g_stack.top();
			g_stack.push(i);
		}

		for (int64_t i = 1; i < block_size; ++i) {
			blocks[block_index][i] = g[i] == -1 ? 0 :
				(blocks[block_index][g[i]] |
				(uint64_t(1) << (uint64_t(g[i]))));
		}
	}

	T getnum(const int64_t index) {
		return v[index];
	}

	//[L,R)の最小値を返す。log2_ceilingがO(1)だとみなすとO(1)
	int getmin(const int64_t L, const int64_t R) {
		if (R <= L)return std::numeric_limits<T>::max();
		const int64_t L_index = L / block_size;
		const int64_t R_index = (R - 1) / block_size;
		if (L_index == R_index) {
			return getmin_micro(L_index, L % block_size, (R - 1) % block_size);
		}
		T answer = macro_array.getmin(L_index + 1, R_index);
		answer = std::min(
			answer, getmin_micro(L_index, L % block_size, block_size - 1));
		answer = std::min(
			answer, getmin_micro(R_index, 0, (R - 1) % block_size));
		return answer;
	}
};

msbとpopcountはAMD64に専用命令がありますが、今回は簡単のため専用命令を使わずに実装しました。（log2_ceiling）

LCP Arrayを利用した完全一致検索

完全一致検索を実際に行うコード例はネット上にありそうであまり無いので以下に記しておきます。以下のコードでは簡単のためSuffix ArrayとLCP Arrayを関数内で都度構築するように書かれていますが、実際に運用したい人はそれらを引数として受け取るように適宜書き換えて下さい。とはいえ完全一致検索を実際に必要とする人はFM-indexなどを使うほうが良い場合が多いと思います。

std::vector<int> SuffixArrayLCPMatch(
	const std::vector<uint8_t>& ref,
	const std::vector<uint8_t>& query) {

	for (const auto x : ref)assert(x != 0);
	for (const auto x : query)assert(x != 0);

	auto ref_ = ref;
	ref_.push_back(0);
	const auto suffix_array = SAISSuffixArrayConstruction(ref_);
	SparseTableMin<int>lcp_array(KasaiLCPArrayConstruction(ref_, suffix_array));

	const auto compare = [&](const int ref_index, const int lcp_lower_bound) {
		//queryと、ref_[k]以降を比較する。
		//ref_[k]以降がquery全体をprefixとして含む場合、queryのほうが若いと判定する。
		//queryには終端文字が付いていないが、あたかも終端文字よりも若い"超終端文字"が付いているかのように判定すると言ってもよい。

		//比較の処理において、先頭の lcp_lower_bound 文字は一致すると仮定する。
		//返り値は、実際に一致した文字数とどちらが先か。

		int i = lcp_lower_bound;

		//ref_.back()は終端文字であり、queryに終端文字は無いので、2つ目の条件は不要である。
		while (i < query.size() &&
			/*ref_index + i < ref_.size() &&*/
			query[i] == ref_[ref_index + i])++i;
		const bool query_is_front = (i == query.size()) || query[i] < ref_[ref_index + i];

		return std::make_pair(i, query_is_front);
	};

	//suffix arrayの最後よりもクエリ文字列のほうが辞書順で後なら、完全一致は存在しない。
	//終端文字はref_にのみ存在し、他のどの文字よりも若いのでそれが言える。
	if (compare(suffix_array.back(), 0).second == false) {
		return std::vector<int>{};
	}

	//suffix array[0]は終端文字列のみの文字列であるから、クエリ文字列よりも辞書順で先である。
	//ゆえに、この時点で、クエリ文字列は辞書順でsuffix arrayのどこかの間に挿入されるはずである。
	//その挿入位置を二分探索で特定する。
	int L = 0, R = suffix_array.size() - 1;
	//これ以降、Lよりも後ろかつRよりも前に挿入されることが判明している。

	//まずはじめに、クエリ文字列と真ん中のsuffixを普通に辞書順比較して、LCPの値と前後関係を求める。
	const int first_mid = (L + R) / 2;
	const auto first_comparison = compare(suffix_array[first_mid], 0);
	int query_ref_lcp = first_comparison.first;
	bool lcp_is_R = first_comparison.second;

	//クエリ文字列が先なら挿入位置は[L,mid]のどこかに絞られる。さもなくば[mid,R]に絞られる。
	if (first_comparison.second)R = first_mid;
	else L = first_mid;

	//以降、「クエリ文字列と片端のsuffixとのLCP」が求まっている状態を保ちつつ二分探索する。
	while (L + 1 < R) {
		const int mid = (L + R) / 2;
		if (lcp_is_R) {
			//この時点で、クエリ文字列とR番目のsuffixとのLCPが求まっていて、query_ref_lcpに格納されている。
			
			//mid番目のsuffixとR番目のsuffixとのLCPを求める。
			//このRMQの実装では半開区間[L,R)の最小値を返す。lcp_array[i]はiとi+1とのLCPなので、第二引数はR+1ではなくRでよい。
			const int mid_R_lcp = lcp_array.getmin(mid, R);

			if (mid_R_lcp < query_ref_lcp) {
				//クエリ文字列とR番目のsuffixが長く一致していて、mid番目のsuffixとR番目のsuffixとが短くしか一致しないのなら、
				//挿入位置は[mid,R]に絞られる。
				//また、[L,mid]内の全suffixに対してクエリは絶対に完全一致しない。
				//なぜなら、LCPがRMQで求まる性質より、任意の正数dに対してLCP(mid,R)>=LCP(mid-d,R)だからである。
				L = mid;
			}
			else if (mid_R_lcp > query_ref_lcp) {
				//クエリ文字列とR番目のsuffixとが短くしか一致せず、mid番目のsuffixとR番目のsuffixとが長く一致するのなら、
				//挿入位置は[L,mid]に絞られる。
				//また、[mid,R]内の全suffixに対してクエリは絶対に完全一致しない。
				//なぜなら、[mid,R]内の全suffixは、query_ref_lcp番目で不一致するという点で共通するからである。
				R = mid;
			}
			else {
				//クエリ文字列とR番目のsuffixとの一致長と、mid番目のsuffixとR番目のsuffixとの一致長が同じなら、
				//クエリ文字列とmid番目のsuffixとを普通に辞書順比較して、LCPの値と前後関係を求める。
				//このとき求まるLCPの値は、R番目のsuffixとのLCPの値以上であることが確定している。
				//ゆえに文字比較を省略でき、最悪時間計算量が削減される。
				const auto comparison = compare(suffix_array[mid], query_ref_lcp);
				query_ref_lcp = comparison.first;
				lcp_is_R = comparison.second;
				if (comparison.second)R = mid;
				else L = mid;
			}
		}
		else {
			//この時点で、クエリ文字列とL番目のsuffixとのLCPが求まっていて、query_ref_lcpに格納されている。
			
			//L番目のsuffixとmid番目のsuffixとのLCPを求める。
			//このRMQの実装では半開区間[L,R)の最小値を返す。lcp_array[i]はiとi+1とのLCPなので、第二引数はmid+1ではなくmidでよい。
			const int L_mid_lcp = lcp_array.getmin(L, mid);

			if (L_mid_lcp < query_ref_lcp) {
				//クエリ文字列とL番目のsuffixとが長く一致していて、L番目のsuffixとmid番目のsuffixとが短くしか一致しないのなら、
				//挿入位置は[L,mid]に絞られる。
				//また、[mid,R]内の全suffixに対してクエリは絶対に完全一致しない。
				//なぜなら、LCPがRMQで求まる性質より、任意の正数dに対してLCP(L,mid)>=LCP(L,mid+d)だからである。
				R = mid;
			}
			else if (L_mid_lcp > query_ref_lcp) {
				//クエリ文字列とL番目のsuffixとが短くしか一致せず、L番目のsuffixとmid番目のsuffixとが長く一致するのなら、
				//挿入位置は[mid,R]に絞られる。
				//また、[L,mid]内の全suffixに対してクエリは絶対に完全一致しない。
				//なぜなら、[L,mid]内の全suffixは、query_ref_lcp番目で不一致するという点で共通するからである。
				L = mid;
			}
			else {
				//クエリ文字列とL番目のsuffixとの一致長と、L番目のsuffixとmid番目のsuffixとの一致長が同じなら、
				//クエリ文字列とmid番目のsuffixとを普通に辞書順比較して、LCPの値と前後関係を求める。
				//このとき求まるLCPの値は、L番目のsuffixとのLCPの値以上であることが確定している。
				//ゆえに文字比較を省略でき、最悪時間計算量が削減される。
				const auto comparison = compare(suffix_array[mid], query_ref_lcp);
				query_ref_lcp = comparison.first;
				lcp_is_R = comparison.second;
				if (comparison.second)R = mid;
				else L = mid;
			}
		}
	}

	//この時点でL+1==Rであり、クエリ文字列は辞書順でR番目のsuffixの直前に位置する。
	//終端文字はref_にのみ存在し、他のどの文字よりも若いので、完全一致するとすればR番目以降である。
	const auto final_comparison = compare(suffix_array[R], 0);

	//クエリ文字列がR番目のsuffixと完全一致しないなら、他のどのsuffixともしないので終了。
	if(final_comparison.first < query.size())return std::vector<int>{};

	std::vector<int>match_start_pos{suffix_array[R]};

	//pos番目のsuffixと完全一致するとき、
	//pos+1番目のsuffixとも完全一致する⇔pos番目のsuffixとpos+1番目のsuffixとのLCPがクエリ文字列長以上である。
	for (int pos = R; pos < suffix_array.size() - 1; pos++) {
		if (lcp_array.getnum(pos) < query.size())break;
		match_start_pos.push_back(suffix_array[pos + 1]);
	}

	sort(match_start_pos.begin(), match_start_pos.end());
	return match_start_pos;
}

suffix array

通常の接尾辞配列(suffix array, SA)について書いておきます。

suffix array とは

短く言うと、全てのsuffixを列挙し、それらを辞書順にソートし、suffixの開始位置の順番を配列に格納したものです。例えば文字列"banana"に対するsuffix arrayは以下の通り{5,3,1,0,4,2}になります。

開始位置	suffix
5	a
3	ana
1	anana
0	banana
4	na
2	nana

suffix arrayが何かについてのより詳しい説明は省きます。

一部の界隈においては、suffix arrayを作る前に入力文字列の末端に終端文字を必ず付ける習慣があります。終端文字は'$'などで表され、他のあらゆる文字より辞書順で若いとされます。その界隈の人が書いたっぽい解説記事では、あらゆる説明の際に文字列の末端に$が付けられています。実際にはべつに$を付けなくても上述のようにsuffix arrayは作れますし使えます。とはいえ実際のところ$を付けることにしたほうが実装がシンプルになるので、この記事においても$を付けることにします。例えば文字列"banana$"に対するsuffix arrayは以下の通り{6,5,3,1,0,4,2}になります。

開始位置	suffix
6	$
5	a$
3	ana$
1	anana$
0	banana$
4	na$
2	nana$

実装がシンプルになることの例

例えばC++で文字列をvector<uint8_t>で持つとします。すると例えば"ana"は{97,110,97}で"anana"は{97,110,97,110,97}になります。"ana"は"anana"の完全なprefixなので、この2つを辞書順で比較するためにはout of rangeを防ぐ処理が必要です。終端文字があればそれが不要になります。

#include <iostream>
#include <vector>
// 終端文字なし
int main(void) {
	const std::vector<uint8_t>ana{ 97,110,97 };
	const std::vector<uint8_t>anana{ 97,110,97,110,97 };
	for (int i = 0;; ++i) {
		if (i == ana.size()) { std::cout << "ana < anana" << std::endl; break; }
		if (i == anana.size()) { std::cout << "ana > anana" << std::endl; break; }
		if (ana[i] < anana[i]) { std::cout << "ana < anana" << std::endl; break; }
		if (ana[i] > anana[i]) { std::cout << "ana > anana" << std::endl; break; }
	}
	return 0;
}

#include <iostream>
#include <vector>
// 終端文字あり
int main(void) {
	const std::vector<uint8_t>ana{ 97,110,97,0 };
	const std::vector<uint8_t>anana{ 97,110,97,110,97,0 };
	for (int i = 0;; ++i) {
		if (ana[i] < anana[i]) { std::cout << "ana < anana" << std::endl; break; }
		if (ana[i] > anana[i]) { std::cout << "ana > anana" << std::endl; break; }
	}
	return 0;
}

ナイーブな構築方法

ナイーブな構築方法では、suffixを実際に全列挙してソートします。時間計算量はです。

std::vector<int> BruteForceSuffixArrayConstruction(const std::vector<uint8_t>& input) {

	//終端文字は0で、終端文字以外は正の数だと仮定する。
	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] > 0);
	assert(input.back() == 0);

	//suffixを全列挙してソートする。
	std::vector<std::vector<uint8_t>>suffix;
	for (int i = 0; i < input.size(); ++i) {
		std::vector<uint8_t> s;
		for (int j = i; j < input.size(); ++j)s.push_back(input[j]);
		suffix.push_back(s);
	}
	sort(suffix.begin(), suffix.end());

	//ソートしたときの順番を記録して返す。ここでは順番をsuffixの長さから求めている。
	std::vector<int>suffix_array;
	for (int i = 0; i < suffix.size(); ++i) {
		suffix_array.push_back(input.size() - suffix[i].size());
	}
	return suffix_array;
}

Larsson-Sadakane法による構築

アルゴリズムの説明は省きます。時間計算量はです。http://mlab.cb.k.u-tokyo.ac.jp/~moris/lecture/2009-winter-experiments/schedule.htmの Larsson_Sadakane.cpp のリンク先に完動するコードが載っています。このコードはLarge-Scale Genome Sequence Processingという本に書かれているコードが元になっていると思われます。

Large-Scale Genome Sequence Processing

SA-IS法による構築

これもアルゴリズムの説明は省きます。時間計算量はです。元論文のAppendix1に完動するコードが載っています。
github.com
下記はこのAppendix1と同様なコードです。GitHubリポジトリのライセンスと同じくCC BY-SA 3.0です。

#include<iostream>
#include<vector>
#include<algorithm>
#include<cassert>

bool TGet(const std::vector<uint8_t>&t, const int index) {
	return (t[index / 8] & (uint8_t(1) << (index % 8))) ? true : false;
}
void TSet(std::vector<uint8_t>&t, const int index, const bool value) {
	t[index / 8] = value ?
		(t[index / 8] | (uint8_t(1) << (index % 8))) :
		(t[index / 8] & (~(uint8_t(1) << (index % 8))));
}
bool IsLMS(const std::vector<uint8_t>&t, const int index) {
	return (index > 0) && TGet(t, index) && (!TGet(t, index - 1));
}
template<typename T>void GetBuckets(
	const T* input,
	std::vector<int>&bucket,
	const int n,
	const int K,
	const bool end) {
	for (int i = 0; i <= K; ++i)bucket[i] = 0;
	for (int i = 0; i < n; ++i)++bucket[input[i]];
	int sum = 0;
	for (int i = 0; i <= K; ++i) {
		sum += bucket[i];
		bucket[i] = end ? sum : (sum - bucket[i]);
	}
}
template<typename T>void InducedSortArrayL(
	const T* input,
	const std::vector<uint8_t>&t,
	int* suffix_array,
	std::vector<int>&bucket,
	const int n,
	const int K,
	const bool end) {
	GetBuckets<T>(input, bucket, n, K, end);
	for (int i = 0; i < n; ++i) {
		const int j = suffix_array[i] - 1;
		if ((j >= 0) && (!TGet(t, j)))suffix_array[bucket[input[j]]++] = j;
	}
}
template<typename T>void InducedSortArrayS(
	const T* input,
	const std::vector<uint8_t>&t,
	int* suffix_array,
	std::vector<int>&bucket,
	const int n,
	const int K,
	const bool end) {
	GetBuckets<T>(input, bucket, n, K, end);
	for (int i = n - 1; i >= 0; --i) {
		const int j = suffix_array[i] - 1;
		if ((j >= 0) && TGet(t, j))suffix_array[--bucket[input[j]]] = j;
	}
}
template<typename T>void SA_IS(const T* input, int* suffix_array, const int n, const int K) {

	std::vector<uint8_t>t((n + 7) / 8, 0);
	TSet(t, n - 2, false);
	TSet(t, n - 1, true);
	for (int i = n - 3; i >= 0; --i) {
		bool sets = (input[i] < input[i + 1]) ||
			(input[i] == input[i + 1] && TGet(t, i + 1));
		TSet(t, i,
			(input[i] < input[i + 1]) ||
			(input[i] == input[i + 1] && TGet(t, i + 1))
		);
	}

	//stage 1
	std::vector<int>bucket(K + 1, 0);
	GetBuckets<T>(input, bucket, n, K, true);
	for (int i = 0; i < n; ++i)suffix_array[i] = -1;
	for (int i = 1; i < n; ++i)if (IsLMS(t, i))suffix_array[--bucket[input[i]]] = i;
	InducedSortArrayL<T>(input, t, suffix_array, bucket, n, K, false);
	InducedSortArrayS<T>(input, t, suffix_array, bucket, n, K, true);
	int n1 = 0;
	for (int i = 0; i < n; ++i) {
		if (IsLMS(t, suffix_array[i]))suffix_array[n1++] = suffix_array[i];
	}
	for (int i = n1; i < n; ++i)suffix_array[i] = -1;
	int name = 0, prev = -1;
	for (int i = 0; i < n1; ++i) {
		int pos = suffix_array[i];
		bool diff = false;
		for (int d = 0; d < n; ++d) {
			if ((prev == -1) ||
				(input[pos + d] != input[prev + d]) ||
				(TGet(t, pos + d) != TGet(t, prev + d))) {
				diff = true;
				break;
			}
			else if (d > 0 && (IsLMS(t, pos + d) || IsLMS(t, prev + d)))break;
		}
		if (diff) {
			name++;
			prev = pos;
		}
		pos = ((pos % 2) == 0) ? (pos / 2) : ((pos - 1) / 2);
		suffix_array[n1 + pos] = name - 1;
	}
	for (int i = n - 1, j = n - 1; i >= n1; --i) {
		if (suffix_array[i] >= 0)suffix_array[j--] = suffix_array[i];
	}

	//stage 2
	int* suffix_array1 = suffix_array;
	int* s1 = suffix_array + n - n1;
	if (name < n1) {
		SA_IS<int>(s1, suffix_array1, n1, name - 1);
	}
	else {
		for (int i = 0; i < n1; ++i)suffix_array1[s1[i]] = i;
	}

	//stage 3
	GetBuckets<T>(input, bucket, n, K, true);
	for (int i = 1, j = 0; i < n; ++i) {
		if (IsLMS(t, i))s1[j++] = i;
	}
	for (int i = 0; i < n1; ++i)suffix_array1[i] = s1[suffix_array1[i]];
	for (int i = n1; i < n; ++i)suffix_array[i] = -1;
	for (int i = n1 - 1; i >= 0; --i) {
		const int j = suffix_array[i];
		suffix_array[i] = -1;
		suffix_array[--bucket[input[j]]] = j;
	}
	InducedSortArrayL<T>(input, t, suffix_array, bucket, n, K, false);
	InducedSortArrayS<T>(input, t, suffix_array, bucket, n, K, true);
}
std::vector<int> SAISSuffixArrayConstruction(const std::vector<uint8_t>& input) {

	for (int i = 0; i < input.size() - 1; ++i)assert(input[i] != 0);
	assert(input.back() == 0);

	std::vector<int>suffix_array(input.size(), 0);
	SA_IS<uint8_t>(&input[0], &suffix_array[0], input.size(), 256);
	return suffix_array;
}
int main(void) {
	std::vector<uint8_t>banana{ 98,97,110,97,110,97,0 };
	const auto suffix_array = SAISSuffixArrayConstruction(banana);
	for (const int x : suffix_array) std::cout << x << std::endl;
	return 0;
}

使いみち

リファレンス文字列のSuffix Arrayがあれば、クエリ文字列との完全一致箇所をsuffix array上の二分探索で求められます。例えば以下のように書けます。以下のコードの返り値は一致箇所の開始点を全列挙したものです。これもCC BY-SA 3.0としておきます。

std::vector<int> SuffixArrayMatch(
	const std::vector<uint8_t>& ref,
	const std::vector<uint8_t>& query) {

	for (const auto x : ref)assert(x != 0);
	for (const auto x : query)assert(x != 0);

	auto ref_ = ref;
	ref_.push_back(0);
	const auto suffix_array = SAISSuffixArrayConstruction(ref_);

	int L = 0, R = suffix_array.size() - 1;
	while (L + 1 < R) {
		const int mid = (L + R) / 2;
		bool query_is_front = true;
		for (int i = 0; i < query.size(); ++i) {
			if (query[i] < ref_[suffix_array[mid] + i]) {
				query_is_front = true;
				break;
			}
			else if (query[i] > ref_[suffix_array[mid] + i]) {
				query_is_front = false;
				break;
			}
		}
		if (query_is_front)R = mid;
		else L = mid;
	}

	std::vector<int>match_start_pos;
	for (int pos = R; pos < suffix_array.size(); pos++) {
		if (suffix_array[pos] > ref_.size() - query.size())break;
		for (int i = 0; i < query.size(); ++i) {
			if (query[i] != ref_[suffix_array[pos] + i]) {
				goto LABEL;
			}
		}
		match_start_pos.push_back(suffix_array[pos]);
	}

LABEL:
	sort(match_start_pos.begin(), match_start_pos.end());
	return match_start_pos;
}

圧縮全文索引について

上述の通り、suffix arrayがあれば全文一致検索クエリを高速に処理できます。しかしsuffix arrayはふつうリファレンス文字列よりも大きい（上のコードだと4倍）ので、非常に長い文字列に対しては使えません。

一般に、全文検索を高速化するための追加データのことを全文索引と呼び、元の文字列よりも小さな全文索引のことを特に圧縮全文索引と呼びます。この記事で扱ったsuffix arrayは全文索引ですが圧縮全文索引ではありません。圧縮全文索引としてFM-index, (FM-index以外の)compressed suffix array(CSA), LZ圧縮ベースのindexingなどが知られています。FM-indexが特に高性能であると言われていてよく用いられます。
d.hatena.ne.jp
例えば上の記事でもそう書かれています。余談ですが、CSAという単語は(FM-index以外の)特定のアルゴリズムを示す場合と、FM-indexを含むある種のコンセプトを示す場合に分かれます。

特定のタスクにおいてCSAのほうが高性能だったと言っている論文もあって、「fm-index vs compressed suffix array」とかでググると見つかります。
Practical aspects of Compressed Suffix Arrays and FM-Index in Searching DNA Sequences - Semantic Scholar