本来打算这篇文章就开始讲SA的应用，可我决定推迟一篇文章，主要原因是：
1. 最近读到的一篇文章中有一个想法我觉得很新颖，这对于希望从事研究类工作的人来说有一定的帮助；
2. 感觉前面的的文章很没有条理，想补充点内容把没说的说完整。

先从SA构造算法的分类开始说起。一般来说，现在有四种思路构造SA，分别如下：
1. Prefix Doubling，这个思路的代表就是MM（倍增）算法；
2. Recursive，代表是DC3算法；
3. Induced Copying，代表是MF算法；
4. 从后缀树构造后缀数组。

其中，Induced Copying目前公认是在实践中最快的方法，但是据我所知目前这类算法实现起来都不太优美，它们大多需要用到Sedgewick和Bentley发明的三路划分字符串快速排序算法，这样代码量较大从而在比赛中不容易推广。
本文要说的算法叫KA（Ko and Aluru），也是Recursive Algorithm中的一种，从论文中找到的数据表明要比DC3算法快且节省内存。我并没有实现该算法，并不是说实现难度很大，而是这个算法开始部分的优美让我陶醉，但后面的部分的笨拙（并非说低效）让我觉得这个算法没能一气呵成，所以我暂时不想去写。
1. 基本算法
首先约定，原串是T。对于T的每个后缀T(i)，把它划分为两个集合：T(i) > T(i+1)的为L集，T(i) < T(i+1)的为S集，其中T(n)属于S。对T按照单个字符基数排序形成数组A，这样T就以字符集E为基础被分成了很多个Bucket，设head(a)表示字母a的Bucket在A中的起始序号，tail(a)表示字母a的Bucket在A中结束序号，如下图所示：

然后，对S集的所有后缀排序，形成数组B（如何做这个后面说明）。
得到了B，从右向左一次扫描。对每一个i，将suf=B[i]放到A中tail[ suf[0] ]的位置，然后tail[ suf[0] ]减少1。这个步骤没什么好说的，实际上就是把B的结果归位到A中。
神奇的是后面。下面先来看一个定理：
定理1.对于两个相同字母开头的后缀suf1和suf2，如果分别来自S集和T集，那么suf1一定最终排在suf2的后面。
其实很容易证明，suf1[0] == suf2[0]，又suf1[0] < suf1[1]，suf2[0] > suf2[1]，所以suf1[1] > suf2[1]，得suf1 > suf2，得证。
于是，利用这个定理，对B数组从左至右扫描，对遇到的每个后缀suf，把该suf放到A数组的tail[ suf[0] ]位置上，并且tail[ suf[0] ]减少1。
后面的思路和DC3类似，要做到O(n)的复杂度，就必须使用S集的排序导出L集的排序，下面首先来看一个过程：
从A数组开始从左至右扫描，对遇到的每个后缀suf，设R[suf]表示这个后缀在T中的起始位置。如果R[suf] – 1是一个T集的后缀，设为suf2，则将其放入head[ suf2[0] ]的位置，并且head[ suf2[0] ]增加1，这个步骤叫做suf1确定suf2。
这个过程一举两得，不但有了L的排序，同时也将S和L归并完成，是KA算法最为让人拍手叫绝的一步（读到此至少我是眼前一亮）。
定理2.以上算法得到的A是T的后缀数组。
要证明它的正确性其实很简单，如下：
1. A的第一个元素肯定是一个S集的后缀。如果是L集的，设为suf1，那么必然存在一个后缀suf2小于suf1。如果suf2仍然是L集的，那么必然有suf3……这个过程不可能永远进行下去，所以必然最小的后缀是来自S集的（是不是感觉和定理1冲突了？）；
2. 对于L中每一个后缀suf，如果要被放到A中，必然是R[suf]+1这个后缀已经被放入A。而R[suf]+1这个后缀除非是S集元素，否则它会要求R[suf]+2的后缀先放入A，直到遇到一个S集的元素，所以，每一个元素都是被A中排在其之前的元素确定的。所以可知算法进行到A[i]时，A[i]必然已赋值；
3. 其次，每一个A[i]的赋值必然是正确的，这一步可以通过数学归纳法来证明。已知第1位A是正确赋值的，设前i位A都能正确赋值，那么第i+1位如果错误，可知必然存在j > i + 1，A[j]应该填充A[i+1]。A[j]必然和A[i+1]的首位字母是一样的（按照算法的流程，不一样是不可能争同一个位置的，对吧？），因此就有R[A[j]]+1的后缀小于R[ A[i+1] ]+1的后缀，那么根据假设和算法的流程，必然会先安排A[j]到位置i+1，因此赋值错误是不可能的。
剩下唯一的难点就是对S进行排序。
方法有二。一是和DC3类似，把所有的S按T中出现的顺序提出来，然后将连续的字符通过排序替换成一个新字符，然后递归形成后缀数组。这样做的难点在于，DC3中只需要将连续的三个字符堪称一个新符号，这里却是不规则的，因此普通的基数排序很难进行，所以算法的作者搞了一个很麻烦的方法，想了解的可以参看算法的原文。
另外一个方法就是I&T算法，是Induced Copying的一种，它直接用前文说的快速字符串排序把S搞定。很直接很暴力，但未必不高效。
最后需要注意的是，其实先将L排好序，再导出S的顺序并形成最终的结果也是可以的，步骤和上面类似，但是放入每个bucket的顺序恰好相反。所以，可以视S和L的数量来定该先排序哪个，这样可以改善效率。
2. 进一步提速
KA算法的导出和归并步骤相比DC3要高效许多，要想进一步提速，就需要减少排序S集的时间。更好的排序思路很难找到，那么就通过给S集瘦身来达到目的吧。
现在，只提取T[i]是S后缀但T[i+1]是L后缀的i形成S*集。将S*按照原有的方法排序，然后依然按原来的方法先吧S*里面的元素分配到A中。
剩下的S集元素可以通过S*的元素导出。方法是，再一次从右向左扫描S*数组，如果R[ B[i] ] – 1是S集的，那么把放到开头字母对应的bucket中最后一个位置上。
这一个改进的证明其实很简单，我不再多说。
其实，这些算法的所有证明其实都和一个性质有关：单调性。参透了这一点就算是对KA算法真正的理解了。

