Common SubstringsDescription
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 0Sample 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 ",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 }
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 ",ans); 129 } 130 return 0; 131 }
记得long long~
2016-09-16 11:47:27
更新