【POJ3415】 Common Substrings（后缀数组|SAM）

Common Substrings

Description

A substring of a string T is defined as:

T(i, k)=TiTi+1...Ti+k-1, 1≤i≤i+k-1≤|T|.

Given two strings A, B and one integer K, we define S, a set of triples (i, j, k):

S = {(i, j, k) | k≥K, A(i, k)=B(j, k)}.

You are to give the value of |S| for specific A, B and K.

Input

The input file contains several blocks of data. For each block, the first line contains one integer K, followed by two lines containing strings A and B, respectively. The input file is ended by K=0.

1 ≤ |A|, |B| ≤ 105
1 ≤ K ≤ min{|A|, |B|}
Characters of A and B are all Latin letters.

Output

For each case, output an integer |S|.

Sample Input

2
aababaa
abaabaa
1
xx
xx
0

Sample Output

22
5


【题意】


　　给两个长度不超过100000的字符串A和B, 求三元组(i, j, k)的数量, 即满足A串从i开始的后缀与B串从j开始的后缀有长度为k的公共前缀, 还要求k不小于某个给你的数K.


 


【分析】


　　如果i后缀与j后缀的LCP长度为L, 在L不小于K的情况下, 它对答案的贡献为L - K + 1. 


　　所以问题就是快速求出两个串的后缀的LCP。


　　就要用到后缀数组，把两个串拼成一个串，中间加一个特殊字符，然后做后缀数组。


　　求出height数组后根据k分组（min大于等于k的分在一组），同一组的LCP一定大于等于k。（height数组的性质啦，这个很基本）


　　


　　接下来就是求出sum｛L-K+1｝，如果直接暴力的话for2遍加一个RMQ也很慢。


　　然后看到大神们说用什么单调栈，哦~~ 单调的优美性质！！！ 


　　


　　可以分成两步求，先求B串在下，A串在上的ans。再反过来求一遍。


　　过程好像比较难说清楚。


　　反正...对于B串上面的一堆A串，因为LCP是求min，所以LCP必定是单调递增，所以像下图那样做：


　　


　　开了三个数组，s表示其单位贡献值，cnt表示有多少个A串是这个贡献值，h为从当前位置到底部的cnt*s的和（h便于计算的）。


　　然后就更新和替换。遇到一个B串就把h[tp]累加到ans中即可。　


　　反过来做也是一样的。


　　单调！！单调！！单调哦！！


　　好像很厉害！！


　　主要看代码~~


 


　　还有，这道题是有大写字母的！！！坑的我RE了巨久！！


　　还有记得long long！！


 


