OI补完计划——Day1后缀数组

Day1

Part 1,Suffix Array

{昨天睡得迟，今天十点才起床，感觉有点想要感冒，千万不要啊啊啊！拼命喝水中。。。。}

今天是计划的第一天，主要是来复习一下之前学过的算法，并总结自己的见解。

Ps：这不是算法讲稿，想学算法的可以看集训队的论文。

所有算法中我自认理解最深的应该是后缀数组了吧，虽然我不会DC3等线性构造方法，但这些算法难于编写和调试，且算法常数较大，实测的速度优势并不明显。（可能说法有些偏颇，但这不是重点）因此，还是以倍增法来主要介绍其构造方法。（深谙其道的可以直接跳到分割线下）

所谓倍增，即通过已解决的范围较小的问题，使用适当的方法，使其结果可以用来解决范围较大的问题，最终得到问题的最终解。因其解题范围往往是成倍增加的，因此得名为倍增法。其应用广泛，RMQ中的ST算法，后缀数组的构造，及市场倍增学等等。。。（扯远了吧。。。）

而后缀数组用到的倍增思想是既然我难以比较所有后缀的大小，那后缀首字母我总能比较吧，这将给我们一个初始序列，这是所有后缀按照首字母排好的顺序。说是首字母但其实是后缀的前缀（以下简称前缀），我们的任务就是不断以后缀的前缀为后缀排序，当前缀长度超过字符串长度时，不就是后缀的顺序吗。而我们又注意到我们可以通过将两个前缀连成一个更长的前缀来扩大我们的解题范围。但问题来了，如同分治思想的能分也能并，我们不仅要扩大问题的规模，同时也要扩大解题的规模。

而我们发现对于字符串S,T(长度相同且已知大小关系),A,B(已知大小关系)，那么我们既可以判断S+A，T+B的大小关系。方法很简单，若S<>T那么高下立判，若S=T则A,B确定了两字符串的关系。

这样回到刚才的问题中，我们可以在O(1)的时间内，用较短前缀的关系来确定较长前缀的关系。这样便使解题的规模不断变大，直至完成整个后缀数组的构造。

其次是辅助后缀数组的一些信息，以此来将后缀数组用与实际应用中，Height数组记录的是，每个后缀与后缀数组中前一个后缀的最长公共前缀(LCP)的长度，通过一定的方法可以在O(n)时间内完成。方法是对于已经求好的两个后缀的LCP，将第二个后缀的首字母去掉则形成了一个新的后缀，在它之前（后缀数组中的顺序）一定有刚才第一个后缀去掉首字母得到的后缀，不过他们不一定相邻。如果相邻则此后缀的Height即为刚才求得的LCP长度减一。如果不相邻呢，中间插入的一定是让此后缀Height在刚才求得的LCP长度减一基础上增加的后缀。这样我们不断利用已求得的解来省略多余的比较，这样就完成了Height数组的求法。

还有一个便是任意两个后缀的LCP(除模式匹配外，多数关于子串的统计问题是用不到他的，而模式匹配应用AC自动机会更好一点吧)求法是这两个后缀之间Height数组的最小值（比较明显吧。。。裸的RMQ）

以上均是对构造方法的学习笔记，毕竟是别人的东西，下面我想谈谈自己在学习并应用它时的一些见解。

首先是模式匹配，通过二分法在后缀数组中找最和模式串相匹配的后缀，即后缀与模式串的公共前缀最长，若长度等于模式串的长度即匹配成功，而我们若关心模式串在母串中出现的次数则可用二分法确定上下界。同样在二分查找过程中有许多比较是不必要的，因此可以用上文提到的任意两个后缀的LCP来优化匹配过程。

------------------------------分呀分呀分割线-------------------------------------

{重点终于来了。。我郁闷了。。我的表达能力太差上面啰啰嗦嗦写了一大堆，才进入今天的正题。。。}

本文的重点在于后缀数组运用与字符串的统计问题。详细分类可见http://wenku.baidu.com/view/228caa45b307e87101f696a8.html。论文中提到了两种统计方法，一种为二分答案后为用Height数组为后缀分组，另一种是用单调栈来实现统计。

