Archive for 程序园

字符串匹配自动机

首先(不严谨地)定义自动机:自动机就是一堆状态的组合,给某个状态的自动机一个输入,它会按照自己的转移规则变到一个新的状态。事实上,可以把状态理解成顶点,把状态的转移理解成上面标着符号的边,在某个状态得到了某个符号的输入时,就转移到那条边指向的状态去。

先介绍单串匹配的自动机。对于M个字符的串S[1..M],一个用这个串进行模式匹配的自动机需要M+1个状态。设状态i输入ch时的转移为auto[i][ch]。状态转移时要满足:若当前在状态i,说明已输入的最后i个字符是S[1..i],若有多个状态满足这个定义,则应该在编号最大的状态。如何保证上述条件始终被满足呢?首先,auto[i][S[i+1]]=i+1,这是明显的。如果ch!=S[i+1],则auto[i][ch]肯定应该转移到某个小于等于i的状态,这个状态j应该满足:它是串S[1..j]是S[1..i]+ch的后缀,且j是最大的。找出j的方法是:计算每个i的“后缀状态”,也就是给自动机输入S[2..i]后自动机应该在的状态suffix[i],则当输入ch时应该转移的状态就是auto[suffix[i]][ch]的转移。计算suffix的方法是:suffix[i+1]=auto[suffix[i]][S[i+1]],suffix[1]=0。单字符串的自动机就是这样构造的。每次转移到状态M时,就找到了一个匹配。

值得注意的是,上述suffix数组所含的信息和KMP算法求出的所谓next数组其实是完全一致的。其实CLRS上也说,KMP算法就是通过少量的预处理实时地均摊O(1)地求出自动机的状态转移(大意)。也就是说,串匹配自动机也可以这样实现:首先通过预处理求出suffix数组,然后在状态i输入了字符ch时,需要转移到状态auto[i][ch],但是具体转移到哪个状态我们并没有预处理求出,我们只知道:若S[i+1]==ch则转移到i+1,否则转移到auto[suffix[i]][ch]。利用这两个条件,前者可以看做递归边界,就可以递归的求出每个需要的转移。可以证明的是这种求的复杂度是均摊O(1)的。

然后来看多字符串匹配的情况,其实和单字符串的时候差别不大,只不过状态之间不再是一个简单的线性结构了,而是一个Trie的结构,要根据刚才那个“后缀状态”的定义求出所有没有在Trie中直接出现的边应该转移到什么地方,同样也可以像KMP算法实时地计算那些转移。具体的方法都跟单串的没什么差别。

SPOJ上有两道题可以练习自动机:PSTRING是利用单串匹配自动机(或者KMP)进行状态转移的DP,WPUZZLES则是直接的多串自动机应用。

关于字符串匹配自动机,更详细的说明和更多应用可以找Maigo(王赟)关于Trie图的论文。

Comments (1)

构造Suffix Array的DC3算法

Suffix Array这个东西的确非常有用,很多关于字符串处理的题目构造一个SA就能迎刃而解。最近在SPOJ上看到了题目SUBST1,一眼看出来是用SA解的,兴冲冲的写了个nlogn的算法,连测试数据都没有出就交上去,结果没有悬念的TLE了(按说n=50000时nlogn算法应该可以不TLE的,比如说给50000个数排序肯定就不会TLE,可见SA的常数之大)……于是没办法,只好开始学SA的线性算法。

找到了这篇论文,似乎就是大家所说的三分法的“Skew Algorithm”的原始出处。作者在文章里说这种算法的基本原理是“Difference Cover”原则,并且在采样(Sample)的时候以3为模,所以这个算法叫DC3(Difference Cover modulo 3)。在本文里这个算法的名字就以原始论文为准了。

