后缀数组——倍增算法

最近在做字符串方面的训练，做到后缀数组时卡了好久，后缀数组的实现有两种方法，而我看的是倍增算法，其实倍增算法的思路还是很简单的（详见白书），只是给的模板代码对于我这种渣渣来说有点艰难，弄了好久才弄懂。

虽说是字符串方面的算法，但是它的实现过程是以数字的形式进行的，也正因如此，该算法同样适用于数字串（字符串应先转换成ASCII码）

在使用倍增算法时应注意的三点：

1，如果是求一个连接串的后缀数组时，连接符应该用一个不可能出现的数字

2，在串的末尾加一个元素零，而且其他元素要保证大于零（比如让其他元素全部自增1）

3，在求sa数组时，第零个元素一定的存的是最后一个序号（因为其他都大于0，所以最后一位的后缀一定是字典序最小的），同理，Rank数组的最后一个元素肯定是0，其实简单点说就是补上的那个零除了用来求sa数组以外，其他大多少情况下都是当作无效位（比如求height数组），数组的具体意思见代码

下面对模板做个详细的分析：

//str是要求的串，sa就是后缀数组，sa[i]存的是字典序排名第i的后缀序号，Rank存的是后缀的排名,Rank[i]存的是第i个后缀的排名，
//height是保存排名相邻的后缀（sa[Rank[i]]和sa[Rank[i]-1]）的最长公共前缀的长度,剩下的是辅助数组
int str[MAXN],height[MAXN],Rank[MAXN],a[MAXN],b[MAXN],c[MAXN];
void getSa(int top,int Max)//top是str的长度，Max传入的是str中元素的所有可能的种类数
{
    int m=Max,*x=a,*y=b;//m和c数组是用来配合做统计的，x可以看成是串的特性数组，一开始是串的原本内容，经过下面的大循环之后存
                        //的是每个后缀的排名，y用来存第一关健字的排名结果
    for(int i=0;i<m;i++)
        c[i]=0;
    for(int i=0;i<top;i++)
        c[x[i]=str[i]]++;
    for(int i=1;i<m;i++)
        c[i]+=c[i-1];
    for(int i=top-1;i>=0;i--)
        sa[--c[x[i]]]=i;
    /************************************************
    上面是对单个元素进行排序，相当于初始化sa数组
    ************************************************/
    for(int k=1;;k=k*2)//k代表第一关健字和第二关键字的间隔
    {
        int p=0;
        for(int i=top-k;i<top;i++)//当i大于等于top-k时，不可能有第二关键字
            y[p++]=i;//y数组记录的是第一关键字的排名
        for(int i=0;i<top;i++)
            if(sa[i]>=k)//当sa[i]>=k时，说明该序号可以作为第二关键字，又因为sa是按字典序排过序的（不是最终排序，有一些在之前的比较中还无法
                        //区分次序，他们的排名是随机的，但是可以区分的都已排序），所以可以直接按顺序存进来
                y[p++]=sa[i]-k;//减掉k就可以获得第一关键字的序号
        /***********************************
        上面是对第二关键字进行排序
        ***********************************/
        for(int i=0;i<m;i++)
            c[i]=0;
        for(int i=0;i<top;i++)
            c[x[y[i]]]++;
        for(int i=1;i<m;i++)
            c[i]+=c[i-1];
        for(int i=top-1;i>=0;i--)
            sa[--c[x[y[i]]]]=y[i];
        /************************************
        上面是对第一关键字进行排序
        ************************************/
        int *t=x;
        x=y;
        y=t;
        /************************************
        交换指针，这时y是特性数组
        ************************************/
        p=1;
        x[sa[0]]=0;
        for(int i=1;i<top;i++)
            x[sa[i]]=(y[sa[i-1]]==y[sa[i]]&&y[sa[i]+k]==y[sa[i-1]+k])?p-1:p++;
        /***********************************************************************
        更新特性数组，y看成是上一阶段的特性数组，用来推新的特性数组赋给x，这个
        代码的意思就是，sa数组是一个不完善的排序结果（有些暂时无法区分先后的只
        能把他们的次序放在一起，然后随机安排他们的次序），利用其大致的有序性，
        再结合上一阶段的特性数组来推特性数组，如果对于上一个阶段的特性数组（y）
        ，第一关键字和第二关键字和前一个相同，就说明还是不能区分，那么他们的特
        性值应该相等，否则就应该比前一个大1，表示排名在其之后
        ***********************************************************************/
        if(p>=top)//若p>=top，表示所有后缀的次序都可以确定了，那么就没必要在循环了，因为sa数组不会再变了
            break;
        m=p;//更新统计量
    }
    return ;
}
void getHeight(int top)
{
    for(int i=1;i<=top;i++)
        Rank[sa[i]]=i;//获得Rank数组
    int t=0;
    for(int i=0;i<top;i++)//从左向右遍历后缀，下一个后缀的最长公共前缀长度一定大于等于当前后缀的值减1，原因很简单，详见白书
    {
        if(t>0)
            t--;//这一步是关键
        int j=sa[Rank[i]-1];
        while(str[i+t]==str[j+t])
            t++;
        height[Rank[i]]=t;
    }
    return ;
}

