SA-IS - omochan's diary

https://local.ugene.unipro.ru/tracker/secure/attachment/12144/Linear+Suffix+Array+Construction+by+Almost+Pure+Induced-Sorting.pdf
elementi-bioinformatica/Two Efficient Algorithms for Linear Time Suffix Array Construction.pdf at master · AlgoLab/elementi-bioinformatica · GitHub

SA-ISの論文2つだが、2つめのほうがわかりやすいし、実装例もあるので圧倒的にこっちをオススメする。

他にも役立ちそうなサイトをあげとく。
https://topcoder.g.hatena.ne.jp/iwiwi/20110720/1311168147
http://shogo82148.github.io/homepage/memo/algorithm/suffix-array/sa-is.html
http://blog.beam2d.net/2011/08/suffix-array-sa-is.html
http://gasin.hatenadiary.jp/entry/2017/06/08/101525

大まかなアルゴリズムの手順は、

1. Induced Sortを使ってLMS-substringをsortする。
2. LMS-substringのindexに関するSAを構築する。
3. 構築したSAを使って求めたい文字列に対するSAを構築する。

まずSuffix Arrayのアルゴリズムに関する前提知識がないとInduced Sortで引っかかる。Induced Sortもかなり高度なことをやっていて論文だけでは到底理解できない。
LMSとついている言葉が大量に出てきて(LMS-substring, LMS-suffix, LMS-prefix)これも混乱を招く要因。

論文を読むのは初めてではないといえどかなり苦戦した。
次からは言葉の定義をちゃんと抑えながら、適宜ググるのも忘れず、読み進めるようにしたい。

以下自分で実装したSA-IS(論文のコードを扱いやすくしただけだが)

namespace SA { //please init, don't use 0 as character
	vi bkt;
	inline bool isLMS(const vi& t, int i) {
		if(i > 0 && t[i] && !t[i - 1]) return true;
		else return false;
	}
	void getBuckets(const vi& s, int K, bool end) {
		int sum = 0, n = (int)s.size();
		fill(bkt.begin(), bkt.begin() + K, 0);
		rep(i, 0, n) bkt[s[i]]++;
		rep(i, 0, K) {
			sum += bkt[i];
			bkt[i] = end ? sum : sum - bkt[i];
		}
	}
	void induceSAl(const vi& t, const vi& s, vi& SA, int K) {
		int n = (int)SA.size();
		getBuckets(s, K, false);
		rep(i, 0, n) {
			int j = SA[i] - 1;
			if(j >= 0 && !t[j]) SA[bkt[s[j]]++] = j;
		}
	}
	void induceSAs(const vi& t, const vi& s, vi& SA, int K) {
		int n = (int)SA.size();
		getBuckets(s, K, true);
		rer(i, n, 0) {
			int j = SA[i] - 1;
			if(j >= 0 && t[j]) SA[--bkt[s[j]]] = j;
		}
	}
	vi get_sa_impl(const vi& s, int K) {
		int n = (int)s.size();
		vi t(n, 0);
		t[n - 1] = 1; t[n - 2] = 0;
		rer(i, n - 2, 0) 
			if(s[i] < s[i + 1] || (s[i] == s[i + 1] && t[i + 1] == 1)) t[i] = 1;
		getBuckets(s, K, true);
		vi SA(n, -1);
		rep(i, 1, n) if(isLMS(t, i)) SA[--bkt[s[i]]] = i;
		induceSAl(t, s, SA, K); induceSAs(t, s, SA, K);
		int n1 = 0;
		rep(i, 0, n) if(isLMS(t, SA[i])) SA[n1++] = SA[i];
		fill(SA.begin() + n1, SA.end(), -1);
		int name = 0, prev = -1;
		rep(i, 0, n1) {
			int pos = SA[i]; bool diff = false;
			rep(d, 0, n) {
				if(prev == -1 || s[pos + d] != s[prev + d] || t[pos + d] != 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;
			SA[n1 + pos] = name - 1;
		}
		vi s1(n1, -1), SA1(n1, -1), p1(n1, -1);
		int at = 0;
		rep(i, n1, n) if(SA[i] >= 0) s1[at++] = SA[i];
		if(name < n1) SA1 = get_sa_impl(s1, name);
		else rep(i, 0, n1) SA1[s1[i]] = i;
		getBuckets(s, K, true);
		at = 0;
		rep(i, 1, n) if(isLMS(t, i)) p1[at++] = i;
		fill(all(SA), -1);
		rer(i, n1, 0) {
			int j = p1[SA1[i]];
			SA[--bkt[s[j]]] = j;
		}
		induceSAl(t, s, SA, K); induceSAs(t, s, SA, K);
		return SA;
	}
	vi get_sa(const vi& sin, int K = 26) {
		int n = (int)sin.size();
		vi s(n + 1, 0);
		rep(i, 0, n) s[i] = sin[i] + 1;
		bkt.resize(max(n + 1, K + 1));
		return get_sa_impl(s, K + 1);
	}
}

無理やり圧縮しまくって81行になった。一応POJ3581でverify済み。
番兵とかも気を遣わず、元の文字列を数字に変換したのをそのままget_saの引数にすればSAが得られるようにした。
namespaceよりもstructにする方がいい気もするが、今のところnamespaceで困らないのでこれで。
(9/2追記)

rer(i, n1, 0) {
	int j = p1[SA1[i]]; //int j = p1[i]としてバグらせていた。
	SA[--bkt[s[j]]] = j;
}
induceSAl(t, s, SA, K); induceSAs(t, s, SA, K);

8/31にあげたコードがバグっていたので修正。なのになぜかPOJ3581は通ったんだよなぁ…。
デバッグの際にgasin氏のSAを活用させてもらった。gasinありがとー!