因为我之前看过一些后缀树的内容，我认为这两种方法实际上都是对后缀树遍历的模拟。关于字符串的统计，后缀树是十分强大的数据结构。通过遍历它几乎可以得到我们想要的所有信息。

先说第二种方法，它是对后缀树的深搜遍历，它利用了栈来模拟实现。首先我们可以发现，后缀数组正是后缀树按深度优先遍历顺序得到的叶子节点，而后缀树中相邻两个叶子节点的最近公共祖先(LCA)正是后缀数组中的Height数组。这样建立一一对应关系后。我们便可以模拟遍历后缀树的过程。先将第一个后缀入栈这相当于后缀树遍历到第一个叶子节点。然后重点在于模拟遍历后缀树时回溯的过程，我们通过Height数组控制当前栈顶节点使其为下一未遍历后缀的父亲结点，即每次都为下一次的访问做好准备（具体实现见文末程序），这样直至模拟遍历完整个后缀树，而我们在遍历的同时就可以获得后缀树中相应节点的信息，其中储存了以其为根的子树的全部信息，如后缀的最大最小起始位置，各串后缀出现的次数等等。这些信息均可由子节点的信息合并而来。这样O(n)时间内即可得到我们想要的所有信息，从而快速的解决问题。

但是问题又出现了，后缀树的高效是建立在其可以快速合并节点信息的基础上的。如果合并节点信息的时间复杂度超过了O(1)，则遍历的时间复杂度就会超出O(n)，使问题无法高效解决。是否存在这样的难以合并的信息呢？答案是肯定的，其实这里的合并节点信息和线段树中的有些相像，线段树难以维护的是该线段内有多少种不同的颜色，对应于后缀树则是子树中有多少种属于不同字符串的后缀。当字符串较少时，我们可以考虑用集合或位运算来实现（此时仍为O(1)），但当字符串十分多(1000或更多)时，合并的时间便大大的增加。我们当然可以用许多数来表示集合（同样使用位运算），但这样速度虽被优化（只是优化常数），但编程复杂度却进一步增加，有可能得不偿失。这不得不启发我们从另一个角度思考问题，既然我们难以合并子节点的信息，何不直接用该节点所包含的后缀来直接更新本节点的信息。这就是刚才提到的第一种方法。

第一种方法是对后缀树的分层遍历，有些类似于广搜。但如果对后缀树中的每一个节点都用其子树中叶子节点直接更新，时间复杂度将会到达O(Hn)。（H为整个后缀树的高度）同样会面临超时的窘境。而我们遇到的此类题目则多为求解最值问题等最优化问题。那么我们可以通过当前层的信息来判定最优解是在当前层的上方还是下方，这样通过二分快速的定位我们所需要寻找的节点来解决问题。它的好处在于，可以在O(n)时间内计算出当前层节点的信息，若符合则看下面层的节点是否符合，不符合则向上，直至解决我们的问题。时间复杂度共O(nlogH)。而如果时间较为宽松的话，我们可以用这个方法在未经路径压缩的后缀树上来搞，编程复杂度会进一步下降，时间复杂度为O(nlogn)，也就是论文中提到的用二分答案的方法给Height数组分组。

当然这还没有完美的解决所有字符串统计的问题（没有最坑爹，只有更坑爹），例如题目中要求求出N个字符串，长度不小于k的不同公共子串的个数，两个公共子串相同当且仅当其在每一个字符串的出现的位置均相同，即当N=2时

AABBAABBA 和 ABBAABBB 不同的公共子串ABB共有4个。

（颜色不同的ABB可自由搭配。）可能表述不太清楚，详细可参见POJ3415，是本问题的简化版只有两个字符串。

那么。。。。。。显而易见，第一种方法不可用（没有最优化，因此不具有可二分性），第二种方法则难于合并子节点的信息，导致问题陷入窘境。

其实不必这么伤心，若题目中的N较大则每个字符串的长度都不会过大，因为如果两者均很大那么会导致将所有的字符串连在一起后长度过长，O(n)都难以实现后缀数组的构造。此时，后缀树的高度会比较小，那么用O(H*len)(len为后缀数组的长度)的方法还是可以一试的。如果N较小，则单调栈的方法也是值得一试的。但这终究是建立在题目是使用后缀数组来解的条件下，或许还有更好的方法，但这超出了我现在已有的水平了，不再做进一步的讨论。

