[ Suffix Array SA-IS Algorithm ] 後綴數組線性SA-IS算法

SA-IS是我目前看過最快的線性後綴數組演算法，但是做為競賽用途而進行簡化後他的效率在某些硬體上會比DC3慢，不過記憶體使用量是DC3的1/3~1/2倍，而最短的實現code也比DC3短很多，因此我認為這是十分優秀的算法

因為某些SA-IS的實作方式會利用傳入的陣列s的空間進行計算，因此傳入的陣列s必須是int範圍，而傳入的字串的尾端字元必須是最小的而且在之前完全沒有出現過，因此必須給字串先做以下處理:

假設要傳入字串 mississippi
先在尾端加入字元'#': mississippi#
算法結束後陣列sa為:11 10 7 4 1 0 9 8 6 3 5 2
將第一個數字11去除後剩下的數字即為字串mississippi的後綴數組:
10 7 4 1 0 9 8 6 3 5 2
關於SA-IS算法的論文請看這裡:
Two Efficient Algorithms for Linear Time Suffix Array Construction

關於SA-IS算法的教學(專利申請文)請看這裡:
后缀数组构造方法  CN 102081673 A

SA-IS的教學投影片:
https://www.cs.helsinki.fi/u/tpkarkka/opetus/11s/spa/lecture11.pdf

SA-IS中文教學:
https://riteme.github.io/blog/2016-6-19/sais.html

如果想對Suffix Array算法進行測試請使用這個Online Judge:
http://www.spoj.com/problems/SARRAY/

以下將會提供一些SA-IS的實做模板

1.台大黑暗code界的黑暗codebook:

	#include<string.h>
	#define MAGIC(XD){\
	memset(sa,0,sizeof(int)*n);\
	memcpy(x,c,sizeof(int)*z);\
	XD\
	memcpy(x+1,c,sizeof(int)*(z-1));\
	for(i=0;i<n;++i){\
	if(sa[i]&&!t[sa[i]-1])sa[x[s[sa[i]-1]]++]=sa[i]-1;\
	}\
	memcpy(x,c,sizeof(int)*z);\
	for(i=n-1;i>=0;--i){\
	if(sa[i]&&t[sa[i]-1])sa[--x[s[sa[i]-1]]]=sa[i]-1;\
	}\
	}
	void sais(int *s,int *sa,int *p,bool *t,int *c,int n,int z){
	bool neq=t[n-1]=1;
	int nn=0,nmxz=-1,*q=sa+n,*ns=s+n,*x=c+z,lst=-1,i,j;
	memset(c,0,sizeof(int)*z);
	for(i=0;i<n;++i)++c[s[i]];
	for(i=0;i<z-1;++i)c[i+1]+=c[i];
	for(i=n-2;i>=0;i--)t[i]=(s[i]==s[i+1]?t[i+1]:s[i]<s[i+1]);
	MAGIC(
	for(i=1;i<n;++i){
	if(t[i]&&!t[i-1])sa[--x[s[i]]]=p[q[i]=nn++]=i;
	}
	);
	for(i=0;i<n;++i)if((j=sa[i])&&t[j]&&!t[j-1]){
	neq=lst<0||memcmp(s+j,s+lst,(p[q[j]+1]-j)*sizeof(int));
	ns[q[lst=j]]=nmxz+=neq;
	}
	if(nmxz==nn-1)for(i=0;i<nn;++i)q[ns[i]]=i;
	else sais(ns,q,p+nn,t+n,c+z,nn,nmxz+1);
	MAGIC(
	for(i=nn-1;i>=0;--i)sa[--x[s[p[q[i]]]]]=p[q[i]];
	);
	}
	#undef MAGIC
	/*****************這些是需要用到的陣列大小**************/
	static const int MXN=10000000;
	int s[MXN*2+5],sa[MXN*2+5],p[MXN+5],c[MXN*2+5];
	bool t[MXN*2+5];

view raw SA-IS black magic.cpp hosted with ❤ by GitHub

其實這段code本來是被壓縮得更短，而且用到非常多的記憶體，不過經過卦長的改良後成功減少記憶體的使用並將他排版成正常人能看得懂的樣子
而這裡的MXN則是字串的最長長度(假設字元數<字串長度)