代码如下：






  1 #include<cstdio>
  2 #include<cstdlib>
  3 #include<cstring>
  4 #include<iostream>
  5 #include<algorithm>
  6 #include<queue>
  7 using namespace std;
  8 #define INF 0xfffffff
  9 #define Maxl 200010
 10 #define LL long long
 11 
 12 int k,la;
 13 char a[Maxl],b[Maxl];
 14 int c[Maxl];
 15 int cl;
 16 
 17 int sa[Maxl],rank[Maxl],Rs[Maxl],wr[Maxl],y[Maxl];
 18 //sa -> 排名第几的是谁
 19 //rank -> i的排名
 20 //Rs数值小于等于i的有多少个
 21 //y -> 第二关键字排名第几的是谁（类似sa）
 22 int height[Maxl];
 23 
 24 void get_sa(int m)
 25 {
 26    memcpy(rank,c,sizeof(rank));
 27     for(int i=0;i<=m;i++) Rs[i]=0;
 28     for(int i=1;i<=cl;i++) Rs[rank[i]]++;
 29     for(int i=1;i<=m;i++) Rs[i]+=Rs[i-1];
 30     for(int i=cl;i>=1;i--) sa[Rs[rank[i]]--]=i;
 31     
 32     int ln=1,p=0;
 33     while(p<cl)
 34     {
 35        int k=0;
 36         for(int i=cl-ln+1;i<=cl;i++) y[++k]=i;
 37         for(int i=1;i<=cl;i++) if(sa[i]>ln) y[++k]=sa[i]-ln;
 38         for(int i=1;i<=cl;i++) wr[i]=rank[y[i]];
 39         
 40         for(int i=0;i<=m;i++) Rs[i]=0;
 41         for(int i=1;i<=cl;i++) Rs[wr[i]]++;
 42         for(int i=1;i<=m;i++) Rs[i]+=Rs[i-1];
 43         for(int i=cl;i>=1;i--) sa[Rs[wr[i]]--]=y[i];
 44         
 45         for(int i=1;i<=cl;i++) wr[i]=rank[i];
 46         for(int i=cl+1;i<=cl+ln;i++) wr[i]=0;
 47         p=1;rank[sa[1]]=1;
 48         for(int i=2;i<=cl;i++)
 49         {
 50             if(wr[sa[i]]!=wr[sa[i-1]]||wr[sa[i]+ln]!=wr[sa[i-1]+ln]) p++;
 51             rank[sa[i]]=p;
 52         }
 53         m=p,ln*=2;
 54     }
 55     sa[0]=rank[0]=0;
 56 }
 57 
 58 void get_he()
 59 {
 60     int kk=0;
 61     for(int i=1;i<=cl;i++)
 62     {
 63         // if(rank[i]==1) break;
 64         int j=sa[rank[i]-1];
 65         if(kk) kk--;
 66         while(c[i+kk]==c[j+kk]&&i+kk<=cl&&j+kk<=cl) kk++;
 67         height[rank[i]]=kk;
 68     }
 69 }
 70 
 71 LL h[Maxl],cnt[Maxl],s[Maxl];
 72 int tp;
 73 void ffind()
 74 {
 75     LL ans=0;tp=0;
 76     s[0]=INF;
 77     //******
 78     for(int i=1;i<=cl-1;i++)
 79     {
 80         if(sa[i]>=la+1) //B串
 81          ans+=h[tp];
 82         if(height[i+1]-k+1<s[tp])//替换
 83         {
 84             LL sum=0;
 85             while(height[i+1]-k+1<=s[tp]&&tp) sum+=cnt[tp--];
 86             s[++tp]=(height[i+1]-k+1),cnt[tp]=sum,h[tp]=h[tp-1]+cnt[tp]*s[tp];
 87         }
 88         else s[++tp]=height[i+1]-k+1,cnt[tp]=0,h[tp]=h[tp-1];
 89         if(sa[i]<=la) cnt[tp]++,h[tp]+=s[tp];//A串
 90         if(height[i+1]<k)
 91         {
 92             tp=0;s[0]=INF;
 93         }
 94     }tp=0;s[tp]=INF;
 95     for(int i=1;i<=cl-1;i++)
 96     {
 97         if(sa[i]<=la) //A串
 98          ans+=h[tp];
 99         if(height[i+1]-k+1<s[tp])//替换
100         {
101             LL sum=0;
102             while(height[i+1]-k+1<=s[tp]&&tp) sum+=cnt[tp--];
103             s[++tp]=(height[i+1]-k+1),cnt[tp]=sum,h[tp]=h[tp-1]+cnt[tp]*s[tp];
104         }
105         else s[++tp]=height[i+1]-k+1,cnt[tp]=0,h[tp]=h[tp-1];
106         if(sa[i]>=la+1) cnt[tp]++,h[tp]+=s[tp];//B串
107         if(height[i+1]<k)
108         {
109             tp=0;s[tp]=INF;
110         }
111     }
112     printf("%I64d\n",ans);
113 }
114 
115 void init()
116 {
117     scanf("%s%s",a,b);
118     int l=strlen(a);cl=0;
119     la=l;
120     for(int i=0;i<l;i++)
121     {
122         if(a[i]>='a'&&a[i]<='z') c[++cl]=a[i]-'a'+1;
123         else c[++cl]=a[i]-'A'+27;
124     }
125     c[++cl]=55;
126     l=strlen(b);
127     for(int i=0;i<l;i++)
128     {
129         if(b[i]>='a'&&b[i]<='z') c[++cl]=b[i]-'a'+1;
130         else c[++cl]=b[i]-'A'+27;
131     }
132 }
133 
134 int main()
135 {
136     while(1)
137     {
138         scanf("%d",&k);
139         if(k==0) break;
140         init();
141         get_sa(55);
142         get_he();
143         ffind();
144     }
145     return 0;
146 }



