关于网上搜查得到的3DC3的基于字符串后缀数组的排序方法的怀疑

最近整理笔试题，发现了在字符串处理中后缀数组的应用，随机翻看了下；网上最多的内容全部是应用一个叫“罗穗骞”的孩子的国家集训论文里面的内容，

其中比较全面讲解了后缀数组的应用问题，但是，就其3DC3的排序方法中我总觉得有些问题，于是手动测试了一下，确实弄出了毛病，先看看他的讲解：

一、DC3算法

         算法分3步：

　　（1）、先将后缀分成两部分，然后对第一部分的后缀排序。

　 　将后缀分成两部分，第一部分是后缀k（k模3不等于0），第二部分是后缀k（k模3等于0）。先对所有起始位置模3不等于0的后缀进行排序，即对 suffix(1)，suffix(2)，suffix(4)，suffix(5)，suffix(7)……进行排序。做法是将suffix(1)和 suffix(2)连接，如果这两个后缀的长度不是3的倍数，那先各自在末尾添0使得长度都变成3的倍数。然后每3个字符为一组，进行基数排序，将每组字符“合并”成一个新的字符。然后用递归的方法求这个新的字符串的后缀数组。如图3所示。在得到新的字符串的sa后，便可以计算出原字符串所有起始位置模3不等于0的后缀的sa。要注意的是，原字符串必须以一个最小的且前面没有出现过的字符结尾，这样才能保证结果正确（请读者思考为什么）。



　　（2）、利用（1）的结果，对第二部分的后缀排序。

　　剩下的后缀是起始位置模3等于0的后缀，而这些后缀都可以看成是一个字符加上一个在（1）中已经求出 rank的后缀，所以只要一次基数排序便可以求出剩下的后缀的sa。

　　（3）、将（1）和（2）的结果合并，即完成对所有后缀排序。

　　这个合并操作跟合并排序中的合并操作一样。每次需要比较两个后缀的大小。分两种情况考虑，第一种情况是suffix(3 * i)和suffix(3 * j + 1)的比较，可以把suffix(3 * i)和suffix(3 * j + 1)表示成：

suffix(3 * i) = r[3 * i] + suffix(3 * i + 1)

suffix(3 * j + 1) = r[3 * j + 1] + suffix(3 * j + 2)

　　其中 suffix(3 * i + 1)和 suffix(3 * j + 2)的比较可以利用（2）的结果快速得到。第二种情况是 suffix(3 * i)和 suffix(3 * j + 2)的比较，可以把 suffix(3 * i)和suffix(3 * j + 2)表示成：

suffix(3 * i) = r[3 * i] + r[3 * i + 1] + suffix(3 * i + 2)　　

suffix(3 * j + 2) = r[3 * j + 2] + r[3 * j + 3] + suffix(3 * (j + 1) + 1)

　　同样的道理，suffix(3 * i + 2)和 suffix(3 * (j + 1) + 1)的比较可以利用（2）的结果快速得到。所以每次的比较都可以高效的完成，这也是之前要每 3个字符合并，而不是每 2个字符合并的原因。

　　具体实现：

    #define F(x) ((x)/3+((x)%3==1?0:tb))

    #define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)

    int wa[maxn],wb[maxn],wv[maxn],ws[maxn];

    int c0(int *r,int a,int b)

    {return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}

    int c12(int k,int *r,int a,int b)

    {if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);

    else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}

    void sort(int *r,int *a,int *b,int n,int m)

    {

        int i;

        for(i=0;i<n;i++) wv[i]=r[a[i]];

        for(i=0;i<m;i++) ws[i]=0;

        for(i=0;i<n;i++) ws[wv[i]]++;

        for(i=1;i<m;i++) ws[i]+=ws[i-1];

        for(i=n-1;i>=0;i--) b[--ws[wv[i]]]=a[i];

        return;

    }

    void dc3(int *r,int *sa,int n,int m)

    {

        int i,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;

        r
=r[n+1]=0;

        for(i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;

        sort(r+2,wa,wb,tbc,m);

        sort(r+1,wb,wa,tbc,m);

        sort(r,wa,wb,tbc,m);

        for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)

        rn[F(wb[i])]=c0(r,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(n%3==1) wb[ta++]=n-1;

        sort(r,wb,wa,ta,m);

        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,r,wa[i],wb[j])?wa[i++]:wb[j++];

        for(;i<ta;p++) sa[p]=wa[i++];

        for(;j<tbc;p++) sa[p]=wb[j++];

        return;

    }

　　各个参数的作用和前面的倍增算法一样，不同的地方是r数组和sa数组的大小都要是3*n，这为了方便下面的递归处理，不用每次都申请新的内存空间。函数中用到的变量：

    int i,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;

　 　rn数组保存的是（1）中要递归处理的新字符串，san数组是新字符串的sa。变量ta表示起始位置模3为0的后缀个数，变量tb表示起始位置模3为1 的后缀个数，已经直接算出。变量tbc表示起始位置模3为1或2的后缀个数。先按（1）中所说的用基数排序把3个字符“合并 ”成一个新的字符。为了方便操作，先将r
和r[n+1]赋值为0。

　　代码：

    r