2.卦長自行實作的模板(記憶體用量少):

	#define pushS(x) sa[--bkt[s[x]]] = x
	#define pushL(x) sa[bkt[s[x]]++] = x
	#define induce_sort(v){\
	fill_n(sa,n,0);\
	copy_n(_bkt,A,bkt);\
	for(i=n1-1;~i;--i)pushS(v[i]);\
	copy_n(_bkt,A-1,bkt+1);\
	for(i=0;i<n;++i)if(sa[i]&&!t[sa[i]-1])pushL(sa[i]-1);\
	copy_n(_bkt,A,bkt);\
	for(i=n-1;~i;--i)if(sa[i]&&t[sa[i]-1])pushS(sa[i]-1);\
	}
	template<typename T>
	void sais(const T s,int n,int *sa,int *_bkt,int *p,bool *t,int A){
	int *rnk=p+n,*s1=p+n/2,*bkt=_bkt+A;
	int n1=0,i,j,x=t[n-1]=1,y=rnk[0]=-1,cnt=-1;
	for(i=n-2;~i;--i)t[i]=(s[i]==s[i+1]?t[i+1]:s[i]<s[i+1]);
	for(i=1;i<n;++i)rnk[i]=t[i]&&!t[i-1]?(p[n1]=i,n1++):-1;
	fill_n(_bkt,A,0);
	for(i=0;i<n;++i)++_bkt[s[i]];
	for(i=1;i<A;++i)_bkt[i]+=_bkt[i-1];
	induce_sort(p);
	for(i=0;i<n;++i)if(~(x=rnk[sa[i]]))
	j=y<0||memcmp(s+p[x],s+p[y],(p[x+1]-p[x])*sizeof(s[0]))
	,s1[y=x]=cnt+=j;
	if(cnt+1<n1)sais(s1,n1,sa,bkt,rnk,t+n,cnt+1);
	else for(i=0;i<n1;++i)sa[s1[i]]=i;
	for(i=0;i<n1;++i)s1[i]=p[sa[i]];
	induce_sort(s1);
	}
	/*****************這些是需要用到的陣列大小**************/
	const int MAXN=10000005,MAXA='z'+1;
	int sa[MAXN],bkt[MAXN+MAXA],p[MAXN*2];
	bool t[MAXN*2];
	char s[MAXN];

view raw SA-IS.cpp hosted with ❤ by GitHub

為了簡化實作方法及減少記憶體的使用，因此將計算後剩下的空間進行重複利用，壓縮後只需要這些記憶空間，而傳入的陣列s可以直接傳入char陣列，因此對使用者來說是非常方便的一份code
這裡的MXN則是字串的最長長度(假設字元數<字串長度)