原文下载： KA Suffix Arrary 构造算法

相隔半个月之久，我又看了许多关于SA的论文，但似乎给我留下深刻映像的并不多。我希望寻找简单而又巧妙的思路也未能成功（不一定复杂度很低），也或许是我已经找到了很好的思考角度但是没有很强的分析能力于是依然无获而终。
本文打算详细说一下DC3算法，早在好几个月前便答应要写。而下一篇将会说一个SA很有意思的应用，来自数据压缩领域（猜到是什么了吗？）。
DC3最早见于论文的名字叫skew算法，不久他们伙同Google的一个人一同撰写了另外一篇文章，内容上和前一篇没有很大区别，但是从更本质的角度阐述了skew算法，即difference cover原理，并提出更一般的DC(v)算法，而原来的skew也就正式更名为DC3算法，因此若只想了解DC3，看哪篇文章都一样。
这个算法和前面不同的是它利用了分治的思想，而该思路最早是由Farach提出并用于Suffix Tree的构造。设想，如果我们能分别得到奇数位置的所有后缀和偶数位置的所有后缀形成的SA，然后再进行一次归并操作，是不是就可以O(nlgn)内构造SA呢？
假设偶数位置的所有后缀形成的SA为SA2，奇数位置为SA1，相应的Rank数组为Rank1和Rank2，原字符串是S（注意，S从1开始计数）。现在SA2已经求出，那么对于所有奇数位置的后缀i，利用文章1中的方法构造二元组< S[i], Rank2[i+1]>，显然对于这个二元组的每一个元，大小比较都有良序定义，则直接基数排序之就得到SA1，避免了一次递归调用。
再来考虑SA2该如何求。既然是分治，因此只要能构造出和原问题类似的子问题即可。令S’是S从位置2开始的后缀。先走题一下，设有两个字符串A和B（为方便讨论，两个都是偶数长度），如果有A < B，那将A和B中的每两个字符看成一个字符后（即每个新的字符是一个二元组，大小比较方法是先比较第一维再比较第二维），依然有A’ < B’，即这样的映射具有保序性。同理，将S’中1、2位置压缩成一个字符，3、4位置压缩，这样S’的所有后缀恰好对应S的所有奇数位置压缩后的后缀。由于有保序性，因此求得S’的SA就得到了S的SA2。
由于每轮只递归一次，|S’| = |S| / 2，因此方程是f(n) = f(n/2) + O(n)，即算法的时间复杂度是O(n)，空间复杂度也是O(n)。
以上的算法很完美，但想必大家都能看出算法的描述和数字3没有任何关系，所以以上所述其实并不是DC3算法。细心的读者或许还能看出另外一个破绽，我没有讲归并如何进行。设p1是来自奇数位的后缀，p2是偶数位，则p1[1]和p2[1]比较后，剩下的部分还不能马上确定大小关系（因为又分别来自SA2和SA1）。实际上，的确是有办法可以做到线性时间内合并，但非常的不trivial，不值得推荐。
而真正的DC3和上述算法的区别就是最后一步归并变得非常简单。试想我们已经有了所有起始位置mod 3 != 0的SA12和Rank12，而SA0和Rank0是位置mod 3 = 0的后缀。归并这两者时，设p0，p1，p2分别是mod 3 = 0,1,2的后缀，那么p0和p1归并可以先比较第一个字符，然后就分别剩下起始mod 3 = 1和mod 3 = 2的后缀，而它们的大小关系通过Rank12就直接得到。同理，比较p0和p2时，先比较前两个字符，又剩下两个p2和p1的后缀，大小关系也能直接判别。
而SA12的求法和上述算法中的SA2是同理的，只不过要提取出mod 3 = 1,2的所有后缀，并且是每三个字符压缩成一个新符号。而唯一不同的地方是，SA12中要先连续放置mod 3 = 1的后缀，然后再放置mod 3 = 2的后缀，这样当符号压缩了以后才是原问题的子问题。
那difference cover究竟是什么呢？定义S是一个DC samples集合，对任意i和j属于[0,n)，存在k使得i+k和j+k都属于S，比如上面的SA12。根据此理论，如果每次压缩4个字符，那么也只要选mod 4 = 2, 3的后缀作为子问题排序即可。当然，更详细的数论的结论建议大家看看论文，简单而有启发。
最后，程序的代码在原始论文的后面有，写得很巧妙，基本上没什么可以改的（我试图自己写一个新的实现，结果始终不如论文作者思路清晰），因此我也不贴这里了。

        很久便说起要写后缀数组（简称SA）的相关内容，结果一直以来未曾有动静，原因是我既忙又懒。最近上网搜了些帖子，发现SA这东西早已经深入人心了，自然感觉自己已经开始慢慢的脱离于时代。