[POJ3415]


 


2016-07-17 15:00:50




 


update：

哈哈学了SAM的我回来更新这道题啦！！

然而很搞笑的是...我打的SAM其实跟后缀数组差不多长！！

而且蒟蒻表示调试了很久！狗带！！SAM太难调了！！



好吧我一开始想不出来..是因为·不太懂那个SAM上某个点表示的串究竟是什么~~

其实是一个连续的区间，记为min(x)~max(x)（这是长度）,

其中min(x)=step[pre[x]]+1,max(x)=step[x]...

为什么呢..其实SAM就是省掉了重复子串吧，我觉得~~


    比如说上面这个，parent最大串一定是他的后缀来的，后面那个b表示的后缀其实有aab、ab、b、xaab....
但是aab、ab、b明显他的parent也有，就不会用后面那个来表示，
所以最短的都有step[pre[x]]+1那么长而且到时如果还要表示aab-----的串的话，其实信息也是保存到parent那里了。（大概，意会一下~）
   
 关于这道题，我们处理的目标是前缀的后缀~（上面用后缀数组的方法处理的是后缀的前缀，都叫后缀xx其实怎么我觉得是反过来的呢~）
　　
用A串建机，B串跑。

跑法有点类似AC自动机上的跑法，跑出来的长度啊都是让后缀匹配最长的！！

我们知道如果有长度大于k的子串s和ss，如过ss是s的子串，而且s合法（就是在A串出现），那么ss也合法，而且ss的答案一定大于等于s的答案。

我一开始就是不知道怎么弄ss多余部分的答案。
　　
    很容易想到就是用rt[x]*(length-k+1),但是直接减的话，你算一些小的串，其实x并不能代表这个串，就是说你现在扫到的位置不一定是这个小串的首位置，那么直接用right数组就会有问题，就会漏解。
　　
    上面说了，x表示的串为min(x)~max（x），我们就在x的位置先算他表示的串，然后在parent那里打标记，用parent的right数组来算，就可以避免这种问题。





　　具体怎么做看GDXB讲解：


先对第一个串构造 SAM,逆拓扑序求出right存入r[]。
某个节点的right集合表示Min(x)~Max(x)这一段子串都出现了r[x]次！
用第二个串对 SAM 做 LCS，当前节点x LCS>=K时，ans+=ans+=r[x]*(len-maxx(k,step[pre[x]]+1)+1);（当前匹配的最长串的子串数）。如果step[pre[x]]>=k,cnt[pre[x]]++;
拓扑序逆序统计不是最长公共子串的状态但是被子串包含的个数