3.論文提供的實作code:

	/*SA-IS Algorithm*/
	const unsigned char mask[]={0x80,0x40,0x20,0x10,0x08,0x04,0x02,0x01};
	#define tget(i) ( (t[(i)/8]&mask[(i)%8])?1:0)
	#define tset(i, b) t[(i)/8]=(b)?(mask[(i)%8]|t[(i)/8]):((~mask[(i)%8])&t[(i)/8])
	#define chr(i) (cs==sizeof(int)?((int*)s)[i]:((unsigned char *)s)[i])
	#define isLMS(i) (i>0&&tget(i)&&!tget(i-1))
	// find the start or end of each bucket
	inline void getBuckets(unsigned char *s,int *bkt,int n,int K,int cs,bool end){
	int i,sum=0;
	// clear all buckets
	for(i=0;i<=K;++i)bkt[i]=0;
	// compute the size of each bucket
	for(i=0;i<n;++i)++bkt[chr(i)];
	for(i=0;i<=K;++i){
	sum+=bkt[i];
	bkt[i]=end?sum:sum-bkt[i];
	}
	}
	// compute SAl
	inline void induceSAl(unsigned char *t,int *SA,unsigned char *s,int *bkt,int n,int K,int cs,bool end){
	int i,j;
	// find starts of buckets
	getBuckets(s,bkt,n,K,cs,end);
	for(i=0;i<n;++i){
	j=SA[i]-1;
	if(j>=0&&!tget(j))SA[bkt[chr(j)]++]=j;
	}
	}
	// compute SAs
	inline void induceSAs(unsigned char *t,int *SA,unsigned char *s,int *bkt,int n,int K,int cs,bool end){
	int i,j;
	// find ends of buckets
	getBuckets(s,bkt,n,K,cs,end);
	for(i=n-1;i>=0;--i){
	j=SA[i]-1;
	if(j>=0&&tget(j))SA[--bkt[chr(j)]]=j;
	}
	}
	// find the suffix array SA of s[0..n-1] in {1..K}?n
	// require s[n-1]=0 (the sentinel!), n>=2
	// use a working space (excluding s and SA) of
	// at most 2.25n+O(1) for a constant alphabet
	void SA_IS(unsigned char *s,int *SA,int n,int K,short cs){
	// LS-type array in bits
	unsigned char *t=(unsigned char *)malloc(n/8+1);
	int i,j;
	// classify the type of each character
	// the sentinel must be in s1, important!!!
	tset(n-2,0);tset(n-1,1);
	for(i=n-3;i>=0;--i)
	tset(i,(chr(i)<chr(i+1)||(chr(i)==chr(i+1)&&tget(i+1)==1))?1:0);
	// stage 1: reduce the problem by at least 1/2
	// sort all the S-substrings
	// bucket array
	int *bkt=(int *)malloc(sizeof(int)*(K+1));
	// find ends of buckets
	getBuckets(s,bkt,n,K,cs,true);
	for(i=0;i<n;++i)SA[i]=-1;
	for(i=1;i<n;++i)if(isLMS(i))SA[--bkt[chr(i)]]=i;
	induceSAl(t,SA,s,bkt,n,K,cs,false);
	induceSAs(t,SA,s,bkt,n,K,cs,true);
	free(bkt);
	// compact all the sorted substrings into
	// the first n1 items of SA
	// 2*n1 must be not larger than n (proveable)
	int n1=0;
	for(i=0;i<n;++i)if(isLMS(SA[i]))SA[n1++]=SA[i];
	// find the lexicographic names of substrings
	// init the name array buffer
	for(i=n1;i<n;++i)SA[i]=-1;
	int name=0,prev=-1;
	for(i=0;i<n1;++i){
	int pos=SA[i];bool diff=false;
	for(int d=0;d<n;++d)
	if(prev==-1||chr(pos+d)!=chr(prev+d)||tget(pos+d)!=tget(prev+d)){
	diff=true;break;
	}else if(d>0&&(isLMS(pos+d)||isLMS(prev+d)))break;
	if(diff){
	++name;prev=pos;
	}
	pos=(pos%2==0)?pos/2:(pos-1)/2;
	SA[n1+pos]=name-1;
	}
	for(i=n-1,j=n-1;i>=n1;--i)if(SA[i]>=0)SA[j--]=SA[i];
	// stage 2: solve the reduced problem
	// recurse if names are not yet unique
	int *SA1=SA,*s1=SA+n-n1;
	if(name<n1)SA_IS((unsigned char*)s1,SA1,n1,name-1,sizeof(int));
	else // generate the suffix array of s1 directly
	for(i=0;i<n1;++i)SA1[s1[i]]=i;
	// stage 3: induce the result for
	// the original problem
	// bucket array
	bkt=(int *)malloc(sizeof(int)*(K+1));
	// put all the LMS characters into their buckets
	// find ends of buckets
	getBuckets(s,bkt,n,K,cs,true);
	for(i=1,j=0;i<n;++i)if(isLMS(i))s1[j++]=i; // get p1
	// get index in s
	for(i=0;i<n1;++i)SA1[i]=s1[SA1[i]];
	// init SA[n1..n-1]
	for(i=n1;i<n;++i)SA[i]=-1;
	for(i=n1-1;i>=0;--i){
	j=SA[i];SA[i]=-1;
	SA[--bkt[chr(j)]]=j;
	}
	induceSAl(t,SA,s,bkt,n,K,cs,false);
	induceSAs(t,SA,s,bkt,n,K,cs,true);
	free(bkt);free(t);
	}

view raw SA-IS paper code.c hosted with ❤ by GitHub

這是其中一個比DC3快的code，而且記憶體使用量是最少的，但是長度很長就是了，不適合在比賽時使用

4.超快記憶體使用超少的模板庫code:
https://gist.github.com/jacky860226/1d33adad858eef71bfe18120d8d69e6d#file-sa-is-very-fast-cpp
因為長度太長所以就直接貼上網址了，沒有人會在比賽時寫這種東西

[ Suffix Array DC3 algoruthm ] 後綴數組線性DC3演算法

這邊提供後綴數組DC3算法的模板，他是一個能在線性時間$\ord N$內求得後綴數組的演算法，注意這邊傳入的陣列s及sa長度必須是字串長的3倍以上，因為DC3會用到大量的記憶體，而且會利用傳入的陣列s的空間進行計算，因此傳入的陣列s必須是int範圍，而傳入的字串的尾端字元必須是最小的而且在之前完全沒有出現過，因此必須給字串先做以下處理:

假設要傳入字串 mississippi
先在尾端加入字元'#': mississippi#
算法結束後陣列sa為:11 10 7 4 1 0 9 8 6 3 5 2
將第一個數字11去除後剩下的數字即為字串mississippi的後綴數組:
10 7 4 1 0 9 8 6 3 5 2

因為是遞迴的關係，所以常數很大，但在大數據的時候比倍增法快
關於DC3演算法的介紹請看這兩篇文章:
http://ching.im/post/algorithm/difference-cover-modulo-algorithm
《后缀数组——处理字符串的有力工具》

如果想對Suffix Array算法進行測試請使用這個Online Judge:
http://www.spoj.com/problems/SARRAY/

在傳入的參數中最後一個A為最大的字元+1通常設為'z'+1

以下提供模板(其實code挺短的):

	#define F(x) (x)/3+((x)%3==1?0:tb)
	#define G(x) (x)<tb?(x)*3+1:((x)-tb)*3+2
	int wa[3*maxn+5],wb[3*maxn+5],wv[3*maxn+5],c[maxn+5];
	inline bool c0(int *s,int a,int b){
	return s[a]==s[b]&&s[a+1]==s[b+1]&&s[a+2]==s[b+2];
	}
	inline bool c12(int k,int *s,int a,int b){
	if(k==2)return s[a]<s[b]||s[a]==s[b]&&c12(1,s,a+1,b+1);
	else return s[a]<s[b]||s[a]==s[b]&&wv[a+1]<wv[b+1];
	}
	inline void radix_sort(int *s,int *a,int *b,int len,int A){
	static int i;
	for(i=0;i<len;++i)wv[i]=s[a[i]];
	for(i=0;i<A;++i)c[i]=0;
	for(i=0;i<len;++i)++c[wv[i]];
	for(i=1;i<A;++i)c[i]+=c[i-1];
	for(i=len-1;i>=0;--i)b[--c[wv[i]]]=a[i];
	}
	void dc3(int *s,int *sa,int len,int A){
	int i,j,*san=sa+len,ta=0,tb=(len+1)/3,tbc=0,p;
	int *rn=s+len;
	s[len]=s[len+1]=0;
	for(i=0;i<len;++i){
	if(i%3!=0)wa[tbc++]=i;
	}
	radix_sort(s+2,wa,wb,tbc,A);
	radix_sort(s+1,wb,wa,tbc,A);
	radix_sort(s,wa,wb,tbc,A);
	for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;++i){
	rn[F(wb[i])]=c0(s,wb[i-1],wb[i])?p-1:p++;
	}
	if(p<tbc)dc3(rn,san,tbc,p);
	else for(i=0;i<tbc;i++)san[rn[i]]=i;
	for(i=0;i<tbc;++i){
	if(san[i]<tb)wb[ta++]=san[i]*3;
	}
	if(len%3==1)wb[ta++]=len-1;
	radix_sort(s,wb,wa,ta,A);
	for(i=0;i<tbc;++i)wv[wb[i]=G(san[i])]=i;
	for(i=0,j=0,p=0;i<ta&&j<tbc;++p){
	sa[p]=c12(wb[j]%3,s,wa[i],wb[j])?wa[i++]:wb[j++];
	}
	for(;i<ta;++p)sa[p]=wa[i++];
	for(;j<tbc;++p)sa[p]=wb[j++];
	}
	#undef F
	#undef G
	/*
	DC3的輸入陣列s必須能儲存0~len的數值，所以設為int
	s及sa陣列的大小必須大於3*maxn
	*/

view raw DC3.cpp hosted with ❤ by GitHub

[ Suffix Array Prefix doubling algorithms ] 後綴數組倍增算法

後綴數組(又稱尾碼陣列)是一個十分強大的字串處理武器，大部分的問題都可以用它來解決，它可以幾乎做到所有後綴樹(Suffix Tree)能做到的事，所以這邊就不介紹後綴樹了