Ps：表述能力奇差，大家将就看一下吧。

【欢迎大家指出文中的不足，我们再作进一步的讨论，如有错误我会及时修改】

【Snow Dancer原创，转载请注明，谢谢】

{POJ3415,文末问题的简化版}
var
  s,sa,sap,num:array[0..200010] of longint;
  r:array[0..1,0..200010] of longint;
  stack:array[0..200010] of record a,b,h:int64;end;
  height:array[0..200010] of int64;
  n,i,j,k,l,size,x,y,di,top,limit:longint;
  ans:int64;
  c:char;
function max(x,y:int64):int64;
  begin if x<y then exit(y) else exit(x); end;
procedure BuildSA;//构造后缀数组
  begin
    size:=255;
    for i:=1 to size do num[i]:=0;
    for i:=1 to n do inc(num[s[i]]);
    for i:=1 to size do inc(num[i],num[i-1]);
    r[1]:=s;
    for i:=n downto 1 do begin
      sa[num[s[i]]]:=i;dec(num[s[i]]);
    end;
    k:=1; x:=1;y:=0;
    repeat
      for l:=1 to k do sap[l]:=n-l+1;
      for i:=1 to n do
        if sa[i]>k then begin
          inc(l);sap[l]:=sa[i]-k;
        end;
      for i:=1 to size do num[i]:=0;
      for i:=1 to n do inc(num[r[x,sap[i]]]);
      for i:=2 to size do inc(num[i],num[i-1]);
      for i:=n downto k do begin
        sa[num[r[x,sap[i]]]]:=sap[i];dec(num[r[x,sap[i]]]);
      end;
      size:=0;
      for i:=1 to n do begin
        if not ((r[x,sa[i]]=r[x,sa[i-1]])and(r[x,sa[i]+k]=r[x,sa[i-1]+k])) then inc(size);
        r[y,sa[i]]:=size;
      end;
      x:=y;y:=y xor 1;k:=k<<1;
    until (k>=n)or(size=n);
  end;
procedure calheight;//计算Height数组
  begin
    k:=0;
    for i:=1 to n do begin
      if k>0 then dec(k);j:=sa[r[x,i]-1];
      while s[i+k]=s[j+k] do inc(k);
      height[r[x,i]]:=k;
    end;
  end;
procedure Traversal;    //深搜遍历过程
  procedure push(k:longint);
    begin
      if height[k+1]>stack[top].h then begin
        inc(top);
        stack[top].h:=height[k+1];
        stack[top].a:=0;
        stack[top].b:=0;
      end;
      if sa[k]<di then
        inc(stack[top].a)     //栈中记录的信息可根据题目变化
      else
        inc(stack[top].b);
    end;
  procedure pop(k:longint);
    begin
      if stack[top].h<limit then begin dec(top);exit;end;
      if height[k+1]<=stack[top-1].h then begin
        inc(ans,stack[top].a*stack[top].b*(stack[top].h-max(limit-1,stack[top-1].h)));
        inc(stack[top-1].a,stack[top].a);
        inc(stack[top-1].b,stack[top].b);
        dec(top);
      end else begin
        inc(ans,stack[top].a*stack[top].b*(stack[top].h-max(height[k+1],limit-1)));
        stack[top].h:=height[k+1];
      end;
    end;
  begin
    for i:=1 to n do begin
      push(i);
      while stack[top].h>height[i+1] do pop(i);
    end;
  end;
begin
  readln(limit);
  while limit<>0 do begin
    n:=0; ans:=0; top:=0;
    fillchar(height,sizeof(height),0);
    fillchar(num,sizeof(num),0);
    fillchar(sa,sizeof(sa),0);
    fillchar(sap,sizeof(sap),0);
    fillchar(s,sizeof(s),0);
    fillchar(r,sizeof(r),0);
    while not eoln do begin
      read(c);inc(n);s[n]:=ord(c);
    end;readln;
    inc(n);s[n]:=124;di:=n;
    while not eoln do begin
      read(c);inc(n);s[n]:=ord(c);
    end;readln;
    BuildSA;
    calheight;
    Traversal;
    writeln(ans);
    readln(limit);
  end;
end.