1 #include <iostream> 2 #include <cstring> 3 #include <cstddef> 4 #include <cstdio> 5 #include <string> 6 #include <algorithm> 7 const int maxn = 100001; 8 int wa[maxn],wb[maxn],wv[maxn],ws[maxn]; 9 int cmp(int *rank, int a,int b,int l) 10 { 11 return rank[a]==rank[b] && rank[a+l]==rank[b+l]; 12 } 13 14 void da(int *r,int *sa,int n,int m) 15 { 16 int i, k, p, *x=wa, *y=wb, *t; 17 18 for(i=0;i<m;i++) ws[i] = 0; 19 for(i=0;i<n;i++) ws[x[i] = r[i]]++; 20 for(i=1;i<m;i++) ws[i] += ws[i-1]; 21 for(i=n-1;i>=0;i--) sa[--ws[x[i]]] = i; 22 23 for(k=1, p=1; p<n; k*=2, m=p) 24 { 25 p=0; for(i=n-k;i<n;i++) y[p++]=i; 26 for(i=0;i<n;i++) if(sa[i]>=k) y[p++]=sa[i]-k; 27 28 for(i=0;i<n;i++) wv[i]=x[y[i]]; 29 30 for(i=0;i<m;i++) ws[i]=0; 31 for(i=0;i<n;i++) ws[wv[i]]++; 32 for(i=1;i<m;i++) ws[i]+=ws[i-1]; 33 for(i=n-1;i>=0;i--) sa[--ws[wv[i]]]=y[i]; 34 35 t=x,x=y,y=t; 36 for(p=1,x[sa[0]]=0,i=1;i<n;i++) 37 x[sa[i]]=cmp(y,sa[i-1],sa[i],k)?p-1:p++; 38 } 39 return; 40 } 41 42 int rank[maxn],height[maxn]; 43 void calheight(int *r,int *sa,int n) 44 { 45 int i,j,k=0; 46 for(i=1;i<=n;i++) rank[sa[i]]=i; 47 for(i=0;i<n;height[rank[i++]]=k) 48 for(k?k--:0,j=sa[rank[i]-1];r[i+k]==r[j+k];k++) 49 ; 50 return; 51 } 52 53 int main() 54 { 55 int r[] = {1, 1, 2, 1, 1, 1, 1, 2}, sa[maxn]; 56 int n = 8; 57 r[8] = 0; 58 da(r,sa,n+1,128); 59 calheight(r,sa,n); 60 }
小清新版
1 #include <iostream> 2 #include <cstring> 3 #include <cstddef> 4 #include <cstdio> 5 #include <string> 6 #include <algorithm> 7 8 const int maxn = 100001; 9 int wa[maxn],wb[maxn],wv[maxn],ws[maxn]; 10 /* 11 wa[]: 本意是保存各个后缀的rank值的,但是这里并没有去存储rank值,因为后续只是涉及wa[]的比较工作, 12 因而这一步可以不用存储真实的rank值,能够反映相对的大小即可。 13 wb[]: 存放的是按第二关键字排序的子串首字符下标 14 wv[]: 存放每个子串的第一关键字 15 ws[]: 存放每个rank值的数目 16 */ 17 int cmp(int *rank, int a,int b,int l) 18 { 19 return rank[a]==rank[b] && rank[a+l]==rank[b+l]; 20 } 21 /* 22 比较函数,合并的子串相同则返回1,不同返回0; 23 就像论文所说,由于末尾填了0,所以如果r[a]==r[b](实际是y[a]==y[b]), 24 说明待合并的两个长为j的字符串,前面那个一定不包含末尾0,因而后面这个的起始位置至多在0的位置,不会再靠后了,因而不会产生数组越界。 25 */ 26 27 /* 28 *r: 数组(字符串) 29 *sa: 后缀数组 30 n: 字符串中字符的个数,注意这里的n里面是包括人为在字符串末尾添加的那个0的 31 m: 字符串中字符的取值范围,是基数排序的一个参数,如果原序列都是字母可以直接取128, 32 如果原序列本身都是整数的话,则m可以取比最大的整数大1的值。 33 */ 34 void da(int *r,int *sa,int n,int m) 35 { 36 int i, k, p, *x=wa, *y=wb, *t; 37 /* 38 x[]里面本意是保存各个后缀的rank值的,但是这里并没有去存储rank值,因为后续只是涉及x[]的比较工作, 39 因而这一步可以不用存储真实的rank值,能够反映相对的大小即可。 40 y[]里存放的是按第二关键字排序的字符串下标 41 *t 作交换指针 42 */ 43 44 45 /*以下四行代码是把长度为1的子串进行基数排序*/ 46 /*如果不理解为什么这样可以达到基数排序的效果,建议自己实际用纸笔模拟一下*/ 47 for(i=0;i<m;i++) ws[i] = 0; 48 for(i=0;i<n;i++) ws[x[i] = r[i]]++; 49 for(i=1;i<m;i++) ws[i] += ws[i-1]; 50 for(i=n-1;i>=0;i--) sa[--ws[x[i]]] = i; 51 /*ws[]数组对原字符串中各字符的数目进行统计 52 第1行清零; 53 第2行上面已经提到x[]保存的是后缀的相对rank值,x[i] = r[i]的意思是将x[i]初始化为各字符的值,字符的ASCII值也就可以代表长度为1的子串的相对顺序 54 第3行的作用即求出最后一个子串i的rank是多少,供第4行使用 55 第4行相当于从后向前得到各子串的sa[]数组,i之所以从n-1开始循环,是为了保证在当字符串中有相等的字符串时,默认靠前的字符串更小一些。*/ 56 for(int i = 0; i < n; i++) 57 { 58 printf("sa[%d]: %d\n", i, sa[i]); 59 } 60 61 /*下面这层循环中p代表rank值不同的字符串的数量,也可以这么理解,所有子串排序后,相等的子串rank值相同,则rank的范围是[1,p]; 62 如果p达到n,即各后缀的rank与sa已全部求出;因为长度不一,所以不可能出现相等的情况; 63 k代表当前待合并的字符串的长度,每次将两个长度为k的字符串合并成一个长度为2k的字符串; 64 m同样代表基数排序的元素的取值范围*/ 65 for(k=1, p=1; p<n; k*=2, m=p) 66 { 67 /*以下两行代码实现对第二关键字的排序*/ 68 /*所谓第二关键字即子串r[i+k, k]与r[j+k, k]的rank比较*/ 69 p=0; for(i=n-k;i<n;i++) y[p++]=i; 70 /*结合论文的插图,我们可以看到位置在第n-k至n的元素的第二关键字都为0,因此如果按第二关键字排序,必然这些元素都是排在前面的(第二关键字为0即无法构成以r[i]为首字符的长度为2k的子串)。*/ 71 /*y[]里存放的是按第二关键字排序的子串首字符下标*/ 72 for(i=0;i<n;i++) if(sa[i]>=k) y[p++]=sa[i]-k; 73 /*结合论文的插图,我们可以看到,下面一行的第二关键字不为0的部分都是根据上一次的排序结果得到的, 74 且上一行中只有sa[i]>=k的第sa[i]个字符串 75 (这里以及后面指的“第?个字符串”不是按字典序排名来的,是按照首字符在字符串中的位置来的) 76 的rank才会作为下一行的第sa[i]-j个字符串的第二关键字,而且显然按sa[i]的顺序rank[sa[i]]是递增的, 77 因此完成了对剩余的元素的第二关键字的排序。*/ 78 79 printf("k: %d\n", k); 80 for(int i = 0; i < p; i++) 81 { 82 printf("y[%d] : %d\n", i, y[i]); 83 } 84 85 for(i=0;i<n;i++) wv[i]=x[y[i]]; 86 /*这里相当于提取出每个字符串的第一关键字(前面说过了x[]是保存rank值的,也就是字符串的第一关键字),放到wv[]里面是方便后面的使用*/ 87 88 /*以下四行代码是按第一关键字进行的基数排序*/ 89 /* 90 wv[]: 存放每个子串的第一关键字,wv[i] = x的含义为按第二关键字第i小的子串的第一关键字的值 91 ws[]: 存放每个rank值的数目 92 */ 93 for(int i = 0; i < n; i++) 94 { 95 printf("x[y[%d]]: %d\n", i, x[y[i]]); 96 } 97 for(i=0;i<m;i++) ws[i]=0; 98 for(i=0;i<n;i++) ws[wv[i]]++; 99 for(i=1;i<m;i++) ws[i]+=ws[i-1]; 100 for(i=n-1;i>=0;i--) sa[--ws[wv[i]]]=y[i]; 101 /*此处巧妙地将第一关键字与第二关键字结合起来了,注意理解*/ 102 103 /*i之所以从n-1开始循环,含义同上,同时注意这里是y[i],因为y[i]里面才存着字符串的下标*/ 104 105 /*下面两行就是计算合并之后的rank值了,而合并之后的rank值应该存在x[]里面,但我们计算的时候又必须用到上一层的rank值, 106 也就是现在x[]里面放的东西,如果我既要从x[]里面拿,又要向x[]里面放,怎么办? 107 当然是先把x[]的东西放到另外一个数组里面,省得乱了。这里就是用交换指针的方式,高效实现了将x[]的东西“复制”到了y[]中。*/ 108 t=x,x=y,y=t; 109 for(p=1,x[sa[0]]=0,i=1;i<n;i++) 110 x[sa[i]]=cmp(y,sa[i-1],sa[i],k)?p-1:p++; 111 112 /* 113 这里就是用x[]存储计算出的各字符串rank的值了,记得我们前面说过,计算sa[]值的时候如果字符串相同是默认前面的更小的, 114 但这里计算rank的时候必须将相同的字符串看作有相同的rank,要不然p==n之后就不会再循环啦 115 p的值表示的是此时关键字不同的串的数量 116 cmp比较函数,合并的子串相同则返回1,不同返回0; 117 注意p和i的初始值需为1,因为循环中存在i-1和p-1,而x[sa[0]]的值也需初始化为0 118 */ 119 } 120 return; 121 } 122 123 //能够线性计算height[]的值的关键在于h[](height[rank[]])的性质,即h[i]>=h[i-1]-1,下面具体分析一下这个不等式的由来。 124 //论文里面证明的部分一开始看得我云里雾里,后来画了一下终于搞明白了,我们先把要证什么放在这:对于第i个后缀,设j=sa[rank[i] - 1],也就是说j是i的按排名来的上一个字符串,按定义来i和j的最长公共前缀就是height[rank[i]],我们现在就是想知道height[rank[i]]至少是多少,而我们要证明的就是至少是height[rank[i-1]]-1。 125 //好啦,现在开始证吧。 126 //首先我们不妨设第i-1个字符串(这里以及后面指的“第?个字符串”不是按字典序排名来的,是按照首字符在字符串中的位置来的)按字典序排名来的前面的那个字符串是第k个字符串,注意k不一定是i-2,因为第k个字符串是按字典序排名来的i-1前面那个,并不是指在原字符串中位置在i-1前面的那个第i-2个字符串。 127 //这时,依据height[]的定义,第k个字符串和第i-1个字符串的公共前缀自然是height[rank[i-1]],现在先讨论一下第k+1个字符串和第i个字符串的关系。 128 //第一种情况,第k个字符串和第i-1个字符串的首字符不同,那么第k+1个字符串的排名既可能在i的前面,也可能在i的后面,但没有关系,因为height[rank[i-1]]就是0了呀,那么无论height[rank[i]]是多少都会有height[rank[i]]>=height[rank[i-1]]-1,也就是h[i]>=h[i-1]-1。 129 //第二种情况,第k个字符串和第i-1个字符串的首字符相同,那么由于第k+1个字符串就是第k个字符串去掉首字符得到的,第i个字符串也是第i-1个字符串去掉首字符得到的,那么显然第k+1个字符串要排在第i个字符串前面,要么就产生矛盾了。同时,第k个字符串和第i-1个字符串的最长公共前缀是height[rank[i-1]],那么自然第k+1个字符串和第i个字符串的最长公共前缀就是height[rank[i-1]]-1。 130 //到此为止,第二种情况的证明还没有完,我们可以试想一下,对于比第i个字符串的字典序排名更靠前的那些字符串,谁和第i个字符串的相似度最高(这里说的相似度是指最长公共前缀的长度)?显然是排名紧邻第i个字符串的那个字符串了呀,即sa[rank[i]-1]。也就是说sa[rank[i]]和sa[rank[i]-1]的最长公共前缀至少是height[rank[i-1]]-1,那么就有height[rank[i]]>=height[rank[i-1]]-1,也即h[i]>=h[i-1]-1。 131 //证明完这些之后,下面的代码也就比较容易看懂了。 132 int rank[maxn],height[maxn]; 133 void calheight(int *r,int *sa,int n) 134 { 135 int i,j,k=0; 136 for(i=1;i<=n;i++) rank[sa[i]]=i; //计算每个字符串的字典序排名 137 for(i=0;i<n;height[rank[i++]]=k) //将计算出来的height[rank[i]]的值,也就是k,赋给height[rank[i]]。i是由0循环到n-1,但实际上height[]计算的顺序是由height[rank[0]]计算到height[rank[n-1]]。 138 for(k?k--:0,j=sa[rank[i]-1];r[i+k]==r[j+k];k++); //上一次的计算结果是k,首先判断一下如果k是0的话,那么k就不用动了,从首字符开始看第i个字符串和第j个字符串前面有多少是相同的,如果k不为0,按我们前面证明的,最长公共前缀的长度至少是k-1,于是从首字符后面k-1个字符开始检查起即可。 139 return; 140 } 141 142 //最后再说明一点,就是关于da和calheight的调用问题,实际上在“小罗”写的源程序里面是如下调用的,这样我们也能清晰的看到da和calheight中的int n不是一个概念,同时height数组的值的有效范围是height[1]~height[n]其中height[1]=0,原因就是sa[0]实际上就是我们补的那个0,所以sa[1]和sa[0]的最长公共前缀自然是0。 143 int main() 144 { 145 int r[] = {1, 1, 2, 1, 1, 1, 1, 2}, sa[maxn]; 146 int n = 8; 147 r[8] = 0; 148 da(r,sa,n+1,128); 149 calheight(r,sa,n); 150 }
【转】带注释版
时间: 2024-11-09 06:15:56