• 【POJ3415】 Common Substrings(后缀数组|SAM)


    Common Substrings

    Description

    A substring of a string T is defined as:

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

    Given two strings AB and one integer K, we define S, a set of triples (ijk):

    S = {(ijk) | kKA(ik)=B(jk)}.

    You are to give the value of |S| for specific AB 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
    ",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
    ",ans);
    129     }
    130     return 0;
    131 }
    [POJ3415]

    记得long long~

    2016-09-16 11:47:27

    更新

     

  • 相关阅读:
    2014025640《嵌入式程序设计》第二周学习总结
    基于Struts2的SpringMVC入门
    2014025640《嵌入式设计》第一周学习总结
    Hadoop综合大作业
    hive基本操作与应用
    用mapreduce 处理气象数据集
    熟悉常用的HBase操作,编写MapReduce作业
    爬虫大作业
    熟悉常用的HDFS操作
    中文词频统计
  • 原文地址:https://www.cnblogs.com/Konjakmoyu/p/5678554.html
Copyright © 2020-2023  润新知