，ans+=cnt[p]*(step[p]- max(K,Min(p)+1)*r[p]，同时维护cnt:cnt[pre[p]]+=cnt[p]。


 







 


 


说一说right数组怎么求，因为我这个一开始都打错了！！


 


right[x]表示x这个点表示的串在原串出现多少次。


 


我们把主链上的right全部清成1，然后按拓扑序把i的答案放到pre[i]上即可。


 


代码如下：(SAM)






  1 #include<cstdio>
  2 #include<cstdlib>
  3 #include<cstring>
  4 #include<iostream>
  5 #include<algorithm>
  6 #include<queue>
  7 #include<cmath>
  8 using namespace std;
  9 #define Maxn 100010
 10 #define LL long long
 11 
 12 char s[Maxn],ss[Maxn];
 13 int len,ml;
 14 LL lim;
 15 
 16 LL ans;
 17 LL mymax(LL x,LL y) {return x>y?x:y;}
 18 
 19 struct node
 20 {
 21     int tot,pre[2*Maxn],son[2*Maxn][60];
 22     LL step[2*Maxn],rt[2*Maxn],cnt[2*Maxn];
 23     int last;
 24     void extend(int k)
 25     {
 26         int np=++tot,p=last;
 27         step[np]=step[last]+1;
 28         rt[np]=1;
 29         while(p&&!son[p][k])
 30         {
 31             son[p][k]=np;
 32             p=pre[p];
 33         }
 34         if(!p) pre[np]=1;
 35         else
 36         {
 37             int q=son[p][k];
 38             if(step[q]==step[p]+1) pre[np]=q;
 39             else
 40             {
 41                 int nq=++tot;
 42                 step[nq]=step[p]+1;
 43                 memcpy(son[nq],son[q],sizeof(son[nq]));
 44                 pre[nq]=pre[q];
 45                 pre[q]=pre[np]=nq;
 46                 while(p&&son[p][k]==q)
 47                 {
 48                     son[p][k]=nq;
 49                     p=pre[p];
 50                 }
 51             }
 52         }
 53         last=np;
 54     }
 55     void init()
 56     {
 57         memset(pre,0,sizeof(pre));
 58         memset(son,0,sizeof(son));
 59         memset(rt,0,sizeof(rt));
 60         memset(cnt,0,sizeof(cnt));
 61     }
 62     void build()
 63     {
 64         step[1]=0;tot=1;
 65         last=1;
 66         for(int i=1;i<=len;i++)
 67         {
 68             int k=s[i]-'a'+1;
 69             if(k<0) k=s[i]-'A'+27;
 70             extend(k);
 71         }
 72     }
 73     int aq[2*Maxn],tp[2*Maxn],Rs[2*Maxn];
 74     void get_tp()
 75     {
 76         for(int i=0;i<=tot;i++) Rs[i]=0;
 77         for(int i=1;i<=tot;i++) Rs[step[i]]++;
 78         for(int i=1;i<=tot;i++) Rs[i]+=Rs[i-1];
 79         for(int i=1;i<=tot;i++) tp[Rs[step[i]]--]=i;
 80     }
 81     void get_rt()
 82     {
 83         int now=1;
 84         for(int i=tot;i>=1;i--) rt[pre[tp[i]]]+=rt[tp[i]]; 
 85     }
 86     void ffind()
 87     {
 88         get_tp();
 89         get_rt();
 90         int now=1,sp=0;
 91         ans=0;
 92         for(int i=1;i<=ml;i++)
 93         {
 94             int k=ss[i]-'a'+1;
 95             if(k<0) k=ss[i]-'A'+27;
 96             while(now&&!son[now][k]) now=pre[now],sp=step[now];
 97             if(son[now][k]) now=son[now][k],sp++;
 98             else now=1,sp=0;
 99             if(sp>=lim) ans+=(sp-mymax(lim,step[pre[now]]+1)+1)*rt[now];
100             if(lim<step[pre[now]]+1) cnt[pre[now]]++;
101         }
102         for(int i=tot;i>=1;i--) 
103             cnt[pre[tp[i]]]+=cnt[tp[i]];
104         for(int i=1;i<=tot;i++) if(step[i]>=lim)
105             ans+=rt[i]*cnt[i]*(step[i]-mymax(step[pre[i]]+1,lim)+1);
106         
107         for(int i=1;i<=tot;i++) pre[i]=0;
108         for(int i=1;i<=tot;i++)
109          for(int j=1;j<=55;j++) son[i][j]=0;
110         for(int i=1;i<=tot;i++) rt[i]=0;
111         for(int i=1;i<=tot;i++) cnt[i]=0;
112     }
113     
114 }suf;
115 
116 int main()
117 {
118     suf.init();
119     while(1)
120     {
121         scanf("%lld",&lim);
122         if(lim==0) break;
123         scanf("%s%s",s+1,ss+1);
124         len=strlen(s+1);
125         ml=strlen(ss+1);
126         suf.build();
127         suf.ffind();
128         printf("%lld\n",ans);
129     }
130     return 0;
131 }



[POJ3415]


 


记得long long~


 


2016-09-16 11:47:27


更新


 






            

            
posted @ 
2016-07-17 14:43 
konjak魔芋 
阅读(758) 
评论(0) 
编辑 
收藏 
举报