DC3算法的基本过程是这样的:首先构造出来一个superstring,这个串的每个字母都是原串的字母拼接上后面的两个字母,但这个串仅包含原串中下标模3不为0的位置。(递归地)求出这个串的SA,设为SA12。显然SA12中串的排列顺序与最终结果SA的排列顺序是相同的,只是SA12包含的串只有原串的2/3。

然后利用这个SA12和原串在线性时间内求出没有包含在SA12中的1/3的串的SA(也就是下标模3为0的那些串),设为SA0。求SA0的方法是:串的初始顺序是它们后面那个串在SA12中的顺序,然后再依据第一个字母进行一次RadixSort就可以了。也就是说SA0中的顺序可以仅用两个关键字来确定:第一关键字是第一个字母,第二关键字是后面的那个串在SA12中的位置。

现在得到了SA0和SA12,只需要把它们合并就可以得到最终的SA了,要想在线性时间内完成合并,需要在常数时间内完成比较。其实也不太难:当模为0的一个串与模为1的一个串比较时,需要比较它们的第一个字母与它们后面的那个串(都在SA12中),当模为0的一个串与模为2的一个串比较时,需要比较它们的前两个字母与它们去掉前两个字母以后的那个串(显然也都在SA12中)。

具体实现有点麻烦,我目前的代码基本上是完全抄的论文后面给的C++程序。不过运行速度真的很理想。因为程序里有一个很强的优化:当superstring中所有的字符都不相同时,可以直接用一遍RadixSort得到SA12,不用递归。经过这样的优化,递归地层数非常少,对于完全随机生成的长度为50000的英文字符串,求一次SA平均只需进行4次调用。所以实际时间效率非常高,常数不算大。

目前我对这个算法的理解还远不够深刻,所以讲的可能不易理解。有能力的话就看上面给出的原始论文吧,只看前三节就可以,后面讲了如何以牺牲一定时间为代价优化空间使用,我们是完全用不着的。

实例代码见上面的论文的附录A吧,我写得跟那个一样。

Comments (6)

最小树形图

最小树形图就是给定一个有向图,求以某个给定顶点为根的有向生成树(也就是说沿着这N-1条有向边可以从根走到任意点),使权和最小。

解决这个问题有一个O(VE)的算法,是这样的:对于除根外的每个顶点,选择一条权最小的入边。如果选出来的V-1条边不构成环,则可以证明这些边就是满足要求的答案。如果存在环,则将环中的边收缩成一个顶点。设x是环中的点,y不是环中的点,v为新的顶点,w0为上一步中x选择的入边的权值,则原来的权值为w的边(x,y)变成权值为w的边(v,y),原来的权值为w的边(y,x)变成权值为w-w0的边(y,v)。这时最小树形图的答案就是环的权值和加上收缩以后的图的答案。

为什么收缩了以后要减小权值呢?很简单,最小树形图需要每个顶点选一条入边,顶点x当前已经有入边了,如果以后再给收缩后的顶点x选择权值为w的(原来属于x的)入边,就意味着已经选择的那条权值为w0的入边可以去掉了,也就是说选择这条边的代价应该是w-w0,所以就做这样的修改。

不过如果不固定根,我还不太知道如何用同样的复杂度解决,谁写过请教教我。

(发现自己最近表达能力很差,这个文章说得清楚些。)

示例程序(TJU2248,用邻接矩阵写的,O(V^3),比较简洁):
tju2248.cpp

Comments (3)

欢迎到ioiforum讨论OI问题

最近大部分人可能上不去OIBH。所以ioiforum的猫欢迎OIers前去讨论OI问题。据我所知,ioiforum是有神牛出没的。所以说若你有疑难问题希望得到指点,到那里发帖应该是个不错的选择。

地址:http://ioiforum.bytecool.com/

(本文为友情广告)

Comments (2)

求最大权二分匹配的KM算法