因為後綴數組可以由後綴樹進行遍歷轉換而來，而構造後綴樹僅需花費線性的時間，所以構造後綴數組的時間可為線性$\ord N$，但是後綴數組本身就是為了減少構造後綴樹的空間與代碼量而被發明出來的，直接由後綴樹轉換是沒意義的
但是仍然有其他線性構造後綴數組的方法，像是DC3、SA-IS等會在下一篇介紹，這次要講的是比較簡單常用的方法-----倍增法
關於後綴數組的使用說明可以參考《后缀数组——处理字符串的有力工具》
關於倍增法的說明可以參考 演算法筆記-SuffixArray 的部分
這邊提供$\ord{N*logN*logN}$及$\ord{NlogN}$的模板

$\ord{N*logN*logN}$:

	#include<algorithm>
	struct CMP{
	int len,k,*rank,a,b;
	inline bool operator()(int i,int j){
	if(rank[i]!=rank[j])return rank[i]<rank[j];
	a=(i+=k)<len?rank[i]:-1;
	b=(j+=k)<len?rank[j]:-1;
	return a<b;
	}
	};
	inline void suffix_array(const char *s,int len,int *sa,int *rank,int *tmp){
	for(int i=0;i<len;++i){
	sa[i]=i;
	rank[i]=s[i];
	}
	CMP cmp={len,1};
	for(;;cmp.k<<=1){
	cmp.rank=rank;
	std::sort(sa,sa+len,cmp);
	tmp[sa[0]]=0;
	for(int i=1;i<len;++i)tmp[sa[i]]=tmp[sa[i-1]]+cmp(sa[i-1],sa[i]);
	if(tmp[sa[len-1]]==len-1)break;
	std::swap(rank,tmp);
	}
	}

view raw doubling algorithm NlogNlogN.cpp hosted with ❤ by GitHub

$\ord{NlogN}$:

	#include<algorithm>
	#define radix_sort(x,y){\
	for(i=0;i<A;++i)c[i]=0;\
	for(i=0;i<len;++i)c[x[y[i]]]++;\
	for(i=1;i<A;++i)c[i]+=c[i-1];\
	for(i=len-1;i>=0;--i)sa[--c[x[y[i]]]]=y[i];\
	}
	#define sac(r,a,b) r[a]!=r[b]||a+k>=len||r[a+k]!=r[b+k]
	inline void suffix_array(const char *s,int len,int *sa,int *rank,int *tmp,int *c){
	int A='z'+1,i,k,id=0;
	for(i=0;i<len;++i)rank[tmp[i]=i]=s[i];
	radix_sort(rank,tmp);
	for(k=1;id<len-1;k<<=1){
	for(id=0,i=len-k;i<len;++i)tmp[id++]=i;
	for(i=0;i<len;++i)if(sa[i]>=k)tmp[id++]=sa[i]-k;
	radix_sort(rank,tmp);
	std::swap(rank,tmp);
	for(rank[sa[0]]=id=0,i=1;i<len;++i)
	rank[sa[i]]=id+=sac(tmp,sa[i-1],sa[i]);
	A=id+1;
	}
	}
	#undef radix_sort
	#undef sac

view raw doubling algorithm NlogN.cpp hosted with ❤ by GitHub

注意此方法必須要在字元集大小為常數的情況下有效，否則必須離散化

所需的陣列長度只需要與字串陣列相同即可
當然$N*logN*logN$的做法會比較值觀，$NlogN$的方法則是利用radix_sort進行的，radix_sort本來在倍增的時候要先排序第一關鍵字跟第二關鍵字，但是第二關鍵字排序的結果可以用已經求好的SA直接求出來
對radix_sort還不了解的人請看這個網頁:https://www.cs.usfca.edu/~galles/visualization/RadixSort.html
如果想對Suffix Array算法進行測試請使用這個Online Judge:
http://www.spoj.com/problems/SARRAY/

英文字母大小寫轉換特殊做法

假設有一個題目是這樣的:
        給定一串英文字母，請將大寫的部分轉成小寫，小寫的部分轉成大寫並輸出

一般我們會用if或是三元運算子做判斷，但是這太麻煩了
經過觀察發現，摁合一個小寫字母ascii與其對應的大寫字母ascii相差皆為32，
而其二進位編碼剛好允許透過 xor 32的方式進行轉換，但是32這個數字不好記，又可以發現ascii 32='空白'，而以下的code就可以將大小寫互換:

1 2 3 4 5 6 7 8

#include<stdio.h> char s[1000005]; int main(){ scanf("%s",s); for(int i=0;s[i];++i)putchar(s[i]^' ');puts(""); for(int i=0;s[i];++i)putchar(s[i]^32);puts(""); return 0; }

兩種寫法會有相同的效果