其实不用把SA看得有多神秘，它只不过是一个字符串的所有后缀排序后的结果。稍微有点常识的人都可以想到，把所有后缀取出来然后跑一个快排就可以拿到SA。没错，如果你认为O(n^2lgn)的复杂度你完全可以接受的话，显然这是一种最快速的编程方法。
自然，无论是竞赛中还是实际的工程项目中，上诉算法都是难以利用的，它复杂度实在太高了。我们尝试改进这个算法，联想到RK算法利用hash的思想将串匹配问题的复杂度降到了O(n)，那么如果找到一个hash函数F，使得对任意的串S1和S2，其中S1 < S2，有F(S1) < F(S2)，即F具备保序性，那么可以先把这些后缀都映射成整数，然后再快排或基排就可以解决问题。
嗯，这个想法到是好，但为了比原来的算法更优，映射操作的复杂度不能太高，所以这个F怎么构造是个大问题。联想到后缀排序和普通排序相比，有一个很大的不同是后缀之间是有关系的，比如第i个位置开始的后缀第i+1个位置开始的后缀相比就多了S(i)而已。为了利用这个关系，自然想到从后往前取后缀的方法。下面介绍一个实现非常简单的自创O(n^2）算法。
设sa(i)表示排名第i位的后缀的起始位置，rank(i)表示起始位置为i的后缀的排名，显然有sa和rank是互为反函数。假设现在处理到第i个后缀，那么构造二元组P(i) = <S(i), rank(i+1)>。对于之前已经处理过的开始位置为k（k > i）的后缀，有P(i) < P(k) 当且仅当S(i) < S(k) 或者 S(i) == S(k) && rank(i+1) < rank(k+1)，由于i+1在i之前已经处理完毕，自然rank(i+1)此时也是有定义的。
使用类似插入排序的过程，使用上面定义的比较大小的方法从sa数组的后面开始找到suffix(i)的位置。在查找的过程中，对于所有P(k) > P(i)，令rank(k) ++，因为suffix(i)跑到了它的前面去（读者可以想想看这样做为什么不会错？）。
基本原理很简单，当然代码也非常简洁，C描述如下：

for ( i = n – 1; i > -1; –i ) {
for ( j = n – 2 – i; j > -1; –j ) {
if ( buf[i] > buf[ SA[j] ] || ( buf[i] == buf[ SA[j] ] && Rank[i+1] > Rank[ SA[j]+1 ] ) ) break;
SA[j+1] = SA[j];
++Rank[ SA[j] ];
}
SA[j+1] = i;
Rank[i] = j+1;
}
}

