JavaScript: 接尾辞配列(Suffix Array)編

概要

これは文字列の接尾辞(開始位置が異なり終端位置が元の文字列と同じ部分文字列)の文字列中の開始位置を要素とする配列を、接尾辞に関して辞書順に並べ替えて得られる単純な配列です。まさに接尾辞木の配列版と言っても過言でしかありませんね。文字列検索等において猛威を振るう事が可能です。
例えばthatの接尾辞はthat, hat, at, tの4個で、これらを辞書順に並べた配列が接尾辞配列となるのです。

始点	接尾辞
0	that
1	hat
2	at
3	t

整列

始点	接尾辞
2	at
1	hat
3	t
0	that

厳密に言うと始点こそが接尾辞配列の正体です。上記の例で言えば[0,1,2,3]から並べ替えられた[2,1,3,0]です。
昇順に並んでいるので文字列の検索は二分探索を使ってウミウシ並に高速に行うことができます。ただし整列処理はかなりの力仕事だという事を忘れてはいけません。
それではこの辺で本題のprogramに話をすり替えていきます。Array.sortを使いこなせば容易です。

function getSA(S){
	for(var i=0,n=S.length,SA=[];i<n;)SA[i]=S.substr(i,n-i++);//接尾辞をかき集めて…
	for(SA.sort();i;)SA[--i]=n-SA[i].length;//整列後に始点に変換
	return SA
}
console.log(getSA("banana"))

原理に則り非常にくそ真面目で馬鹿力任せな雰囲気丸出しですね。惚れ惚れしちゃいます。でも超簡単お手軽なのはいいですが、これで巨大文字列を処理しようとすると最悪の場合死に至ります。空間消費量が核爆発的に増えていくのは明らかで、ほとんど使い物になりません。
全面的に修正です。文字列ではなく最初から文字符号に変換しておきます。文字列の複製は行わず、整列関数に特製比較関数を食わせておけば前述の問題は解消します。

function getSA(S,n){
	return new Uint32Array(n=S.length).map((a,b)=>b).sort((a,b,c,d)=>{
		for(c=a<b?n-b:n-a;c--;)if(d=S[a++]-S[b++])return d;return b-a
	})
}
console.log(getSA([98,97,110,97,110,97]))

記述量削減版

function getSA(S,n){
	return[...Array(n=S.length).keys()].sort((a,b,c,d)=>{
		for(c=a<b?n-b:n-a;c--;)if(d=S[a++]-S[b++])return d;return b-a
	})
}

少ない文字列を連続で処理するなら非常に役立つでしょう。しかし巨大文字列を処理するにはまだ心許ないのです。
それはともかく、何やら1文字の変数名だらけで、可読性が皆無のような… いえ、気のせいです。気にしてはいけません。この先ずっとこんな調子です。

multikey quicksortによる高速化

手抜きはこの辺にして、御登場頂くのはmultikey quicksort様です。ただのquick sortと違い、これは文字単位で比較を行うので接尾辞の整列に適しているのです。