=r[n+1]=0;

    for(i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;

    sort(r+2,wa,wb,tbc,m);

    sort(r+1,wb,wa,tbc,m);

    sort(r,wa,wb,tbc,m);

　　其中sort函数的作用是进行基数排序。代码：

    void sort(int *r,int *a,int *b,int n,int m)

    {

        int i;

        for(i=0;i<n;i++) wv[i]=r[a[i]];

        for(i=0;i<m;i++) ws[i]=0;

        for(i=0;i<n;i++) ws[wv[i]]++;

        for(i=1;i<m;i++) ws[i]+=ws[i-1];

        for(i=n-1;i>=0;i--) b[--ws[wv[i]]]=a[i];

        return;

    }

　　基数排序结束后，新的字符的排名保存在 wb数组中。

　　跟倍增算法一样，在基数排序以后，求新的字符串时要判断两个字符组是否完全相同。代码：

    for(p=1,rn[F(wb[0])]=0,i=1; i<tbc;i++)

    rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;

　　其中 F(x)是计算出原字符串的 suffix(x)在新的字符串中的起始位置，c0函数是比较是否完全相同，在开头加一段代码：

    #define F(x) ((x)/3+((x)%3==1?0:tb))

    inline int c0(int *r,int a,int b)

    {return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}

　　接下来是递归处理新的字符串，这里和倍增算法一样，可以加一个小优化，如果p等于tbc，那么说明在新的字符串中没有相同的字符，这样可以直接求出san数组，并不用递归处理。代码：

    if(p<tbc) dc3(rn,san,tbc,p);

    else for(i=0;i<tbc;i++) san[rn[i]]=i;

　　然后是第（2）步，将所有起始位置模3等于0的后缀进行排序。其中对第二关键字的排序结果可以由新字符串的sa直接计算得到，没有必要再排一次。代码：

    for(i=0;i<tbc;i++) if(san[i]<tb) wb[ta++]=san[i]*3;

    if(n%3==1) wb[ta++]=n-1;

    sort(r,wb,wa,ta,m);

    for(i=0;i<tbc;i++) wv[wb[i]=G(san[i])]=i;

　　要注意的是，如果n%3==1，要特殊处理suffix(n-1)，因为在san数组里并没有suffix(n)。G(x)是计算新字符串的suffix(x)在原字符串中的位置，和F(x)为互逆运算。在开头加一段：

    #define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)。

　　最后是第（3）步，合并所有后缀的排序结果，保存在sa数组中。代码：

    for(i=0,j=0,p=0;i<ta && j<tbc;p++)

    sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];

    for(;i<ta;p++) sa[p]=wa[i++];

    for(;j<tbc;p++) sa[p]=wb[j++];

   

　　其中c12函数是按（3）中所说的比较后缀大小的函数，k=1是第一种情况，k=2是第二种情况。代码：

    int c12(int k,int *r,int a,int b)

    {if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);

    else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}

　　整个DC3算法基本写好，代码大约40行。

　　算法分析：

　　假设这个算法的时间复杂度为f(n)。容易看出第（1）步排序的时间为O(n)（一般来说，m比较小，这里忽略不计），新的字符串的长度不超过2n/3，求新字符串的 sa的时间为f(2n/3)，第（2）和第（3）步的时间都是 O(n)。所以

f(n) ＝ O(n) + f(2n/3)

f(n) ≤ c×n + f(2n/3)

f(n) ≤ c×n + c×(2n/3) + c×(4n/9) + c×(8n/27) + …… ≤ 3c×n

所以 f(n) = O(n)

以上是他论文中的讲解，但是我用了一个程序测了一组数据，发现存在问题：

#include<stdio.h>
#define maxn 20
#define F(x) ((x)/3+((x)%3==1?0:tb))
#define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)
int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
int c0(int *r,int a,int b)
{return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}
int c12(int k,int *r,int a,int b)
{if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}
void sort(int *r,int *a,int *b,int n,int m)
{
int i;
for(i=0;i<n;i++) wv[i]=r[a[i]];
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[wv[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) b[--ws[wv[i]]]=a[i];
return;
}
void dc3(int *r,int *sa,int n,int m)
{
int i,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;
r
=r[n+1]=0;
for(i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;
sort(r+2,wa,wb,tbc,m);
sort(r+1,wb,wa,tbc,m);
sort(r,wa,wb,tbc,m);
for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)
rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;
if(p<tbc) {
for (int f = 0; f < tbc; f++)
{
printf("%d ",rn[f]);
}
printf("\n");
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(n%3==1) wb[ta++]=n-1;
sort(r,wb,wa,ta,m);
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,r,wa[i],wb[j])?wa[i++]:wb[j++];
for(;i<ta;p++) sa[p]=wa[i++];
for(;j<tbc;p++) sa[p]=wb[j++];
return;
}

int main(int agrc, char **argv)
{
int my_testcase[60] = {4,2,5,5,2,4,3,2,1,2,1,3,4,1,2,1,0,0};
int my_sa[60];
dc3(my_testcase,my_sa,16,6);
for (int i = 0; i < 16; i++) {
printf("%d ", my_sa[i]);
}
return 0;
}