MP法

MP法を勉強したのでメモ。

文字列の頭良い感じの線形アルゴリズムたち - あなたは嘘つきですかと聞かれたら「YES」と答えるブログ

これを参考にしました。

A[i]の定義はS[0:i-1]の接頭辞と接尾辞が最大何文字一致しているかである。
[0:i-1]が一致しているからそんなの常にiになるんじゃねと思うが、接頭辞と接尾辞は文字列自身を含まないのでこの定義でOK。

これを求めて何がいいかというと、例えば文字列Tの中でSがあるか判定するプログラムで、S[0:i-1]で一致しているが、S[i]で不一致の時、次はSの何文字目から見ればいいかがわかる。

具体的には、
T=AGCAGCA…
S=AGCAGCG 
を先頭から比較していくと、i=6で不一致が発生する。T[3,6],S[0,3]は一致しているので、次見るべきはi=4からということになる。
(実際のMP法とは少し違うが)

iまでAが求まっているとする。またA[i]=kとする。A[i+1]を求めるために、まずS[0:k],S[i-k:i]が一致しているか見る。
しかしA[i]=kからS[0:k-1],S[i-k:i-1]が一致することは既知なので実際はS[k]とS[i]を比較すればよい。
もし一致していればA[i]=k+1だが一致していなければ次見るところはA[k]になる。なぜならS[0:k]でうまくいかない以上、次候補となるのは、二番目に最初からのsubstrと後ろからのsubstrが一致しているところである。それはまさしくA[k]であるからだ。

KMPを使って、蟻本p327「禁止文字列」の前処理partを行ってみる。
S[i]=AGCT[j]なら単にnext[i][j]=i+1とすれば良くて
S[i]≠AGCT[j]ならk=A[i]としてmp法の要領でS[k]とAGCT[j]が一致している場所を探す。

しかしこのアルゴリズムでは最悪計算量O(K^2)である。Aばかりの文字列を考えてもらってみるとわかると思う。
next[i]['G']などを求めるのにいちいちO(K)かかる。

蟻本にはO(K)でできる方法があると書いてあるが、これはおそらくKMP法を使ったら達成できるので、朝になったらやってみます。

vector<int> MP(string S) { //it doesn't include itself
	int n = (int)S.size();
	vector<int> A(n + 1);
	A[0] = -1;
	rep(i, 0, n) {
		int j = A[i];
		while (j >= 0 && S[i] != S[j]) j = A[j];
		j++;
		A[i + 1] = j;
	}
	return A;
}

void solve() {
	N = 3; S = "AT"; K = (int)S.size();
	vector<int> A = MP(S);
	debug(A);
	rep(i, 0, K) {
		rep(j, 0, 4) {
			if(S[i] == AGCT[j]) nex[i][j] = i + 1;
			else {
				int k = A[i];
				while(k >= 0 && S[k] != AGCT[j]) k = A[k];
				k++;
				nex[i][j] = k;
			}
		}
	}

	dp[0][0] = 1;
	rep(i, 1, K) dp[0][i] = 0;

	rep(t, 0, N) {
		rep(i, 0, K) dp[t + 1][i] = 0;

		rep(i, 0, K) {
			rep(j, 0, 4) {
				int ti = nex[i][j];
				if(ti == K) continue;
				ADD(dp[t + 1][ti], dp[t][i]);
			}
		}
	}

	int ans = 0;
	rep(i, 0, K) ADD(ans, dp[N][i]);
	cout << ans << "\n";
}