//subsort用比較関数
function compare(T,a,b,i,n,r){
	for(n-=a<b?b:a;i<n&&!(r=T[a+i]-T[b+i]);)i++;return r||b-a
}
//小区間用shellsort
function subsort(T,A,l,h,d,n){
	if(10<h-l)
		for(var i=l+2,j,m=i,t;i<h;A[j]=t)
			for(t=A[j=++i];m<j&&compare(T,t,A[j-3],d,n)<0;)A[j]=A[j-=3];
	for(i=l;i<h;A[j]=t)
		for(t=A[j=++i];l<j&&compare(T,t,A[j-1],d,n)<0;)A[j--]=A[j]
}
//中央値献上
function median(T,A,a,b,c,d,z){
	var r=A[a]+d,s=A[b]+d;
	r=r<z?T[r]:-1,s=s<z?T[s]:-1,d+=A[c],d=d<z?T[d]:-1;
	return r<s?s<d?b:r<d?c:a:r<d?a:s<d?c:b
}
//枢軸選択
function pivot(T,A,r,l,h,d,z){
	if(r<65)var m=l+(r>>1),p;
	else r>>=3,
		l=median(T,A,l,p=l+r,p+=r,d,z),
		m=median(T,A,p+=r,p+=r,p+=r,d,z),
		h=median(T,A,p+=r,p+=r,h,d,z);
	p=A[m=median(T,A,l,m,h,d,z)];A[m]=A[l];A[l]=p;
	return T[p+d]
}
//multikey quicksort
function mqsort(T,A,l,h,d,n){
	for(var a,b,i,j,o,p,r,s=0,t,L=[],H=[],D=[];;){
		if(h-l<20){
			l<h&&subsort(T,A,l,h,d,n);
			if(!s)break;
			l=L[--s];h=H[s];d=D[s];
			continue}
		p=pivot(T,A,h-l+1,l,h,d,n);
		for(i=a=l,j=b=h;;A[j--]=t){
			for(;i<=j&&(r=((t=A[i]+d)<n?T[t]:-1))<=p;i++)
				r^p||(A[i]=A[a],A[a++]=t-d);
			for(;i<=j&&(r=((t=A[j]+d)<n?T[t]:-1))>=p;j--)
				r^p||(A[j]=A[b],A[b--]=t-d);
			if(i>j)break;
			t=A[i],A[i++]=A[j]}
		if((p=a-l)<(r=i-a))r=p;o=j;
		for(p=l;r--;A[o--]=t)t=A[p],A[p++]=A[o];
		if((p=h-b)<(r=b-j))r=p;o=h;
		for(p=i;r--;A[o--]=t)t=A[p],A[p++]=A[o];
		i+=l-a;j-=~h+b;
		if(j<=h)L[s]=j,H[s]=h,D[s++]=d;
		if(i<j)L[s]=i,H[s]=j-1,D[s++]=d+1;
		h=i-1
	}
}
/*
@T: Array/Uint8Array. 整列元の数値配列
@A: Array/Uint32Array/Int32Array. 整列後の配列
@n: 配列Tの要素数. 省略可
*/
function getSA(T,A,n){
	for(var i=n=n||T.length,SA=A||new Uint32Array(n);i;)SA[--i]=i;
	mqsort(T,SA,0,n-1,0,n);
	return SA
}

mqsortは再帰しないようにしました。他にも工夫として分割範囲が狭い場合はshellsort君を召喚しています。枢軸も範囲が狭ければ3か所、広ければ9か所から厳選。見事なものです。これを元に更に高速化してみましょう。
最初に接尾辞の先頭1～2文字分だけ分布数え整列するのです。そしてmqsortでとどめを刺します。

//先頭1文字で整列
function o1dsort(T,A,C,n,cs){
	for(var i=n;i;)++C[T[A[--i]=i]];
	for(;i<cs;)C[i+1]+=C[i++];
	for(;n;)A[--C[T[--n]]]=n
}
//先頭2文字で整列
function o2dsort(T,A,C,n,cs){
	for(var i=n,c;i;)++C[c|(c=T[A[--i]=i])<<8];
	for(;i<cs;)C[i+1]+=C[i++];
	for(c=0;n;)A[--C[c|(c=T[--n])<<8]]=n
}
/*
@T: Array/Uint8Array. 整列元の数値配列
@A: Array/Uint32Array/Int32Array. 整列後の配列
@n: 配列Tの要素数. 省略可
@return: 配列TのSuffix Array
*/
function getSA(T,A,n){n=n||T.length;
	var SA=A||new Int32Array(n);
	if(n<2)return SA;
	if(n<3)return SA[n=0|T[0]<T[1]]=1,SA[1-n]=0,SA;
	if(n<300)return mqsort(T,SA,0,n-1,0,n),SA;
	var a,b=0,c=n>>14?65536:256,d=c>>9?2:1,i=0,C=new Uint32Array(c+1);
	(c>>9?o2dsort:o1dsort)(T,SA,C,n,c);
	for(;b<n;b-a>1&&mqsort(T,SA,a,b-1,d,n))a=C[i++],b=C[i]; //区間ごとに整列. 累積頻度が配列Tの大きさに達したら終了
	return SA
}

上記の例では配列Tの大きさが300未満ならmqsortを1回呼んで終了、16384未満なら先頭1文字、それ以上では先頭2文字で分布数え整列を施してからmqsortを連呼します。入力配列が小さいと分布数え整列の労力も無視できないので、そんな感じに分けた方が良いです。mqsort等の実装は前述の通り。
これでかなり高速化しますが、弱点があります。同じ文字が広範囲に続いたり、繰り返しの多い文字列の処理に非常に時間がかかるのです。それを克服しないと実用性ほぼ皆無のごみです。
そこで通称Larsson Sadakane法と呼ばれる有名な手口を利用します。