最大权二分匹配问题就是给二分图的每条边一个权值,选择若干不相交的边,得到的总权值最大。解决这个问题可以用KM算法。理解KM算法需要首先理解“可行顶标”的概念。可行顶标是指关于二分图两边的每个点的一个值lx[i]或ly[j],保证对于每条边w[i][j]都有lx[i]+ly[j]-w[i][j]>=0。如果所有满足lx[i]+ly[j]==w[i][j]的边组成的导出子图中存在一个完美匹配,那么这个完美匹配肯定就是原图中的最大权匹配。理由很简单:这个匹配的权值之和恰等于所有顶标的和,由于上面的那个不等式,另外的任何匹配方案的权值和都不会大于所有顶标的和。

但问题是,对于当前的顶标的导出子图并不一定存在完美匹配。这时,可以用某种方法对顶标进行调整。调整的方法是:根据最后一次不成功的寻找交错路的DFS,取所有i被访问到而j没被访问到的边(i,j)的lx[i]+ly[j]-w[i][j]的最小值d。将交错树中的所有左端点的顶标减小d,右端点的顶标增加d。经过这样的调整以后:原本在导出子图里面的边,两边的顶标都变了,不等式的等号仍然成立,仍然在导出子图里面;原本不在导出子图里面的边,它的左端点的顶标减小了,右端点的顶标没有变,而且由于d的定义,不等式仍然成立,所以他就可能进入了导出子图里。

初始时随便指定一个可行顶标,比如说lx[i]=max{w[i][j]|j是右边的点},ly[i]=0。然后对每个顶点进行类似Hungary算法的find过程,如果某次find没有成功,则按照这次find访问到的点对可行顶标进行上述调整。这样就可以逐步找到完美匹配了。

值得注意的一点是,按照上述d的定义去求d的话需要O(N^2)的时间,因为d需要被求O(N^2)次,这就成了算法的瓶颈。可以这样优化:设slack[j]表示右边的点j的所有不在导出子图的边对应的lx[i]+ly[j]-w[i][j]的最小值,在find过程中,若某条边不在导出子图中就用它对相应的slack值进行更新。然后求d只要用O(N)的时间找到slack中的最小值就可以了。

如果是求最小权匹配,只需要把那个不等式反一下就行了。算法需要作出的改变是:lx的初值为所有临界边中的最小值,find中t反号。

示例程序(Ural 1076):
ural1076.cpp

Comments (4)

如此优雅的C++

当在程序中使用滚动数组之类的优化手段时,一般我们会先写出不用滚动数组的程序,用小数据调试、验证后再将它改成用滚动数组实现。这时,你会如何实现这样的改动呢?好吧,请看看,使用C++,可以将滚动数组实现得多么优雅,如何将一个只改声明部分、不改实现部分就引入滚动数组的优化。

程序请见:http://ddsgu.yo2.cn/183.html

Comments (2)

实时离散化:NOI2003 editor另解

NOI2003的editor(文本编辑器)一题的标准解法是“块状链表”。这个东西我以前没有写过,怎么想都觉得无法将这个东西妥帖地实现。找到的一个这道题的C++程序竟然有三百多行,让我看得着实很恶心。这道题能不能不用这种很不优美的数据结构来解决呢?答案是肯定的!

我把我解决这个问题时使用的技术称为“实时离散化”。首先解释离散化,一种合适的理解就是将需要处理的序列中的某些连续的区间当作一个元素来操作,以减少总的元素数量。对于这道题而言,如果我们事先能够知道“在哪些字母的前后会插入一个区间”“在哪些字母的前后会有区间被删除”,把这些字母称为“关键字母”,用一个关键字母来代表它与它后面的所有非关键字母,就可以大大的降低时间复杂度了。因为INSERT和DELETE的操作总数不超过4000次,“关键字母”的总数也不会超过4000*2个,这样就相当于所有的插入删除操作都在最多8000个(而非2000000个)元素组成的串上进行,当然就会很快了。问题是,我们无法在一个字母被插入时判断出它是否会在将来成为“关键字母”,也就是说,无法像某些题目一样只做一次初始的离散化,我们的“离散化”必须是“实时”的,也就是,必须实时的维护可以被一同处理的区间。

