您的位置:首页 > 其它

后缀数组(长度不小于k的公共子串的个数)

2014-07-04 18:56 656 查看
把两个串拼接起来,中间用神奇地‘$'字符隔开。然后用倍增算法求sa数组、rank数组、height数组,然后利用height数组统计第一个串称A串的所有后缀与第二个串B串的所有后缀所共有的长度大于k的子串。
上面的逻辑过程容易理解,用aa和aa就能推出为什么答案是5,子串a有2个,aa有1个,2*2 + 1*1 = 5.问题的关键是要找出各个子串总匹配数(子串相等称为匹配)。height数组表示的是和其他串所有最长公共前缀中的最长公共前缀,也是和前一个串的公共前缀,用这个数组就可以求出总匹配数。height[i]-k+1代表和前一个子串相比符合条件的匹配数,如果height[i]大于它属于组内(根据height[i]是否大于k分组,height数组是波浪形,很容易分组)前面的每一个height[j]的话,就可以一次一次和前面的height[j]比较如果是大于,就增加height[i]-k+1,这就代表从sa[i]开始的后缀与sa[j]开始的后缀有height[i]-k+1个相等的串,一旦小于,就应该做些处理,把多出来的那部分去掉,因为这部分以后再也不会有相等的串了,没什么价值了。但这样,复杂度是O(n^2),可用单调栈优化,具体实现见代码。

A串对B串,B串对A串的统计过程是一样的,分析下A串对B串然后以此类推到B串对A串即可。先说下统计公共子串的原理,其实是利用A串的后缀与排名比它高的B串的后缀的最长公共前缀l进行统计,如果l大于k,就增加l-k+1个公共子串。
然后要说到这个最长公共前缀,子串b与当前串(排名比他低的A串后缀)的最长公共前缀是它的高一位排名到当前后缀之间最小的那个height,我们要统计前面所有B串的子串和当前A串后缀所共有的公共子串只需要用到几个height值就好,而这些height值是递增的。为什么是递增的,比如0,1,2,3,1,2,3,那么统计到第一个1的时候,0的和它的最长公共前缀为1,统计到2的时候前两个的最长公共前缀是1,2,统计到3的时候是1,2,3。统计到1的时候则有一个很大的变动,因为前4个的最长公共前缀都可以是这个1,所有只用到这个1就好...
正因为上面一段的分析,我们可以用一个栈和一个值total来动态维护前面所有后缀和当前串的最长公共前缀长度。栈里存的是height值,total其实就是栈里递增的所有height值和当前串所共有的大于k的公共子串。当我们的某个更小的height值进栈时,需要维护这个total值,因为栈顶的几个height相对大些,表示的公共子串多些,当遇到更小的height值,最长公共前缀长度应该是更小的那个,所以要地减去那一部分多出来的,然后用个num[to]表示栈顶的这个height代表几个大灯与它的height。

#include<iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <string>
#include <queue>
#include <algorithm>
#include <map>
#include <iomanip>
#define INF 99999999
typedef long long LL;
using namespace std;

const int MAX=2*100000+10;
int *rank,r[MAX],sa[MAX],height[MAX];
int wa[MAX],wb[MAX],wm[MAX];
char s[MAX];

bool cmp(int *r,int a,int b,int l)
{
    return r[a] == r[b] && r[a+l] == r[b+l];
}

void makesa(int *r,int *sa,int n,int m)
{
    int *x=wa,*y=wb,*t;
    for(int i=0; i<m; ++i)wm[i]=0;
    for(int i=0; i<n; ++i)wm[x[i]=r[i]]++;
    for(int i=1; i<m; ++i)wm[i]+=wm[i-1];
    for(int i=n-1; i>=0; --i)sa[--wm[x[i]]]=i;
    for(int i=0,j=1,p=0; p<n; j=j*2,m=p)
    {
        for(p=0,i=n-j; i<n; ++i)y[p++]=i;
        for(i=0; i<n; ++i)if(sa[i]>=j)y[p++]=sa[i]-j;
        for(i=0; i<m; ++i)wm[i]=0;
        for(i=0; i<n; ++i)wm[x[y[i]]]++;
        for(i=1; i<m; ++i)wm[i]+=wm[i-1];
        for(i=n-1; i>=0; --i)sa[--wm[x[y[i]]]]=y[i];
        for(t=x,x=y,y=t,i=p=1,x[sa[0]]=0; i<n; ++i)
        {
            x[sa[i]]=cmp(y,sa[i],sa[i-1],j)?p-1:p++;
        }
    }
    rank=x;
}

void calheight(int *r,int *sa,int n)
{
    for(int i=0,j=0,k=0; i<n; height[rank[i++]]=k)
    {
        for(k?--k:0,j=sa[rank[i]-1]; r[i+k] == r[j+k]; ++k);
    }
}

LL calculate(int n,int len,int k)
{
    int *mark=wb,*ans=wm,Top=0;//num[1],num[2]分别表示字符串A,B,suffix(0~i-1)和suffix(i)的最长公共子串>=k的总个数
    LL sum=0,num[3]= {0};
    for(int i=1; i<=n; ++i)
    {
        if(height[i]<k)
        {
            Top=num[1]=num[2]=0;
        }
        else
        {
            for(int size=Top; size && ans[size]>height[i]-k+1; --size) //维护单调栈,ans记录的是suffix(j)和suffix(i-1)>=k的最长公共子串的个数,个数越多表示height[j]越大
            {
                num[mark[size]]+=height[i]-k+1-ans[size];//suffix(j)和suffix(i)>=k的最长公共子串只能是长度为k~height[i],所以需要减去(ans[size]-(height[i]-k+1))
                ans[size]=height[i]-k+1;//更新个数(更新单调栈,使栈里面元素非递减)
            }
            ans[++Top]=height[i]-k+1;
            if(sa[i-1]<len)mark[Top]=1;//由于num新增加的结果是suffix(i-1)和suffix(i)的结果,所以是判断sa[i-1]
            if(sa[i-1]>len)mark[Top]=2;
            num[mark[Top]]+=height[i]-k+1;//增加由suffix(i-1)和suffix(i)产生的结果
            if(sa[i]<len)sum+=num[2];//表示和suffix(i)产生的结果新增加B串的suffix(0~i-1)和suffix(i)>=k的个数
            if(sa[i]>len)sum+=num[1];//表示和suffix(i)产生的结果新增加A串的suffix(0~i-1)和suffix(i)>=k的个数
        }
    }
    return sum;
}

int main()
{
    int k,n,len;
    while(~scanf("%d",&k),k)
    {
        scanf("%s",s);
        for(n=0; s
 != '\0'; ++n)r
=s
;
        r[len=n]='#';
        scanf("%s",s+n+1);
        for(++n; s
 != '\0'; ++n)r
=s
;
        r
=0;
        makesa(r,sa,n+1,256);
        calheight(r,sa,n);
        cout<<calculate(n,len,k)<<endl;
    }
    return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: 
相关文章推荐