[/code]

这个算法最大的有点是编码容易，空间复杂度低，而且常数也低，同时也避免了之前算法的最坏情况。但是，O(n^2)的复杂度也的确很难将其应用到大数据量的场合，充其量可以用于快速建模或者Topcoder这种地方。不过，至少我们打开了思路，觉得SA并非intractable的东西，而且，文中对于已有结果的利用的思路也是可圈可点。
再进一步改进。现在的算法是借用的插入排序的模式，虽然后缀比较的复杂度从O(n)降到了O(1)，但是排序本身的复杂度上升了，所以总的改进不算大。那么，能不能将以上后缀比较的方法和某个高效的排序算法（如快排、基排）结合产生更快的算法呢?
为了回答这个问题，先来探究上面比较方法的适用条件。显然，P(i)的两维都必须是可比较的，即对于rank的比较必须有意义。这一条件自然是对快排的沉重打击，因为快排在进行到最后一个元素之前其顺序都是混乱的。那基排呢，看得出在排序完成之前也不满足有序这一条件。所以，上面的算法看来是不能再改进了。
但是，基排有一个绝佳的性质就是部分有序，即在进行第k+1位排序时，若是MSB基排已经按前k位有序了，若是LSB基排则已经按后n-k-1位有序了。假设适用MSB，当后缀按照前k位时，构造二元组P(i) = < rank(i), rank(i+k) >，显然此时P(i)的比较是有意义的（因为rank(i)和rank(i+k)有定义），所以可以直接得到在2k下的排序结果。
这样，每次排序的长度增大了一倍而不是1位，使得基排只需要O(lgn)趟就可以完成，总的复杂度为O(nlgn)。这就是发表在第一篇后缀数组论文中著名的Manber算法，04年国家集训队中讲述的倍增算法就是说它。
懒了，不想写了，还是把这篇文章分两部分吧，下次再说其它更好的算法。

我弱了，一直以为SAP算法是说的Edmonds-Karp算法（EK），现在才知道这个专指Ahuja和Orlin发明的一个非常高效的增广路系的最大流算法，根据实测，比之前常用的Dinic算法还要好一些，而且两者在代码上几乎无区别，所以非常容易实现。

网络流的基本知识我就不介绍了，直接对比Dinic说下SAP算法。Dinic利用了层次图的概念，使得每次找增广路不用从头开始，从而巧妙的使程序的趟数降到了n，复杂度下降为O(n^2m)。SAP则是利用了preflow的距离标号思想，进一步的免去了建层次图的过程，同时preflow的Gap启发也可以利用上，所以即使复杂度上限仍然是O(n^2m)，但实际跑起来要快许多。