听上去有点玄妙,其实很简单。用一个很长的字符数组保存曾经被插入的所有字符串,当前正在编辑的文本用一个区间的数组来表示,每个区间代表它在前述字符数组中的位置。当需要在某个区间的中间插入时,只需按照这个被插入的位置将这个区间分裂成两个,再在两个新形成的区间的中间插入就可以了。同样,删除时也只需要按照需要被删除的文本两端将它们所在的区间分裂成两个,然后删除连续的若干个区间就可以了。设修改操作一共做了M次,则任何时候的区间总数不会超过2M个,所以每次操作的时间复杂度都是O(M),所有修改操作的总时间复杂度是O(M^2),而且常数很小,由于M<=4000,是完全可以承受的。INSERT和DELETE都实现了,GET也很好实现,挨个输出就行了,PREV和NEXT更是常数时间而已,最多可达到50000次的MOVE呢?如果每个MOVE的复杂度也都是O(M)似乎会有超时的危险。注意到这样一个明显的事实:若有连续的MOVE操作,只需要做最后一次的操作即可。经过这样的处理,需要实际做的MOVE操作大约也只是O(M)的级别。

对于这类在序列上进行添加、删除、反序等修改的题目来说,设修改操作有M次,序列的最大长度为N,块状链表的复杂度是O(M*N^0.5),“实时离散化”的复杂度是O(M^2)或者O(M^2+N),可见它在很多情况下是很有代替块状链表的可能的。对于某些题目而言,我认为它是相对于块状链表的一个更高层的抽象,所以编程要简单优雅很多(我的堆砌着大堆短函数的拖沓程序用了156行,其中核心的部分仅有不到100行),也不失效率。对于本题来说,我的程序可以P4 2.oG的电脑上最大数据1.6s左右,是可以AC的。但是,“实时离散化”还是不能完全代替块状链表,比如说当M与N同阶的时候就会退化成O(N^2)(因为已经失去了离散化的意义了),而块状链表仍然保持着O(N^1.5)。

示例程序:
editor.cpp

Leave a Comment

求最大流的Relabel-to-Front算法

除了用各种方法在剩余网络中不断找增广路(augmenting)的Ford-Fulkerson系的算法外,还有一种求最大流的算法被称为压入与重标记(Push-Relabel)算法。它的基本操作有:压入,作用于一条边,将边的始点的预流尽可能多的压向终点;重标记,作用于一个点,将它的高度(也就是label)设为所有邻接点的高度的最小值加一。Push-Relabel系的算法普遍要比Ford-Fulkerson系的算法快,但是缺点是相对难以理解。

Relabel-to-Front使用一个链表保存溢出顶点,用Discharge操作不断使溢出顶点不再溢出。Discharge的操作过程是:若找不到可被压入的临边,则重标记,否则对临边压入,直至点不再溢出。算法的主过程是:首先将源点出发的所有边充满,然后将除源和汇外的所有顶点保存在一个链表里,从链表头开始进行Discharge,如果完成后顶点的高度有所增加,则将这个顶点置于链表的头部,对下一个顶点开始Discharge。

Relabel-to-Front算法的时间复杂度是O(V^3),还有一个叫Highest Label Preflow Push的算法复杂度据说是O(V^2*E^0.5)。我研究了一下HLPP,感觉它和Relabel-to-Front本质上没有区别,因为Relabel-to-Front每次前移的都是高度最高的顶点,所以也相当于每次选择最高的标号进行更新。还有一个感觉也会很好实现的算法是使用队列维护溢出顶点,每次对pop出来的顶点discharge,出现了新的溢出顶点时入队。

