JavaScript: 拡張接尾辞配列

文字列を扱う上で拡張接尾辞配列というものがとても役立ちます。その実装は[esaxx](https://code.google.com/archive/p/esaxx)等で丸出しとなっています。本記事ではJavaScriptで書き換えて、web browser上でしゃぶり尽く方法を少しばかり紹介。
拡張接尾辞配列は接尾辞木を配列で構築したようなものです。接尾辞配列(SA)と最長共通接頭辞配列(LCP)を組み合わせる事で実装し、文字列検索等に絶大な効果を発揮します。
## 実装編
まず基本の接尾辞配列召喚program。線形時間で構築する事で有名な[SAIS](http://wpage.cocolog-nifty.com/wp/2010/08/sais-241.html)のJS版
```js:SAIS
/*
sais(T,SA[,n,k])
@T : input
@SA: output of suffix
@n : |T|
@k : max symbol in T[0..n-1]
* almost the same as saisSort(T,SA,0,n,k,0)

saisSort(T,SA,fs,n,k,isbwt)
@T : input
@SA: output of suffix
@fs: starting 0
@n : |T|
@k : max symbol in T[0..n-1]
@isbwt : true gives BWT array, false gives suffix array

sais.BWTe(T[,U,k])
@T : input
@U : output of BWT
@k : max symbol in T[0..n-1]
@return : BWT array if !U is true, otherwise primary index
*/
// find the start or end of each bucket
function getCounts(T,C,n,k,o){
	if(o)for(;k;)C[--k]=0;for(;n;)++C[T[--n]]}

function getBuckets(C,B,n,k,e){
	var i=0,s=0;
	if(e)for(;i<k;)B[i]=s+=C[i++];
	else for(;i<k;s+=e)e=C[i],B[i++]=s
}
// sort all type LMS suffixes
function LMSsort1(T,SA,C,B,n,k){
	C==B&&getCounts(T,C,n,k,1);
	getBuckets(C,B,n,k);
	var i=0,j=n-1,c0,c1=T[j],b=B[c1];
	for(SA[b++]=T[--j]<c1?~j:j;i<n;++i)
		if(0<(j=SA[i])){
			if((c0=T[j])!==c1)B[c1]=b,b=B[c1=c0];
			SA[b++]=T[--j]<c1?~j:j,SA[i]=0}
		else if(j<0)SA[i]=~j;
	C==B&&getCounts(T,C,n,k,1);
	getBuckets(C,B,n,k,1);
	for(b=B[c1=0];i;)
		if(0<(j=SA[--i])){
			if((c0=T[j])!==c1)B[c1]=b,b=B[c1=c0];
			SA[--b]=T[--j]>c1?~j-1:j,SA[i]=0}}

function LMSpostproc1(T,SA,n,m){
	for(var i=0,j,p,q=n,pl,ql=0,name=0,c0,c1,d;(p=SA[i])<0;)SA[i++]=~p;
	if(i<m)for(j=i;;)
		if((p=SA[++i])<0){
			SA[j++]=~p,SA[i]=0;
			if(j===m)break
	}
	for(c0=T[i=j=n-1];c1=c0,i--&&(c0=T[i])>=c1;);
	for(;~i;){
		for(;c1=c0,i--&&(c0=T[i])<=c1;);
		if(~i){
			SA[m+(i+1>>1)]=j-i;
			for(j=i+1;c1=c0,i--&&(c0=T[i])>=c1;);
		}
	}
	for(;++i<m;SA[j]=name){
		p=SA[i],d=pl=SA[j=m+(p>>1)];
		if(pl===ql&&q+pl<n){
			for(;d--&&T[p+d]===T[q+d];);++d
		}
		if(d)++name,q=p,ql=pl
	}return name
}
function LMSsort2(T,SA,C,B,D,n,k){
	getBuckets(C,B,n,k);
	var i=0,j=n-1,d=0,c0,c1=T[j],b=B[c1],t=T[--j]<c1;j+=n;
	for(SA[b++]=t&1?~j:j;i<n;++i)
		if(0<(j=SA[i])){
			if(n<=j)d++,j-=n;
			if((c0=T[j])!==c1)B[c1]=b,b=B[c1=c0];
			if(D[t=c0<<1|T[--j]<c1]!==d)j+=n,D[t]=d;
			SA[b++]=t&1?~j:j,SA[i]=0
		}else if(j<0)SA[i]=~j;
	for(;i;)if(0<SA[--i]&&SA[i]<n){
		for(SA[j=i]+=n;SA[--j]<n;);SA[i=j]-=n
	}
	getBuckets(C,B,n,k,1);
	for(i=n,d++,b=B[c1=0];i;)
		if(0<(j=SA[--i])){
			if(n<=j)d++,j-=n;
			if((c0=T[j])!==c1)B[c1]=b,b=B[c1=c0];
			if(D[t=c0<<1|T[--j]>c1]!==d)j+=n,D[t]=d;
			SA[--b]=t&1?~j-1:j,SA[i]=0
		}
}
function LMSpostproc2(SA,n,m){
	for(var i=0,j,d,name=0;(j=SA[i])<0;n>j||name++)SA[i++]=j=~j;
	if(i<m)for(d=i;;)
		if((j=SA[++i])<0){
			SA[d++]=j=~j;SA[i]=0;n>j||name++;
			if(d===m)break
		}
	i=m;
	if(name<m){
		for(d=name+1;i;SA[m+(j>>1)]=d)
			if(n<=(j=SA[--i]))j-=n,--d
	}else for(;i;)if(n<=(j=SA[--i]))SA[i]=j-n;
	return name
}
//compute SA and BWT
function induceSA(T,SA,C,B,n,k){
	C==B&&getCounts(T,C,n,k,1);
	getBuckets(C,B,n,k);
	var i=0,j=n-1,c0,c1=T[j],b=B[c1];
	for(SA[b++]=0<j&&T[j-1]<c1?~j:j;i<n;){
		j=SA[i],SA[i++]=~j;
		if(0<j){
			if((c0=T[--j])!==c1)B[c1]=b,b=B[c1=c0];
			SA[b++]=0<j&&T[j-1]<c1?~j:j
		}
	}
	C==B&&getCounts(T,C,n,k,1);
	getBuckets(C,B,n,k,1);
	for(b=B[c1=0];i;)
		if(0<(j=SA[--i])){
			if((c0=T[--j])!==c1)B[c1]=b,b=B[c1=c0];
			SA[--b]=!j||T[j-1]>c1?~j:j
		}else SA[i]=~j
}
function computeBWT(T,SA,C,B,n,k){
	C==B&&getCounts(T,C,n,k,1);
	getBuckets(C,B,n,k);
	var i=0,j=n-1,pidx=-1,c0,c1=T[j],b=B[c1];
	for(SA[b++]=0<j&&T[j-1]<c1?~j:j;i<n;++i)
		if(0<(j=SA[i])){
			SA[i]=~(c0=T[--j]);
			if(c0!==c1)B[c1]=b,b=B[c1=c0];
			SA[b++]=0<j&&T[j-1]<c1?~j:j
		}else if(j)SA[i]=~j;
	C==B&&getCounts(T,C,n,k,1);
	getBuckets(C,B,n,k,1);
	for(b=B[c1=0];i;)
		if(0<(j=SA[--i])){
			SA[i]=c0=T[--j];
			if(c0!==c1)B[c1]=b,b=B[c1=c0];
			SA[--b]=0<j&&T[j-1]>c1?~T[j-1]:j
		}else j?SA[i]=~j:pidx=i;
	return pidx
}
//find the suffix array SA of T[0..n-1]in {0..255}^n
function saisSort(T,SA,fs,n,k,isbwt){
	var C,B,D,b,i=0,j=n,m=0,p,q,name,c0=T[n-1],c1,f;
	if(q=k<257){
		C=new Int32Array(k);
		if(k>fs)B=new Int32Array(k),f=3;
		else B=SA.subarray(p=n+fs-k),f=1
	}else if(q=k>fs)C=B=new Int32Array(k),f=12;
	else{
		C=SA.subarray(n+fs-k);
		if(k<=fs-k)B=SA.subarray(p=n+fs-k-k),f=0;
		else if(k<1025)B=new Int32Array(k),f=2;
		else B=C,f=8
	}
	if(n<0x40000000&&k*2<=n)f|=f&1?2*k>fs-k?16:32:!f&&4*k<=fs&&32;
	getCounts(T,C,n,k,!q);
	for(getBuckets(C,B,n,k,1);i<n;)SA[i++]=0;i--;
	for(;c1=c0,i--&&(c0=T[i])>=c1;);
	for(;~i;){
		for(;c1=c0,i--&&(c0=T[i])<=c1;);
		if(~i){
			SA[b]=j,b=--B[c1],++m;
			for(j=i;c1=c0,i--&&(c0=T[i])>=c1;);
		}
	}
	if(m>1){
		if(f&48){
			D=f&16?new Int32Array(k*2):SA.subarray(p-k*2);
			B[T[j+1]]++;
			for(i=j=0;i<k;D[i]=D[i+++k]=0)
				if(B[i]!==(j+=C[i]))SA[B[i]]+=n;
			LMSsort2(T,SA,C,B,D,n,k),name=LMSpostproc2(SA,n,m)
		}else LMSsort1(T,SA,C,B,n,k),name=LMSpostproc1(T,SA,n,m)
	}else if(name=m)SA[b]=j+1;
	if(b=name<m){
		b=n+fs-m*2;
		if(f&4)C=null;if(f&2)B=null;
		if(!(f&13))k+name>b?f|=8:b-=k;
		i=m+(n>>1);D=SA.subarray(m+b);
		for(j=m;m<i;)if(c0=SA[--i])D[--j]=c0-1;
		saisSort(D,SA,b,m,name,0);
		for(c0=T[i=n-1];c1=c0,i--&&(c0=T[i])>=c1;);
		for(j=m;~i;){
			for(;c1=c0,i--&&(c0=T[i])<=c1;);
			if(0<=i)for(D[--j]=i+1;c1=c0,i--&&(c0=T[i])>=c1;);
		}
		for(;++i<m;)SA[i]=D[SA[i]];
		if(b=f&4)C=B=new Int32Array(k);
		if(f&2)B=new Int32Array(k)
	}
	f&8&&getCounts(T,C,n,k,!b);
	if(m>1){
		getBuckets(C,B,n,k,1),j=n,c1=T[p=SA[i=m-1]];
		do{
			for(q=B[c0=c1];q<j;)SA[--j]=0;
			for(;SA[--j]=p,i--&&(c1=T[p=SA[i]])===c0;);
		}while(~i);
		for(;j;)SA[--j]=0
	}
	if(isbwt)return computeBWT(T,SA,C,B,n,k);
	induceSA(T,SA,C,B,n,k)
}
function sais(T,SA,n,k){
	n=n>>>0||T.length;k=k>>>0||255;
	if(!T||!SA||!n)return-1;
	if(n<2)return n?SA[0]=0:0;
	if(n<3)return SA[n=0|T[0]<T[1]]=1,SA[1-n]=0;
	return saisSort(T,SA,0,n,k+1,0)
}
sais.BWTe=function(T,U,k){k=k>>>0||255;
	var n=T.length,i=n<257?1:n<65537?2:n<16777217?3:4,S=U||new Uint8Array(n+i),SA=new Int32Array(n);
	if(!T||!n)return[0,0];
	if(n<2){if(n)S[0]=T[0];return U?n:(S[1]=0,S)}
	k=saisSort(T,SA,0,i=n,++k,1);
	S[0]=T[n-1],S.set(SA.subarray(0,k++),1),S.set(SA.subarray(k),k);
	if(U)return k;
	for(k--;i--;k>>>=8)S[n++]=k&255;
	return S
};
sais.BWTd=function(A,p,k){k=k>>>0||255;
	var a=0,b=++k,c,B=A.length,s=0,C=[],P=[];
	if(!(p>-1))for(c=B<258?1:B<65539?2:B<16777220?3:4;c--;A.length--)p=(p<<8|A[--B])>>>0;
	for(;b;)C[--b]=0;
	for(;b<B;)P[b]=C[A[b++]]++;
	for(B=Array(b);s<b;s+=c)c=C[a],C[a++]=s;
	for(c=0;b;c<p&&c++)c=P[c]+C[B[--b]=A[c]];
	return B
};
```
そして次が拡張接尾辞配列の構築。空間消費量は最低でも入力の20倍となっております
```js:esaxx
/* eSA(T,SA,L,R,D,n,k)
@brief Build an enhanced suffix array of a given string in linear time.
For an input text T, eSA() builds an enhancd suffix array in linear time.
i-th internal node is represented as a triple(L[i],R[i],D[i]);
L[i] and R[i] is the left/right boundary of the suffix array as SA[L[i]....R[i]-1]
D[i] is the depth of the internal node.
The number of internal node is at most N-1 and return the actual number by
@param T[0...n-1]The input string.(random access iterator)
@param SA[0...n-1]The output suffix array(random access iterator)
@param L[0...n-1]The output left boundary of internal node(random access iterator)
@param R[0...n-1]The output right boundary of internal node(random access iterator)
@param D[0...n-1]The output depth of internal node(random access iterator)
@param n The length of the input string
@param k The alphabet size
@return The number of internal node
*/
function eSA(T,SA,L,R,D,n,k){
	if(n<0||k<1)return-1;
	if(saisSort(T,SA,0,n,k,0))return-2;
	return suffixtree(T,SA,L,R,D,n)
}
function suffixtree(T,SA,L,R,D,n){
	for(var h=0,i=n,l=n+1,p,r=0,s,t,S=[[-1,-1]];i--;)L[SA[i]]=i;
	for(R[0]=R[n]=i;++i<n;)if(s=L[i]){
		for(p=SA[s-1];T[i+h]===T[p+h];)++h;// i+h<n&&p+h<n
		R[s]=h&&h--
	}
	for(s=S[i=p=0];;s=S[++p]=[i,l-SA[i++]]){
		for(h=[i,t=R[i]];s[1]>t;s=S[--p]){
			if(i-s[0]>1)L[r]=s[0],R[r]=i,D[r++]=s[1];
			h[0]=s[0]
		}
		if(s[1]<t)S[++p]=h;
		if(i==n)break
	}return r
}
```
```js:利用例
let S="abracadabra";
let T=Uint8Array.from(S,a=>a.charCodeAt());
let n=T.length, SA=new Int32Array(n);
let L=new Int32Array(n), R=new Int32Array(n), D=new Int32Array(n);

eSA(T,SA,L,R,D,n,256);

for(let a=0,d;d=D[a];a++){
	let l=L[a], r=R[a], i=SA[l];
	document.write("pos: ",i,", length: ",d,", hits: ",r-l,"<pre><b>",S.substr(i,d),"</b></pre><hr>")
}
```
などといった感じで部分文字列を列挙できます。配列Dは部分文字列の長さのようなものが格納されていて、正数であれば部分文字列とみなせます。0が出現した時点で処理を打ち切って良い。
配列Lと配列Rには接尾辞配列における部分文字列の区間(`SA[L[i]...R[i]-1]`)が格納されています。つまり**L[i]-R[i]=共通部分文字列の出現数**とみなせます。
## 空間消費量削減版
前述の関数`suffixtree`を分割して配列の数を減らす。
```js
/* 最長共通接頭辞(Longest Common Prefix)構築
@T: Array/TypedArray. 文字列を数値化した配列
@SA: Array/Int32Array. 接尾辞配列(Suffix Array)
@H: Array/Int32Array. 高さ配列(n+1の大きさ)
@n: T.length
@return: H
*/
function makeLCP(T,SA,H,n){
	var h=0,i=n,p,s,ISA=new Int32Array(n),R=H||new Int32Array(n+1);
	for(;i--;)ISA[SA[i]]=i;
	for(R[0]=i;++i<n;)if(s=ISA[i]){
		for(p=SA[s-1];T[i+h]===T[p+h];)++h;// i+h<n&&p+h<n
		R[s]=h&&h--
	}return R
}
/* 接尾辞配列巡回
@T: String/Array/TypedArray. 文字列 or 文字列を数値化した配列
@SA: Array/Int32Array. 接尾辞配列(Suffix Array)
@H: Array/Int32Array. 高さ配列(T.length+1の大きさ)
@fn: 部分文字列列挙用関数
	fn(T,SA,rp,lp,len)
	@T: 部分文字列
	@SA: 接尾辞配列
	@rp: SA上の最後の出現位置
	@lp: SA上の最初の出現位置. rp-lp=出現数
	@len: 文字列長
*/
function traverseSA(T,SA,H,fn){
	for(let S=[[-1,-1]],h,i=0,n=SA.length,l=n+1,p=0,s=S[0],t;;s=S[++p]=[i,l-SA[i++]]){
		for(h=[i,t=H[i]];s[1]>t;s=S[--p]){
			i-s[0]>1&&fn(T,SA,i,s[0],s[1]);
			h[0]=s[0]
		}
		if(s[1]<t)S[++p]=h;
		if(i===n)break
	}
}
```
```js:for example
function output(T,SA,rp,lp,len){
	let i=SA[lp];
	document.write("pos: ",i,", length: ",len,", count: ",rp-lp,"<pre><b>",T.substr(i,len),"</b></pre><hr>")
}
let S="abracadabra", T=Uint8Array.from(S,a=>a.charCodeAt());
let n=T.length, SA=new Int32Array(n), H=new Int32Array(n+1);

saisSort(T,SA,0,n,256,0);
makeLCP(T,SA,H,n);
traverseSA(S,SA,H,output)
```
## 重複回避
このままだと部分文字列の両端が同一部分文字列と位置が重なっていても、出現回数が加算されます。例えば**ababa**には部分文字列`aba`が2箇所出現していますが、重なっている部分を除外すると1箇所となります(`"ababa".match(/aba/g).length==1`)。
それを回避したいなら一手間必要となります。素朴で馬鹿力任せな方法は、部分文字列の位置を整列して1つずつ重なっている分を間引く事。上述の`output`関数を以下のように修正。整列関数も追加。
```js
//挿入整列
function isort(A,a,z){
	for(var b=a,c,d,e;++a<z;A[d]=c)
		if((c=A[d=e=a])<A[b])for(;d>b;)A[d]=A[--d];
		else for(;c<A[--e];)A[d]=A[d=e]
}
//ternary quick sort
function tqsort(A,a,z){
	if(z-a<10)return isort(A,a,z+1);
	var c=A[a],d,i=A[z],j=a,k=z,p=A[a+z>>1];
	c>i?c<=p?p=c:p>i||(p=i):i<=p?p=i:c>p&&(p=c);
	for(i=a;j<=k;)
		if((c=A[j])>p)A[j]=A[k],A[k--]=c;
		else{if(c>p)d=A[i],A[i++]=c,A[j]=d;j++}
	a<--i&&tqsort(A,a,i);
	++k<z&&tqsort(A,k,z)
}
function output(T,SA,rp,lp,len){
	let i=SA[lp],hit=rp-lp,p;
	if(hit<3){
		p=i-SA[lp+1];
		if(p<0)p=-p;
		p<len&&hit-- //2地点の差が文字長未満=重なっている
	}else{
		let n=hit,I=[],c;
		for(p=0;hit;)I[--hit]=SA[lp+hit];//出現位置をcopy
		if(n<4){//3箇所の場合
			if(I[1]>I[2])c=I[1],I[1]=I[2],I[2]=c;
			if(I[0]>I[1]){c=I[0],I[0]=I[1],I[1]=c;
			if(c>I[2])I[1]=I[2],I[2]=c}
		}else n<10?isort(I,0,n):tqsort(I,0,n-1);//4箇所以上なら整列関数召喚
		//重なっている分を間引く
		for(;p<n;hit++)
			for(c=len+I[p];c>I[++p];);
	}
	document.write("pos: ",i,", length: ",len,", count: ",hit,"<pre><b>",T.substr(i,len),"</b></pre><hr>")
}
```
実験
<p class=codepen data-height="450" data-default-tab="result" data-slug-hash="GgJwZyO" data-pen-title="enumerate substrings by SA+LCP" data-user="xezz" style="height:300px;box-sizing:border-box;display:flex;align-items:center;justify-content:center;border:2px solid;margin:1em 0;padding:1em"> <span>See the Pen <a href="https://codepen.io/xezz/pen/GgJwZyO">enumerate substrings by SA+LCP</a> by xezz (<a href="https://codepen.io/xezz">@xezz</a>) on <a href="https://codepen.io">CodePen</a>.</span></p><script async src="https://public.codepenassets.com/embed/index.js"></script>