这个算法怎么实现呢？回忆Dinic算法找阻塞流的过程，每次增广后，当前点退回到那条关键边的发出点w上，从w开始找新的可行弧。在这里，如果w没有可行弧，则退回到w的前一个点，而SAP在这一步上做了改进，先利用preflow中提高标号的方法（即将w的标号变为与w在residual network上相连的所有点中最小的一个加1），再继续寻找可行弧增广。

这样，有了提升标号过程，就不用重建层次图了。当源点的标号>=n时，程序就结束，或者当有Gap启发时，一旦发现某个高度出现断层，也可立即退出。

就这样，终于把POJ那道用我“优化”过无数次的Dinic却仍然过不到的Core Dual Cpu给过了。给我的教训：以后不能太自以为是，没了解过的算法就别太想当然。

最后，推荐大家尝试此算法。

本打算介绍Hopcroft-Karp算法（简称HK）和Unit Capacity Network Flow两种算法的，但无赖后者牵涉到的东西太多，就像和网络流有关的问题一样，主要讨论不同实现和启发对性能的影响。鉴于尔后有集中复习网络流的打算，这里便只说说第一种算法。

关于二分图匹配的匈牙利算法，现在一般作为图算法的入门学习内容，其原因就在于其编程简单，适合不求甚解式的强迫记忆，致使大多数人搞不懂如何从匈牙利树和增广轨入手设计算法（ms扯远了….）。匈牙利算法的理论复杂度是O(VE)，但实践中我发现一般是跑得很快的，特别是再贪心一个初始匹配后。不过，我们似乎更关心如何降低这个算法的理论上界，于是有了HK算法。

回忆匈牙利算法，算法一共进行V轮，每一轮都从左部的一个未检查点开始查找增广轨。改进的思想是一次性找到多条增广轨，这样就可以减少算法进行的轮数。那么，这样是不是真的能改进效率，先给出算法描述如下：

设初始匹配M为空；

每一轮用当前图中的所有未匹配点求出相对于M的不相交最短增广轨集合P，如果P为空，则M此时已经是最大匹配，否则用P种的每条增广轨扩展M。

那么，这个算法为什么高效？我的分析如下：

设S1^S2表示求S1和S2两个集合的对称差。

1. 设M是一个匹配，p是关于M的一条最短增广轨，p’是关于M^p的一条最短增广轨，那么 p’ >= p + 2；

   证明：几乎是废话，一条增广轨的长度是奇数，两个奇数当然至少相差2。

2. 一共只用进行O( sqrt(V) )轮。

   证明：设当前的最短增广轨长度L=sqrt(n)，显然到达这个长度只需要进行L/2个phase。设此时的匹配为M，而最大匹配为M*，则|M*| – |M| <= n / (sqrt(n) + 2) <= sqrt(n)。如果后续每阶段只共现一条增广轨，那么也只要sqrt(n)次就好，所以一共需要3/2sqrt(n)轮增广。

上面的证明用到了一个基本定理，即对于任意两个匹配M1和M2，如果|M1| > |M2|，则M1^M2有至少|M1| – |M2|条不相交的增广轨。

现在，HK算法的复杂度也很显然了，就是O( sqrt(V) * E )。从实践来看，的确效率也较匈牙利有所提高，虽然又引入了更大的常数。另外，这个算法的实现也相当简单，每次只要把所有未配点一起bfs一次就好，大概只用写50行左右。不过相对dfs实现的匈牙利，显然编程复杂度提高了很多，而且从目前经验来看比赛中一般没有什么题能卡贪心初始的匈牙利的，所以推荐想开阔视野的朋友看看这个算法。

最后我还想请教下懂这个算法的人。很多书上书该算法和Dinic算法的思想一致，我一直没想明白，有谁能提示一下么？

ps：字符串算法的东西过段时间会继续，这两天被无聊的事情折腾得要死，太久没更新blog也很不爽，所以先放一篇上来。