Push-Relabel类的算法有一个名为gap heuristic的优化,就是当存在一个整数0<k<V,没有任何顶点满足h[v]=k时,对所有h[v]>k的顶点v做更新,若它小于V+1就置为V+1。不过这个我还没实现,因为现在算法的效率我已经很满意了,以后有空再实现吧。

今天实现的Relabel-to-Front比昨天的距离标号最短增广路又要高效一倍左右。最后感叹一下,网络流真是有趣啊……甚至都想买那本专门讲网络流的书来看看。

示例程序(还是USACO Training Ditch,不过都是自己出数据测的):
ditch(preflow).cpp

Comments (5)

求最大流的使用距离标号的最短增广路算法

求最大流有一种经典的算法,就是每次找增广路时用BFS找,保证找到的增广路是弧数最少的,也就是所谓的Edmonds-Karp算法。可以证明的是在使用最短路增广时增广过程不超过V*E次,每次BFS的时间都是O(E),所以Edmonds-Karp的时间复杂度就是O(V*E^2)。

如果能让每次寻找增广路时的时间复杂度降下来,那么就能提高算法效率了,使用距离标号的最短增广路算法就是这样的。所谓距离标号,就是某个点到汇点的最少的弧的数量(另外一种距离标号是从源点到该点的最少的弧的数量,本质上没什么区别)。设点i的标号为D[i],那么如果将满足D[i]=D[j]+1的弧(i,j)叫做允许弧,且增广时只走允许弧,那么就可以达到“怎么走都是最短路”的效果。每个点的初始标号可以在一开始用一次从汇点沿所有反向边的BFS求出,问题就是如何在增广过程中维护这个距离标号。

维护距离标号的方法是这样的:当找增广路过程中发现某点出发没有允许弧时,将这个点的距离标号设为由它出发的所有弧的终点的距离标号的最小值加一。这种维护距离标号的方法的正确性我就不证了。

由于距离标号的存在,由于“怎么走都是最短路”,所以就可以采用DFS找增广路,用一个栈保存当前路径的弧即可。当某个点的距离标号被改变时,栈中指向它的那条弧肯定已经不是允许弧了,所以就让它出栈,并继续用栈顶的弧的端点增广。为了使每次找增广路的时间变成均摊O(V),还有一个重要的优化是对于每个点保存“当前弧”:初始时当前弧是邻接表的第一条弧;在邻接表中查找时从当前弧开始查找,找到了一条允许弧,就把这条弧设为当前弧;改变距离标号时,把当前弧重新设为邻接表的第一条弧,还有一种在常数上有所优化的写法是改变距离标号时把当前弧设为那条提供了最小标号的弧。当前弧的写法之所以正确就在于任何时候我们都能保证在邻接表中当前弧的前面肯定不存在允许弧。

还有一个常数优化是在每次找到路径并增广完毕之后不要将路径中所有的顶点退栈,而是只将瓶颈边以及之后的边退栈,这是借鉴了Dinic算法的思想。注意任何时候待增广的“当前点”都应该是栈顶的点的终点。这的确只是一个常数优化,由于当前边结构的存在,我们肯定可以在O(n)的时间内复原路径中瓶颈边之前的所有边。

我程序做了很多优化(甚至包括动态内存分配的时候使用placement new进行优化),最终能做到比wxs写的dinic快了大约一倍,但程序只有一百二十多行,我对这个结果很满意。距离标号增广的确是一个非常优秀非常实用的算法。

示例程序(USACO Training Ditch):
ditch.cpp

Comments (3)

“彩球游戏”题解

郑重声明:本人并不自认为“会”此题,只是AC了,说说自己程序的思路而已。

原题:http://hi.baidu.com/astar/blog/item/fe22ab18c3c54b0635fa4192.htm

首先要对输入的数据进行预处理,按照“有连续的X个颜色为Y的小球”的格式,将连续的同色小球一起考虑。s[i]表示第i组小球的数量,type[i]表示第i组小球的类型。注意,所有连续的大于等于M个小球,都可以等价地当作M-1个小球来处理。