function compare(R,a,b,c,n,r,i){
	for(i=c,n-=a<b?b:a;i<n&&!(r=R[a+i]-R[b+i]);)i+=c;return r||b-a
}
function subsort(R,A,l,h,d,n){
	if(10<h-l)
		for(var i=l+2,j,m=i,t;i<h;A[j]=t)
			for(t=A[j=++i];m<j&&compare(R,t,A[j-3],d,n)<0;)A[j]=A[j-=3];
	for(i=l;i<h;A[j]=t)
		for(t=A[j=++i];l<j&&compare(R,t,A[j-1],d,n)<0;)A[j--]=A[j]
}
function median(R,A,a,b,c,d,z){
	var r=A[a]+d,s=A[b]+d;
	r=r<z?R[r]:-1,s=s<z?R[s]:-1,d+=A[c],d=d<z?R[d]:-1;
	return r<s?s<d?b:r<d?c:a:r<d?a:s<d?c:b
}
function pivot(R,A,r,l,h,d,z){
	if(r<65)var m=l+(r>>1),p;
	else r>>=3,
		l=median(R,A,l,p=l+r,p+=r,d,z),
		m=median(R,A,p+=r,p+=r,p+=r,d,z),
		h=median(R,A,p+=r,p+=r,h,d,z);
	p=A[m=median(R,A,l,m,h,d,z)];//A[m]=A[l];A[l]=p;
	return R[p+d]
}
function mqsort(R,I,l,h,d,n){
	if(h-l<20){
		for(l<h&&subsort(R,I,l,h,d,n);l<=h;)I[R[I[l]]=l++]=-1;return}
	for(var a=l,b=h,i=l,j=h,k,t,o,p=pivot(R,I,h-l+1,l,h,d,n);;I[j--]=k){
		for(;i<=j&&(k=((t=I[i]+d)<n?R[t]-p:-1))<1;i++)
			k||(I[i]=I[a],I[a++]=t-d);
		for(;i<=j&&(k=((t=I[j]+d)<n?R[t]-p:-1))>=0;j--)
			k||(I[j]=I[b],I[b--]=t-d);
		if(i>j)break;
		k=I[i],I[i++]=I[j]}
	if((p=a-l)<(k=i-a))k=p;
	for(p=l,o=j;k--;I[o--]=t)t=I[p],I[p++]=I[o];
	if((p=h-b)<(k=b-j))k=p;
	for(p=i,o=h;k--;I[o--]=t)t=I[p],I[p++]=I[o];
	i+=l-a;j-=~h+b;
	l<i&&mqsort(R,I,l,i-1,d,n);
	if(j>i)if(j-i<2)I[R[I[i]]=i]=-1;
		else for(p=j-1;i<j;)R[I[i++]]=p;
	j>h||mqsort(R,I,j,h,d,n)
}
function o1dsort(T,C,F,I,R,n,l){
	for(var c,i=n;i;)++F[T[--i]];
	for(i=0;i<l;F[i]+=c)C[i+1]=c=F[i++];
	for(i=n;i;)R[I[--F[c=T[--i]]]=i]=C[++c]-1
}
function o2dsort(T,C,F,I,R,n,l){
	for(var c,i=n;i;)F[c|(c=T[--i])<<8]++;
	for(;i<l;F[i]+=c)C[i+1]=c=F[i++];
	for(c=0,i=n;i;)R[I[--F[n=c|(c=T[--i])<<8]]=i]=C[++n]-1
}
/* Larsson Sadakane法
@T: Array/Uint8Array. 整列元の数値配列
@A: Array/Int32Array. 整列後の配列
@n: 配列Tの要素数. 省略可
@return: 配列TのSuffix Array
*/
function LSsort(T,A,n){n=n||T.length;
	var I=A||new Int32Array(n);
	if(n<2)return n?I[0]=0:0,I;
	if(n<3)return I[n=0|T[0]<T[1]]=1,I[1-n]=0,I;
	var c,f,i,l=n>>14?65536:256,F=new Uint32Array(l),C=new Uint32Array(l+1),R=new Int32Array(n);
	(l>>9?o2dsort:o1dsort)(T,C,F,I,R,n,l);
	for(l=l>>9?2:1;f=l<n;l*=2){ //doubling technique
		for(i=0;i<n;i<n&&mqsort(R,I,i,f=R[I[i]],l,n,i=f+1,f=0)){
			for(c=i;c<n&&I[c]<0;)c-=I[c];
			if(i<c)I[i]=i-(i=c)}
		if(f)break}
	for(;I[R[--n]]=n;);
	return I
}

これでもう心配ゴム用、実用的なprogramの完成です。入力の大きさをnとした時、空間消費9n、計算量はO(n log2 n)程度となります。しかし理論上は接尾辞配列はO(n)で構築できます。この程度で満足してはいけません。