构思这篇文章是一段非常痛苦的过程，主要在于苦于无法用简洁的语言来表述我想说的话。KMP算法的介绍网络上可谓一搜一大把，目前引用最广泛的还谓Matrix67的那篇，不过就我当年初学时的记忆而言，还是读得晕乎乎的，但似乎还算有较强的可读性。

KMP本质上是构造了DFA并进行了模拟（DFA和NFA是什么我就不多说了，参考自动机理论的书籍），因此很显然一旦从模版T构造了自动机D，用D去匹配主串S的过程就是线性的。KMP最引人入胜的地方就在于构造D的自匹配过程，它充分利用了D是一个DAG的性质，使得构造过程也是线性的。

我在第一次学习KMP时就被告知，KMP的速度在平均情况下还不如brute force来得快，原因就在于随机字符串的自重复性并不高。尔后，又听说一些算法实际速度是KMP的好几倍，在平均情况下它们的复杂度是亚线性的，比如Horspool和Sunday算法，于是就倍感神奇，原来终日崇拜的KMP并非最强大的。是的，线性并非是串匹配算法的复杂度下界，我们完全可以做得更好，而后面我将会说到很多的最优算法。

那么，能不能将KMP变得更快呢？回忆下，KMP是DFA，如果能找到其对应的NFA，然后用惯用的位并行的手段来模拟NFA的执行，那么复杂度大约就是O(n/w)，其中w位机器字长。这一思路是可行的，形式化如下：

设D为NFA，则D(i)为1表示T(1),T(2)……T(i)是S(1),S(2)……S(j)的后缀。状态转移表的构造也很简单，设该表为P，则如下伪码：

For i = 1 to m

P[ T(i) ] |= 1 << i

End

假设从j进行状态转移到j+1只要如下一个位运算：

D = ( ( D << 1 ) & P[ S(j+1) ] ) | 1

于是，这样就将KMP的大大提高了速度，不仅是理论上的，也是实践中的。

KMP的另外一个研究方向是Extend KMP（以下简称EK），它是说求得T与所有的S(i)的最长公共前缀（LCP），当然，要控制复杂度在线性以内。

EK我第一次听说是07年baidu校园招聘的笔试题中，它当时的题目是求最长回文子串，当然这是一个耳熟能详，路人皆知可以用Suffix Array很好解决的问题。事后听一个同学说他写了三个算法：Suffix Array，Suffix Tree和EK，当时就不明白EK是什么东西，但又没当面问他，于是这个东西就搁置了很久。知道后来北大的月赛一道题说可以用EK来做，我才终于从03年林希德的文章中开始认识到它，就像KMP一样，这个算法也一下就吸引了我。

设Q(i)表示T和S(i )的后缀的LCP，P(i)表示T和T(i)的后缀的LCP，那么和KMP一样，我们试图用P来求得Q，而P可以用自匹配求得，并且和求Q的过程相似。

我以求P为例简要说明一下。P(2)就直接匹配即可，从i = 3开始，如下：

设k < i，E(k) = k + P(k) – 1，且对所有j < i，有E(k) >= E(j)。

那么，当E(k) >= i时，便可以推知T(i) = T( i – k + 1 )，于是如果P( i – k + 1 ) < E(k) – i + 1，那么P(i) = P( i – k + 1)，否则P(i) >= P( i – k + 1 )，从E(k)开始向后匹配到E(i)，有P(i) = E(i) – i + 1，并且更新 k = i；

还有就是E(k) < i，肯定有E(k) = i – 1，不过这个不重要，重要的是直接从i开始做暴力匹配即可得到E(i)，则P(i) = E(i) – i + 1，更新k = i。

希望我把EK说清楚了，不过这种东西还是自己推导一下有意思，而且记忆周期更长。

最后来罗列下题目。KMP的经典题目是POJ 2185，是要找最小覆盖矩形，如果你认为懂KMP了就去尝试它。EK的经典题是POJ 3376，有一定挑战；当然还有就是上面说的最长回文子串，提醒下用分治＋EK来做是其中一种方法。

嗯，下次打算说下Suffix Array，主要是它那个传说中的线性DC3算法，不过我现在还没把握能不能简单的把它说清楚，姑且认为可以吧。