设计状态为v[i,j,k1,k2],表示第i组小球到第j组小球的左边多了k1个与s[i]同色的小球右边多了k2个与s[j]同色的小球时被完全消去的最小代价。平凡的状态包括:i>j时,状态的值为0;s[i]+k1>=M时,状态的值为v[i+1,j,0,k2];s[j]+k2>=M时,状态的值为v[i,j-1,k1,0];i==j时,状态的值为M-(s[i]+k1)。

对于非平凡的状态,只需要考虑最左边的一组小球或最右边的一组小球如何处理。以左边的这组小球为例,可取的策略有两种:使用M-(s[i]+k1)个小球将它们单独消去;设s[i]与s[x]同色,在[i+1..x-1]这个序列被先行消除后,将这组小球与s[x]合并后再消除。右边的小球也可以同样定义这两种策略。还有一种策略就是,当s[i]与s[j]同色时,将[i+1..j-1]这个序列先行消除,再将s[i]、s[j]以及k1+k2个小球同时消除。状态转移方程如下:

v[i,j,k1,k2]=min{
M-(s[i]+k2) + v[i+1,j,0,k2];
M-(s[j]+k2)+v[i,j-1,k1,0];
v[i+1,x,0,0]+v[x,j,s[i]+k1,k2] (if type[x]=type[i]);
v[i,x,k1,s[j]+k2]+v[x+1,j-1,0,0] (if type[x]=type[j]);
v[i+1,j-1,0,0] (if type[i]=type[j]&&s[i]+s[j]+k1+k2>=M);
v[i+1,j-1,0,0]+(M-(s[i]+s[j]+k1+k2)) (if type[i]=type[j]&&s[i]+s[j]+k1+k2<M);
}(i<x<j)

很麻烦吧……我也觉得很麻烦。 如果谁能简化一下那真是太好了。

这样的话,一共有O(N^2*M^2)个 状态,每组状态的转移需要O(N)的时间,总复杂度是O(N^3*M^2),不就严重TLE了么……话是这么说,也许这道题就能告诉我们Big O这种关于复杂度上界的分析有时是会不靠谱的。

具体实现时一定要采用自顶向下的DP实现,也就是所谓的记忆化搜索了。这时你会发现,实际需要计算出来的状态只占总状态数的一小部分。我对这题的4320个数据做了简单的统计,需要计算的非平凡状态的数量平均为2057.766个,其中此数量为0~1000的有2244个,1001~5000的有1541个,5001~10000的有393个, 10001~20000的有139个,20000以上的只有3个,最多的需要计算26621个状态。

我想这个题应该可以得到一个比O(N^2*M^2)更好的需要计算的非平凡状态数的上界,可本人算法分析的能力有限,希望有高人能够在这一点上指教一二。

最后说一下空间问题,如果开一个200×200x20×20的四维数组的话似乎太浪费空间了,另外,由于“内存分页”等我也不太理解的原因,程序的空间用量对实际运行时间(并非“用户时间”或CPU时间)是有影响的。我采用的方法是每次动态分配一个N*N*M*M的一维数组,用来模拟四维数组。例如状态v[i,j,k1,k2]在这个数组中的下标就应该是((i*N+j)*M+k1)*M+k2。当然,也可以将它理解成一种没有冲突的hash——将取值有限制的四元组映射为一个整数。对于本题的数据来说,我的程序最多时用了4128KB的空间。

于是,这道题的4320个数据就被比较完美的解决了。注意我并不是说这道题被完美的解决了,因为算法的理论复杂度仍然很高,不能保证无法设计出达到理论复杂度上界的数据。

欢迎批评指正。

zuma.rar 这是我做的cena评测包,内有我的程序,数据由陈启峰提供,共10组输入输出文件,每个输入文件中含432组数据。

Comments (5)