下面上一道模板题，注意main函数

poj3729：http://acm.hust.edu.cn/vjudge/problem/viewProblem.action?id=15681

#include<stdio.h>
#include<algorithm>
using namespace std;
const int MAXN=100005;
int str[MAXN],sa[MAXN],height[MAXN],Rank[MAXN],a[MAXN],b[MAXN],c[MAXN];
int n,m,k;
void getSa(int top,int Max)
{
    int m=Max,*x=a,*y=b;
    for(int i=0;i<m;i++)
        c[i]=0;
    for(int i=0;i<top;i++)
        c[x[i]=str[i]]++;
    for(int i=1;i<m;i++)
        c[i]+=c[i-1];
    for(int i=top-1;i>=0;i--)
        sa[--c[x[i]]]=i;
    for(int k=1;;k=k*2)
    {
        int p=0;
        for(int i=top-k;i<top;i++)
            y[p++]=i;
        for(int i=0;i<top;i++)
            if(sa[i]>=k)
                y[p++]=sa[i]-k;
        for(int i=0;i<m;i++)
            c[i]=0;
        for(int i=0;i<top;i++)
            c[x[y[i]]]++;
        for(int i=1;i<m;i++)
            c[i]+=c[i-1];
        for(int i=top-1;i>=0;i--)
            sa[--c[x[y[i]]]]=y[i];
        int *t=x;
        x=y;
        y=t;
        p=1;
        x[sa[0]]=0;
        for(int i=1;i<top;i++)
            x[sa[i]]=(y[sa[i-1]]==y[sa[i]]&&y[sa[i]+k]==y[sa[i-1]+k])?p-1:p++;
        if(p>=top)
            break;
        m=p;
    }
    return ;
}
void getHeight(int top)
{
    for(int i=1;i<=top;i++)
        Rank[sa[i]]=i;
    int t=0;
    for(int i=0;i<top;i++)
    {
        if(t>0)
            t--;
        int j=sa[Rank[i]-1];
        while(str[i+t]==str[j+t])
        {
            t++;
        }
        height[Rank[i]]=t;
    }
    return ;
}
__int64 getAns(int mid,int top)
{
    __int64 t1=0,t2=0;
    __int64 ans=0;
    for (int i=1;i<=top;i++)
    {
        if (height[i]<mid)
        {
            if (t2>0) ans+=t1;
            t1=t2=0;
            if (sa[i]<n) t1++;
            if (sa[i]>n) t2++;
        }
        else
        {
            if (sa[i]<n) t1++;
            if (sa[i]>n) t2++;
        }
    }
    return ans;
}
int main()
{
    while(scanf("%d%d%d",&n,&m,&k)!=EOF)
    {
        int top=0;
        for(int i=0;i<n;i++)
        {
            scanf("%d",&str[top]);
            str[top++]++;//自增1
        }
        str[top++]=10002;//用极限值加1来连接
        for(int i=0;i<m;i++)
        {
            scanf("%d",&str[top]);
            str[top++]++;//自增1
        }
        str[top]=0;//末尾补零
        getSa(top+1,10003);
        getHeight(top);
        __int64 ans=getAns(k,top)-getAns(k+1,top);//大于等于k的个数减去大于等于k+1的个数就是结果
        printf("%I64d\n",ans);
    }
    return 0;
}