Induced Sorting

2008年頃にGe Nong、Sen Zhang、Wai Hong Chanらが発明した画期的な手法です。原作だとC言語による記述で約100行という有り様です。単純ながらも非常に巧妙な手口で小さくまとまっています。
しかしこのままだとbinary fileに対応できないので、ほぼ使い物になりません。少し工夫して対応するとともに、無駄な空間消費となる配列t(LS-type array)も消し去ってやりましょう。

//find the start or end of each bucket
function getBuckets(T,B,n,K,end,o){
	var i=K,s=n;
	if(o)for(;i;)B[--i]=0;
	for(i=0;s;)B[T[--s]]++;//compute the size of each bucket
	if(end)for(;s<n;)B[i]=s+=B[i++];
	else for(;s<n;s+=K)K=B[i],B[i++]=s
}
//compute SA
function induceSA(SA,T,B,n,K){
	//compute SAl
	getBuckets(T,B,n,K,0,1);//find starts of buckets
	SA[B[T[n-1]]++]=n-1;
	for(var i=0,j;i<n;j<1||T[j]>T[--j]||(SA[B[T[j]]++]=j))j=SA[i++];
	//compute SAs
	getBuckets(T,B,n,K,1,1);//find ends of buckets
	for(;i;j<1||T[j]<T[--j]||(SA[--B[T[j]]]=j))j=SA[--i]
}
/*
@T: Array/Uint8Array. 整列元の数値配列
@A: Array/Int32Array. 整列後の配列
@n: 配列Tの要素数
@K: 配列Tの要素の最大値
*/
//find the suffix array SA of T[0..n-1]in {0..K}^n. use a working space(excluding T and SA)of at most 2n+O(1)for a constant alphabet
function sais(T,SA,n,K){
	//stage 1:reduce the problem by at least 1/2. sort all the S-substrings
	var i=0,j,c=0,p,n1=0,name=0,m=n,l,h=0,B=new Int32Array(K+1);
	//find ends of buckets
	for(getBuckets(T,B,n,K,1);i<n;)SA[i++]=0;
	for(i--;i;)
		if(T[i--]+c>T[i])c=1;
		else if(c)SA[--B[T[i+1]]]=i+1,c=0;
	induceSA(SA,T,B,n,K);
	//compact all the sorted substrings into the first n1 items of SA. 2*n1 must be not larger than n(proveable)
	for(;i<n;)if(p=SA[i++],0<p&&T[p-1]>(c=T[p])){
		for(j=p+1;j<n&&T[j]===c;)++j;
		if(j<n&&c<T[j])SA[n1++]=p
	}
	for(c=0;i>n1;)SA[--i]=0;//init the name array buffer
	//store the length of all substrings
	for(i=j=n-1;i;)
		if(T[i--]+c>T[i])c=1;
		else if(c)SA[n1+(i+1>>1)]=j-i,j=i,c=0;
	//find the lexicographic names of all substrings
	for(;i<n1;SA[j]=name){
		p=SA[i++],c=l=SA[j=n1+(p>>1)];
		if(l===h){for(;c--&&T[p+c]===T[m+c];);++c}
		if(c)++name,m=p,h=l;
	}
	//stage 2:solve the reduced problem
	//recurse if names are not yet unique
	if(m=name<i){
		for(i=j=n;j>n1;)if(p=SA[--j])SA[--i]=p-1;
		B=null;sais(p=SA.subarray(n-j,n),SA,j,name-1);
		for(i=n-1,c=0;i;)
			if(T[i--]+c>T[i])c=1;
			else if(c)p[--j]=i+1,c=0;
		for(B=new Int32Array(K+1);i<n1;)SA[i]=p[SA[i++]]//get index in T
	}
	//stage 3:induce the result for the original problem
	//put all left-most p characters into their buckets
	for(getBuckets(T,B,n,K,1,!m);i<n;)SA[i++]=0;
	for(i=n1;i;SA[--B[T[j]]]=j)j=SA[--i],SA[i]=0;
	induceSA(SA,T,B,n,K)
}

こんなものでしょう。ついでにBucketを解放したり再確保したりする時機も、再帰呼び出しが発生する時に限定。他にもJavaScriptの仕様上、getBucketsでBucketの初期化を省略できる時は省略。

もっと高速な奴

べらぼうに高速なImproved two-stage suffix sorting + rank sort
超高速なlibsaisのJavaScript移植版