其实这个方法我去年11月就发现了，但当时觉得比赛已经结束，发布也没有意义了，就权当私下研究。不过转眼间几个月就去了，想想在这个圈子呆了一些时候，也得到了他人的帮助，想做些有意思的事情回赠，觉得这个还不错，所以就翻翻账。

实际上，想想Sparse Table（简称ST）算法和线段树算法（Segment Tree），就会发现，问题的解决之道就是对区间采取长短不同的覆盖策略，然后统计、合并子区间的结果。那么，Binary Index Tree（简称BIT，点击下载参考文献1）策略其实也是一种不错的区间划分方法，为什么不利用呢？

先约定下后续的符号。T是记录数据的数组，C是BIT的索引数组，用于存储与计算的子区间结果。现在的问题是求Range Minimum Query。

C[k]表示BIT中的第k个节点，它的值是[k – low(k) + 1, k]这个区间中所有元素的最小值，T[k]表示数组的第k个元素，low(k) = k & -k。那么k为奇数时就有C[k] = T[k]；k为偶数时，C[k]可以再被细分为lg(low(k)) 个子区间和T[k]元素本身。下面的代码就可以求得C数组的值：

[code:cpp]

void preprocess()
{
  int i;
  int y, k;  

  for ( i = 1; i <= n; ++i ) {
    y = LOW(i);
    C[i] = T[i];
    for ( k = 1; k < y; k <<= 1 ) C[i] <?= C[i-k];
  }
}

[/code]

代码非常简洁，同时很容易证明该预处理的时间复杂度是O(n)，常数是2，并且没有过多的操作，所以速度非常快。

查询也很简单。对一个查询区间[A, B]，如果B–low(B) < A，我们就访问T[B]的值，然后转移到B所管辖的子区间中去，即令B = B–1；否则访问C[B]，令B = B-low(B)，如此反复，直到B < A为止，代码如下：

[code:cpp]

int query( int A, int B )
{
   int answer = 0x7fffffff;   

   A = A – 1;
   while ( 1 )  {
        for (; B – low(B) >= A; B -= low(B) ) answer <?= C[B];
        if ( B == A ) break;
        anwser <?= T[B];
        B = B - 1;
    }

    return answer;
}

[/code]

代码很简洁，但是这个操作的复杂度却比较高，最坏情况为O(lg^2(n))，系数也是2。

这个算法有个常数级优化，观察到一个区间k的后半段一般都是零散的小区间，所以如果A < B – low(B) / 2，那么，实际上我们花了很多时间去访问后面的那些小区间，所以预处理时额外用个数组作记录，这样可以加快操作。代价是多花费了空间，以及代码变得不是很清晰。

总结一下。这个算法的时间是O(n) – O(lg^2(n))，空间是O(n)，而且时空上的常数都是2，代码也一如BIT其它操作一样简练，所以在实际中使用比起ST（内存量太高）和Segment Tree（编程较复杂，空间也比较多）在某些情况下是有优势的（部分测试数据见下表），大家不妨尝试一下。

 

本文算法和ST算法性能对比

硬件环境：

CPU: Centrino Mobile 1.6G

Memory: 512MB

Hard Drive: 5400RPM

符号：

n = 数组中的数据量

m = 查询数量

BIT = 本文算法

ST = Sparse Table算法

Time = 耗时

Memory = 内存用量

结果：

n = 100000,  m = 10000

BIT:  Time =         77.6ms, Memory =       1890.4KB

ST:    Time =         96ms, Memory =       7812KB

 

n = 100000, m = 10000000

BIT:  Time =         20589ms, Memory =       1892.8KB

ST:    Time =         11852.8ms, Memory =       7812.8KB

 

n = 1000000, m = 100000

BIT:  Time =         1003.2ms, Memory =       8920KB

ST:    Time =         1148.8ms, Memory =       79308KB

 

n = 1000000, m = 1000000

BIT:  Time =         4915.2ms, Memory =       8932KB

ST:    Time =         2415.2ms, Memory =       79308KB

 

n = 3000000, m = 1000000

BIT:  Time =         6733.6ms, Memory =       24552KB

ST:    Time =         4376.8ms, Memory =       254296KB

 

无向图最小割的相关算法及其延伸

这个题目实在太大了，因为这是一个研究了许多年发表了无数篇文章的问题，我自然也无法在短期内很好的诠释它。但为了与网上的其它文章有所区别，我打算把几个内容都揉到文章中去，但是这些都似乎在ICPC中尚无用武之地，就权当是爱好吧。
网络流算法自然是可以解决该问题，但是我们不作讨论。还有一类方法是90年代初期发展起来以缩图（Contraction algorithm）为基础的算法，目前最简单最高效的stoer-wagner算法是1997年发表的，它的主要思想如下：

定理1 对于图中任意两点s和t，它们要么属于最小割的两个不同集中，要么属于同一个集。
如果是后者，那么合并s和t后并不影响最小割。基于这么个思想，如果每次能求出图中某两点之间的最小割，然后更新答案后合并它们再继续求最小割，就得到最终答案了。
怎么完成求图中任意两点之间的最小割呢？不能借助网络流算法，否则我们说那么多就完全没有意义了。索性我们来个革命性的举动，对G做Maximal Adjacent Search，具体过程如下：
1.任选一个点s到集合A中，定义W(A, p）为A中的所有点到A外一点p的权总和；
2.对刚才选定的s，更新W(A,p)（该值递增）；
3.选出A外一点p，且W(A,p)最大的作为新的s，弱A != G(V)，则继续2。

我们把最后进入A的两点记为s和t，那么W(A,t)就是s和t的最小割的值了（后面为方便，把W(A,t)记为C）。做完这一步，就把s 、t合并了，于是图中又少了一点，然后重复，直到图中只有一个点位置，所以一共做了n – 1个phase，每做一次花费O(n^2)，所以一共是O(n^3)（当然，这个过程很象是个prim，所以对于稀疏图，可以用fibnacci堆把搞到O( m + nlgn )，而我是假设它是一个密图）。
还是继承我一贯的风格，对于这种论文上有的东西我不会给予证明，大家可以去搜stoer-wagner算法的原文来参考。
接下来，既然已经找到s-t最小割了，那能否快速构造出以它们为源汇的最大流呢？观察到实际上上面求得的最小割就是把所有于t相连的边都割掉，也就是说这些边最后都会被充满；s-t之间的边和为c，于是先将最大流变为C = C-c，去掉这些边，再按下面的步骤来做：
1.将图中每个点v的编号映射到它进入A的名次，如t -> n，s -> n – 1；
2.把图中每条边有向化：对于边e的两端点u、v，如果u > v，那么e化为 u -> v；
3.对图中每个点v，将v的所有出边按另一个端点的编号从大到小排序；
4.把t的所有边预先充满；对s边，按照刚才排序的顺序，尽量把C分配上去；
5.算法从v = n – 2开始执行下面的过程；
6.对于v的所有入边，记Fs为进入v的流量，Ft为流出v的流量，如果Ft >= Fs，那么到7继续，否则将Fs – Ft的流量像源一样把分配到出边中；
7.v = v – 1，如果 v = 0 那么到过程结束，否则到6继续。
上面算法的证明非常巧妙，几乎就和stoer-wanger算法一样有创造性，非常值得学习。而该算法的复杂度仅仅是O(m)，所以每次求了s-t割后再附带求个最大流并不会改变算法的总复杂度，非常优秀。

在最小割算法研究中，其实Monte-Carlo 算法占了相当的比例，Karger–Stein算法是其中的代表。它与stoer-wagner算法的区别就是不作那个最大领接搜索，仅仅依赖于定理1，每次从边集中随机选一条来合并，直到图中只有两个点为止。这样一个思路是很简单的，但它们对于该算法的正确性分析却是令人兴奋的，首先有如下定理：

定理2 任何一个图的最小割个数不超过C(n, 2)。

也就是说该算法能成功的概率是1/n^2，因此只要随机次数达到O(n^2)，成功率就相当可观了。
值得一提的是，这是随机算法中一个经典设计的案例，因为定理1最早就是由Karger提出并证明的，而后在95年提出了上面的算法，在97年它们又对该算法做了改进（思路不复杂，但是空间耗费比较大，所以略过），而此时stoer-wagner刚刚问世（据我所知该论文94年在欧洲某会议发表了，但是97年在ACM发表后才扩大了其影响力）。
另外，其实很多有效的方法来提高找min-cut的效率的，比如有名的PR Heuristic系列就是用来进行启发式判断，这些大家可以去网上找找。
就介绍这么多，相信本文除了提醒你有许多神奇的东西外，你是没有学到什么的，所以还是看论文吧。

<!– @page { size: 21cm 29.7cm; margin: 2cm } P { margin-bottom: 0.21cm本来打算把题目范围扩大到K Constrained Path的，但是最后还是放弃了，因为考虑到目前的能力有限，无法研究一个很深入的问题，所以选了一个比较常用且有代表性的议题来分析。

Loopless的意思就是要求无环路，一般来说，图不存在非正环则第一次求出来的路就是无环的。那么，一个图的Loopless S-T路有哪些比较显然的性质呢？

1. 不满足OP，即Kth Loopless Path 并非由 Jth Loopless Path 组成（J < K ），就是说这条路中的一段可能是这两点间的>Kth的路（例子很好找，也很容易想通，等看我推荐的论文也行） 。但是，如果对每个点维护一个已访问的节点表，那么这样一种压缩的状态就是满足OP的；

2. S –> T的Loopless Path只有有限条，其实这是无环的直接推论。

由性质1，我们已经发现传统的标号修订、设定以及上次说的Path Deletion算法都已经不中用了，那么，偏离路算法呢？看下面的定理：

定理 将yen算法中求从偏离点v出来到终点的与以前不重复的代价最小路改为求出从偏离点出来与之前不重复且不经过P(S -> v)所有点的代价最小，即可求出Loopless Path。

定理证明仍然是反证，略过，请参照我以后给出的论文。

既然是个现成的算法，那是不是直接套用以前的程序呢？就朴素的实现来说，因为要去掉已经历点求最短路，因此对每条选出的路径求其所有的偏离路的代价最少（按稀疏图计算）是O(V*E*lgV)，那么总就是O(K*V*E*lgV*lgK)，这显然是比较高的，因此，我们需要探究一些优化方案。

首先，原来的文章中提出的只在优先队列中保留K项的优化仍然适用，因此，优先队列还是选择平衡树；其次，这次必须要构造Pseudo Tree，因为我们无法预处理删除图中某些点后位于T点的SPT；再次，因为构造Pseudo Tree的缘故，所以就不应该将某点的所有偏离路都发展出来，而只能挑其中最小的一支构造，否则将会导致建树的复杂度太高。

最重要的是，有没有办法可以有效降低计算偏离路径的代价呢？下面我举例介绍一篇论文提出的方法：

设当前选出的路是1 -> 2 -> 3 -> 4 -> 5，偏离点是2，T = 5，S = 1。以前的方法是从2开始一直到4构造偏离路，现在我们反着来，从4 -> 2。计算4时，去掉图中1、2、3三点，计算一次位于5的SPT。该3了，这时不要重新计算一次SPT，而只要把4加入到刚才计算好的SPT中，松弛几次，就可以得到新的SPT了。就这样，一直到2，每次都把上次的点加入到图中，然后做一次SPFA。

这样做的确是有效果的，因为它基本上相当于在原图上做一次SPFA的代价，即O(EV)。但它的改进还不是本质的，不过我也暂时没发现其它优化，同时该方法的实现也很简单，只需要将Pseudo Tree存成Parental Link的形式（就是从叶子向根走的方式），而实际上这只是一个顺便的行为，因为Parental Link的树本来就更易于管理。

还有要注意的是偏离点可能已经向外延伸了好几条路线，因此再延伸时就要选与前面那些不同的，而其它点的处理相对容易些，因为它们都只延伸了一条边。不过，由于任何点都在某时可能成为偏离点，实际上也就无所谓一般与特殊。解决这个问题一般有两种，方案，一是采用位码表示哪些后继已经选过了，二是用链表把串起来。前者在图规模不大时有奇效，后者一般用于密图。

好了，算法就这些，这样做的时间复杂度是O(K*V*E*lgK)，空间O(K*N^2)。

最后，还是推荐做POJ 3137，很多东西还是要动手才知道有难度。需要提醒的是该题还要求比较路径的alphabetical序，因此在优先队列的处理上就不同于以前提出的合并同值项了，因为不再有两条路在该比较下相等。而关于alphabetical更多的内容，可以在上次那篇文章中提到的链接中找到，也欢迎回复讨论。

上一次讲述的算法非常容易实现以及扩展，但缺点也是显而易见的，就是时空效率很低，到最后感觉就是一个盲目的扩展搜索，毫无一点智能可言，倒是启发式搜索在此处有非常高实用价值，只是内存用量不容易控制。那么，是否还有更好的办法呢，或者说是求解该问题的思路？答案是肯定的，这就是我们要介绍的偏离路算法（该算法在LRJ黑书中提到过，号称主页上有讲，可惜我没发现）。

我写文章的习惯是绝对不会原封不动的把原算法呈现给读者，要么提出自己的理解，要么就是在原有基础上进行改进。偏离路算法系中有个很著名的就是Yen算法（不知道他是不是就是该派的创始人？），后面我将围绕它展开讨论，并提出自己的改进措施。当然你如果只想知道更多的关于Yen算法本身的内容，也可以参照ImLazy给出的http://imlazy.yculblog.com/post.1956603.html这篇文章及上一姐妹篇，可能会找到根多你想要的东西。

Deviation Foundations

还是照例介绍概念和基本理论，然后再引入算法。
偏离路的概念其实很形象，设想有两条S -> T路P1 和 P2，那么它们必然存在一点V，使得两者在S -> V这段完全重合，V -> T后开始分道扬镳，那么我们就说P2是相对于P1的偏离路，V就是偏离点，V -> T路段就是分离路。唯一需要注意的是我们定义P1的分离点就是S。

如果已经求出了S -> T之间的P1,P2,……,Pk，现在我们用它们来建立一棵树叫Pseudo Tree(PT)。这是个非常Amazing的想法，我感觉它的用途和性质都还没有发掘完，因此还应该有很大的应用空间（离题了？）。我想用下面边的图示来说明什么是Pseudo Tree，我语言匮乏，没办法用语言来描述它的伟大（用绿色标示的就是分离点）。

需要注意的一点就是，同一个点在PT中出现，你应该把理解成是多个不同的点，就是分层图那样。

这棵树究竟有什么用？可能目前还没什么发现，但是至少告诉了我们前K短路之间是有某种联系的，换句话说，就是似乎能由前K短路推导出前K+1短路。恩，我们的确找到了新思路，剩下的就是要证明一些有用的性质能辅助我们的推导了。

定理 我们用Tk表示由前K条最短路组成的Pseudo Tree，用P(k+1,V)表示第K+1条最短路相对于Tk的偏离路的从Vk下一个节点Vknext开始的路段，那么P(k+1,V) -> T 一定是从Vknext到T的最短路。

怎么证明？我提示下用反证就可以了，很简单，具体还是去看论文吧。
有了这个定理，Yen算法的框架也就出来了，用白话描述一下，再给伪码：
1.先求得S -> T最短路;
2.每次选取一条S -> T路径，这条路就是Pk;
3.对于Pk，把其从偏离点开始直到T之前的每一点想像成和某条新路径P’的分离点，因此分别把它们看作分离点，求出在该点作偏离的意义下不与已有路径重复的新S -> T路，那么P’因该作为Pk+1的一个候选.
4.如果k < K，那么回到 2。

恩，说起来就这么简单，不过可以看出来实现可能要用到优先队列，我用X表示，而对任意点U，用phi(U)表示U -> T的最短路，AT(x)表示x在ST上的那些边的集合：
[code:cpp]
k = 1;
put the S -> T shortest path into X

while ( k < K && X != empty ) {

P = gather the first path in X

for each node v From Vk to T

compute the arc (v,y) such that cost(v,y) + phi(y) is minimized over E(v) – AT(v)

construct this path and Insert into X as a candidate to (k+1)’th path

k = k + 1

}
[/code]

看明白了吗？或许没有，如果你属于此类读者，那么我强烈建议你去读读上面我提到的那个blog上的文章，然后再回过来继续阅读。

Implementation

实现部分应该就是本文的特色了，现在我姑且认为你已经明白了该算法，如果没有就请先不要继续 阅读下去了，回过去先仔细弄明白吧。

反正，如果想直接裸实现上面的算法，就我的实验看来，实际时空用量和编码难度都不十分理想。同时，还有重边的问题上面没有处理，看来都还是一件没有完成的艺术品。

定理 从任意一条路P改造而来的新的候选路径不可能比P的权值小。
推论 如果有不少于K条S -> T路径比某条路P的权值小，那么不需要再继续扩展P。

以上两条看来是很显然的，但同时也很好用：X中只需要保留K条路径，这样可以有效限制盲目扩展，保证时空效率。但这对X提出了新要求，要快速插入新元素，快速查找最大（判断是否X中最大的元素也比待插入的元素小）和最小（作为优先队列的基本用途）的元素，显然普通的二叉堆是不能胜任了，看来得请平衡树出山才能摆平。

还有个麻烦的问题就是偏离点老显得有点与众不同，因为从它出发在树上已扩展的边不止一条，如果为它专门编写代码似乎显得程序累赘。我采用了一个笨（一说迂回）办法来避免此问题，就是每次某点扩展新路径时把所有的边都扩展出来发到X中（不是每次都去选最小的那条边），这样每次我们取出新路径的时候都只需要从偏离点的下一点开始，于是程序自然就显得简洁多了。

而之所以敢如此改造，就是因为我们的X中只有K个项，自然可以抵挡很多无用的扩展，而没有这一性质的帮忙是绝不可行的（因为优先队列的空间开销将曾O(M)的速度增长）。

再来看内存。因为我们说了每条从偏离点下方开始的路都是一条最短路，所以我们也不需要记录这条路径，预处理一次就可以，空间复杂度是O(K)的，已经很优秀了，和Dijkstra算法同阶。

而时间呢？按照我改进的实现，应该是很比较精确的上限O( K*|E|*lgK)，而K是在不断减小，自然程序运行中这个上限是很难靠近的，故有很好的实际表现。

Further Improvement

在我阅读到Yen算法的文章中，作者还提到了一个他的改进，叫MPS算法，其主要还是利用重赋权技术。不过鉴于他的改造不能和我的一起使用（遗憾？不完全是，保留K项的方法还是可以用的，其实，这一思路非常重要，你还将再我以后的文章中读到），所以有兴趣的读者可以参考。

嗯，今天的内容就说完了。真正的偏离路思想的确是很微妙的，并非我几句就能道出，不过我说一些自己的理解总能起一些作用，也很欢迎在下面留言参与讨论。当然，POJ 2449 是个绝好的练习资源，不容错过。

开篇
以前写了两篇关于K短路问题的文章，其主要思想都是利用A*算法来解决问题。在第二篇文章中，我引入了分层图的思想，应该算一个非常容易理解和实现的模型。最近再次研究了这一个问题，但这次是读别人的论文而非自己研究，学到了很多新方法，现准备汇总于此。计划该文一共有五篇，其主题分别是：

1. Label Correcting algorithm
2. Deviation algorithm
3. Path deletion algorithm
4. Loopless path
5. Conclusion, preview and code

因为以前关于Label Setting算法的讨论还算比较彻底，因此这回我就不讨论它了（实际上，后面将讨论的新算法都是从这个思想衍生来的），直接从新内容开始。还有一个我曾经提到的很强的Eppstein算法，由于我还没有搞明白，所以不在此讨论。

有个前提是，前3篇文章都是说的Unconstrained的问题，并且是基于有向图的。
下面介绍后文要用到的一些符号：

S: 源点
T: 终点
P(u,v): 一条从u到v的路径（顶点可重复）
P_k(u,v): 一条从u到v的第k长的路径（顶点可重复）
(u,v): 一条u->v边
K: 欲求的第K长路径

Label Correcting Basics
有两个非常熟悉的算法是属于这一类型的，分别是Bellman-Ford和Floyd，其实它们都是可以推广来求K短路的。
回忆一下这两个算法的基本原理，就是所谓的最优化原理(Optimality Principle)，或叫最优子结构，复述如下：

定理1(OP)：如果P(S,T)是最短路，那么对任意t属于P(S,T)，有P(S,t)和P(t,T)都是最短路。

学DP和随机过程的时候都会介绍它，我就不多说了。

定理2(最短路存在条件)：有向图G中，S、T两点间的最短路若存在，当且仅当不存在一条S、T路径经过负环。

定理3(K短路的最优性原理)：令P_k(S,T)=P_x(S,u)+P_y(u,v)+P_z(v,T)，则x,y,z<=k。

这个定理证明很简单，我就不写了，可以参考我后面给出的论文链接，如果要简单理解，其实可以想想我曾经说过的层次图。不要看这个定理形式简单，它却明确的告诉我们所有的k短路都是由小于等于k的路段建立的。

Implementation
好了，还是说说Bellman-Ford是怎么推广到求Pk的吧。
对每个点x，定义F(x,i)表示到x的第i条路的长度。如果将F(x)排序，那F(x,i)就是i短路了。又由定理3，我们只要保留每个点的前K短路就行了。用P(x,i)表示(x,i)是否在队列中，这样避免多次松弛。
Q表示一个Queue，F(x).size表示取这个向量的长度，伪代码如下：

[code:cpp]

Q={(S,0)}

对所有(x,i),P(x,i)=false

while (Q!=empty) {

    (x,i) is extracted From Q

    for each edge goes From x,(x,y)

    L=F(x,i)+dist(x,y)

    if F(y).size<K  add a new label to F(y) and add this to Q

    else {

        (y,j) is largest label among all in F(y)

        if L<F(y,j) {

            F(y,j)=L

            if P(y,j)=false then Insert into Q

        }

    }

}

[/code]

答案就是F(T)中最大的一个，复杂度是O(Knm)。

是不是很自然？那么关于Floyd的推广自然也很快能写出: (F(x,y,i)表示x->y的第i条路的长度)

[code:cpp]

for x=1 to n

    for y=1 to n

        for z=1 to n

            对每一对F(y,x,i)和F(x,z,j)，取出来像上面一样更新F(y,z)

[/code]

读者很快发现，上面算法的朴素实现复杂度是O(K^3*n^3)，不过Clever一点的话是可以实现O(K^2*n^3)的。不过一般可以用扩展Floyd的场合K都不算大，而后一个实现其实引入了很大的常数，所以实际效果不一定优秀。

好了，两个算法就是如此简单，比起A*算法更容易实现。在ICPC比赛中，其实写这两个算法也是可以考虑的，毕竟编程复杂度是个不可忽视的问题。这一节的内容还比较简单，因为大家基本上对这两个算法还是很熟悉的。下一次我就重点说说偏离算法，这是个非常有趣的思想，时空性能好而且实现也简单，思路非常值得借鉴。

先说几个概念：Transient State (TS)与 Recurrent State (RS)
设从状态 i开始，经过若干次转移后还能回到i的概率为Pi，那么若Pi<1，则称i为Transient State，若Pi=1，则为Recurrent State。
因此，一个Ergodic Markov Chain（不熟悉此概念的请参照上次的文章）中的所有状态都是Recurrent。
Absorbing State (AS)是指进去后就不能出来的状态，即P(i,i) ＝ 1，可见AS为一个RS。
听上去有点像自动机？没错，其实Markov链本身就是一个Stochastical的自动机，因此不要忘了，在一个链中，AS和RS都可能有多个。
很自然的，我们所关心的问题形式化为：从任意一个给定的TS开始，最后进入某个AS的概率。
对任何一个带AS状态转移矩阵，首先通过行列交换把它变成标准式（分块矩阵）：

	  
          Q   R
P   =
	  O   I

即把SA全部放到下面。那Q就是TS -> TS的转移矩阵，R是TS -> AS的转移矩阵。
容易知道，恰好第n步从TS -> AS的概率为F(n) = (Q^(n-1)) * R，因此对所有的n从1到无穷求和，就可以知道转移到AS的概率为：
P=(I+Q^1+Q^2+Q^3….+Q^n)*R
因为Q中全为TS，所以Pi < 1，最后Lim( Q^n, n -> ∞ ) = 0，因此可以知道：
I+Q^1+Q^2+Q^3….+Q^n = W = Inverse(I-Q)
其中Inverse表示求矩阵的逆。假设所要求的目标矩阵为F，其中F(i, j)表示从第i个TS转移到第j个AS的概率，那么F = W*R。
OK，上面的推导省略了很多，读者可以自行演算。既然数学基础有了，那么就该考虑如何编程求解了，看下面一个例题。
从网上找到一个OI（据说是GHY出的）的题目，下面结合题目分析下。

猴子是非常聪明的孩子，更为可贵的是，她非常喜欢玩飞行棋。但是猴子的智商太高了，这种棋类游戏在她眼中“已经无法适应社会生产力的发展了”，所以，她决定加强飞行棋来锻炼一下思维。
她加强的飞行棋规则是这样的：
棋面有N个格子，一开始棋子在1号格子中。有一个色子，六个面分别涂上了0〜5这六个数字。每个格子上有一个数字，表示棋子到达当前位置之后向后移动的位数（当格子上的数字为0时表示棋子原地不动。当格子上的数字为负数时表示向前移动。保证1号格子中的数字为0。一轮必须移动到不能再移动为止）。每轮打一次色子，然后把棋子向后移动色子想上的面的数字。要使这个游戏胜利也很简单，只要跳出棋盘就可以了。（但是这个游戏有可能永远不会结束，此时就 Game Over）
由于猴子玩遍飞行棋未尝碰到敌手，阳阳常常非常惨烈的认输。现在阳阳想知道，对于一个棋盘，有多大的几率永远都不会结束。

【输入文件】
输入文件第一行包含两个整数N（1≤N≤100）。接下来给出N个数字，表示棋面上的格子数。输入文件末尾可能会有多余信息，请忽略。

【输出文件】
输出一个实数，表示最后失败的几率。精确到小数点后5位。

【样例输入】
4
0 1 -1 0

【样例输出】
0.40000

【样例说明】
当棋子落在2、3号格子中就会因为无法结束游戏而失败。
第一轮棋子落在0号格子中的几率为1/6，落在2、3号格子中的几率为1/3。
第二轮棋子再次落在0号格子中的几率为(1/6)2，落在2、3号格子中的几率为(1/6)*(1/3)
第三轮落在2、3号格子中的概率为(1/6)2*(1/3)……
总概率为1/3+(1/6)*(1/3)+(1/6)2*(1/3)……
因此根据无穷等比数列求和公式，失败的概率为(1/3)/(1-1/6)=0.4

此题是典型的Markov问题。
将棋盘按顺序编号，并在最后追加一个格子，代表棋盘外，那么可以写出如下的转移矩阵：

1/6 1/6 1/6 1/6 1/3
0 0 1 0 0
0 1 0 0 0
0 0 0 1/6 5/6
0 0 0 0 1

很幸运，该矩阵已经化成标准形式。我们尝试用高斯消元来做：

 F=Inverse(I-Q) * R  ⇔  (I-Q)*F = R

如果把R和F的每个对应列向量拿出来写上面的等式，那么每个等式都是标准的解一次线性方程组。而由于系数矩阵是固定的，所以这个过程可以只用一次消元来完成即可。
对数学敏感的读者已经注意到，就这题来说（其实这一现象很普遍），上面这个I-Q是不可逆的，但是可以不用管它（Why? 因为如果缩图以后不就行了），还是用全选主元来做，遇到零行跳过即可。那么，既然如此，其始就可以直接用I-P来消元编程不更简单（当然，如果AS很多，这样做的复杂度明显较高）？
因此，O(n^3)算法便计算出了，F是个列向量，F中所有元素的和便是答案。

好了，再来看一个稍微复杂一点但是和上面类似的题目。题目大意是两个人分别写一串数字，然后用一个随机发生器来产生数字序列，如果哪个序列被产生出来了，那么游戏结束，对应的人胜利。
此题是今年南京赛区的题目，比赛的时候只有少数人通过，据说用Markov来做的都因为精度过不了，而要用具体数学上面介绍的一个优美结论。但这不是本文重点，我还是来说说算法。
游戏取胜是指某人写的序列是随机序列的后缀。那么，自然想到把两人所写的序列分拆成若干个前缀的集合，那么每个前缀就是一个状态，而序列本身就是吸收态。设F(S)表示从S开始到达某个吸收态的概率，转移自然如下：

F(S)= 0.1*F(lcp(S+’0’)) + 0.1*F(lcp(S+’1’)) + 0.1*F(lcp(S+’2’)) + …..

其中lcp是找到状态集合中与S+’k’有最长公共前缀的状态。
再往后，就是求到达吸收态的问题了，完全可以用上面一样的方法解决，不再赘述。
至此，关于随机过程的学习先告一个段落，以后再深入研究，:P。

很早以前写过关于A*求k短路的文章，那时候还不明白为什么还可以把所有点重复的放入堆中，只知道那样求出来的就是对的。知其然不知其所以然是件容易引发伤痛的事啊，前天一次pku的比赛就得到了应验，因此下决心要把这一类的问题好好想一想。

搞了一天，总算有点成果，在此就总结一下这几类问题：

1. 两点间的最短路的条数；
2. 多关键字的极短路问题；
3. 给定长度的路的条数；
4. K短路及其条数；
5. 程序设计与实现。

在此，先说明几个符号。

label 是每次出堆的一个标号, 定义如下：

struct label {

    int vertex; // 该标号所代表的顶点

    int val;    // 标号的值
};

phi(x)：顶点x的标号，即源点到该点的最短路长度。
f(x)：源点到x的路径条数
S和T：分别是源和汇

—————————————–分割线———————————————

问题1，当用顶点s来作更新时，设对某条长为v的边s->t，我们做如下更新 ：
if phi(s)+ v == phi(t) then

    f(t) += f(s)
else if phi(t)>phi(s)+ v then

    f(t) = f(s)

不用再多说了吧，这不是今天的重点，略过。

 

问题2，多关键字的问题很具有实际意义的，就以双关键字为例子，其它的可以扩展。比如交通网络中的路一般都有收费和用时两个参量，在这种情况 下，就不存在一个最的定义，只有极值的问题。具体来说，若用(cost,time)来表示一个标号，那么显然标号不再是一个全序而是偏序。在堆中比较两个标号可以采用第一、 第二关键字分别比较的方法，因此每次出堆的一定是一个极小的标号。由于每个顶点可以存在多个标号，因此对每个点维护一个标号的集合。如果把这个集合排序的 话，那么一定是按cost严格递增，同时按time严格递减。

 

问题3，若要判断是否存在长度为L的path，设L’ = phi(T)，那么delta = L-L’,即对图中的每个点的标号来说，可以维护一个delta的冗余度（即某个标号比最小标号大delta也是可取的），在这样的长度下都是有希望最后到达T的长度是L。这只是一个必要条件，能否进一步限定？设T到 图中任意点x的最短距离为lambda(x)，那么对某个点x，只有标号val同时满足：
1. val – phi(x) <= delta;
2. val  + lambda(x) <= L 

时才有希望。进一步发现，式子中的lambda就是用作A*的启发函数，因此我们把问题3就归结到了问题4，下面一起处理。

 

问题4，这是本文的核心，我将详细说明。

我们知道A*对空间的要求是很高的，如果用扩展dijkstra的话则解收敛的更慢，因此更无法接受。如何优化扩展的顶点的数量就成了亟待解决的问题。

定义： 将上面的delta我们称为冗余度，把每个点x复制成独立的 x + delta 个顶点，表示为(x,0), (x,1), …….。这些点形成一个新图用G’。G’的形态是怎样的呢？假想把原图G复制delta份，第k份只由(x,k)的点组成，边连接的方式不变，即原来有x->y，在第k层就有(x,k)->(y,k)。相应 扩展其他概念，用phi(x,k)表示从(S,0)点到(x,k)点的最短距离，f(x,k)表示到达(s,k)的最短路数量。

现在我们来用(S,0)的最短路树在G’中是什么样的。假设我们直接对G’跑一个dijkstra算法得到phi和f。对某条长v的边(s,k)->(t,k)来说，若phi(s,k) + v == phi(t,k)，即该边连接的点在同层，保留该边。否则如果有phi(s,k) + v == phi(t,k’)，k’> k，即顺着这条边走到达了高层。那么我们删除 (s,k)->(t,k) ，建立(s,k)->(t,k’)。剩下的情况就是该边不属于最短路树，直接删除。

有了这个概念，不难知道求k短路条数其实就是求f(T,delta)，即从(S,0)->(T,delta)的最短路的数量，化为我们已解决的问题1。但这还没完，层次图毕竟是个概念，要真的建出来也不现实。在实现中要怎么解决呢？继续看下文。

 

问题5，最后就是编程的问题了，我用代码来表示比较清晰。

Run dijkstra in the level 0 graph   // 现在phi(x,0)就是已知的了。
while ( heap is not empty ) { 
 Get a minimal label (s, val) from min-heap
 For every out-going edge s->t of length v
 length = val + v; 
 k = val - phi(s,0); 
 theta=length - phi(t,0); 
 if ( theta <= delta ) {
 f(t,theta) += f(s,k);
 if (t,theta) is not in heap then Insert it into heap 
 }
}

从代码中我们也看出，其实编程中并没用到phi(x,k)，k>0，所以其实保留phi(x, 0)就可以了。

最后计算一下复杂度。|V’| <= delta *|V|，|E’| <= delta*(|E|+|V|^2)，所以总复杂度O(delta*V^2*(lgV+lg(delta))，这实际上也是在以前的文章中遗留的用A*求K短路的复杂度上限。

好了，问题都解决了。可是仍然不够圆满，因为新问题又来了，如果delta比较大该怎么做呢？还是有办法，可以利用专门的K短路算法解决（见我的K短路文章）。

以前学DP的时候，每本书都讲要能够DP的状态设计必须要满足无后效性，现在才知道这个后效性原来也有它的真名，就是Markov效应（或者叫Optimality Principle）。当然，我说这句话只是为了引出今天的主题，不过也可以窥探出他们之间存在某些联系。

简单来说，Markov决策就是带概率转移的递推，稍微形式化点可以写如下的方程：

F(S) = F(S’)*0.5 + F(S’’)*0.5

其中S’和S’’表示S可以通过两种不同的决策到达它们，选取各个策略的概率都是50%。当然，每次一说到递推，我们便会联想到用矩阵来描述它的状态转移，因此，再形式化一点，可以得到如下的表示方法：

对一个状态集合S，Fn为一行向量，Fn(Si)表示第n步转移后处在状态Si的概率，那么有：

Fn= Fn-1 * A， 其中A(i,j)表示 从状态i转移到 j 的概率。

对图熟悉的话，其实上面的状态转移可以看成是在一个有向图（这个图以后就叫Markov链）上进行的遍历，边上的权就是转移的概率，A不就是该图的邻接矩阵？

如果你理解不到我说的，那说明你可能需要补充点关于图的知识，我就不多说了。下面将介绍几个概念和定理，然后说下具体怎样应用到ACM/ICPC中。

定义一. 如果在Markov Chain中，每个状态都能在一步或多步走到其它任何一个状态，那么该Chain就叫可一遍遍历的(Ergodic)，或者叫不可约(irreducible)。

定义二. 如果存在m>=0，对所有n >= m，使A^n中的每个元素都为正，那么该Chain称为规则的(Regular)。

再简单点说，满足定义一说明该图的传递闭包的底图是一个完全图，定义二说不仅要是完全图，还表示从任何一个顶点出发，走m步以上一定能到达其它任何一个顶点。因此，每个规则的Chain都是可行遍的，但反过来不为真。

有了这些定义，马上就介绍两个非常神奇的定理，定理的证明我就不说了，比较简单（看来有时有想法才是关键）：

如果矩阵A是规则的Markov Chain，那么lim A^n（n->∞）= W。W中每一个行向量都相同，并且每个元素大于零，其元素和为1。

设W为上述的极限矩阵，w为W的一行向量，c为一全1列向量，那么有：

1. w*P=w，且所有满足x*P=x的向量x都有x=k*w；

2. P*c=c, 且所有满足 P*x=x的向量x都有x=k*c。

好了，现在我们把目光集中到w上。我们怎么能快速求出w呢？显然，模拟绝对是一个好方法，因为P^n收敛很快，因此取n=1000或者更大，那么就可以求出W，进而求出w。

这一点可以通过快速矩阵幂在n^3 * lgn以内完成。

还有一个方法。观察w*P=w，我们把两边同时取转置有 Trans(P)*Trans(w)=Trans(w)。将右边移到左边再合并，有(Trans(P)-I)*Trans(w)=0。这不是线性方程组么？

但是，先不要急，这还是个齐次方程，高斯消元的解显然是零解，但我们还有一个方程没用:

w1+w2+….+wn=1，如果把它带入并替换A中的某一个方程，这个时候再用高斯消元就可以得到解了，复杂度是O(n^3)，比上面的方法要速度快且精度高，但是编程要复杂一些。

补充一句，上面的w即是在矩阵A下的稳态（这个定义很重要，请区别我们以后要说到带吸收态的矩阵）。

最后来看个题把，POJ 1926，题目意思是图中每个节点都有一个权值，然后每一个时间段后，这个顶点都会把它的权植平均分配到其相邻的点上，经过若干分钟后，这个过程就稳定下来，即每个顶点的权值固定了，求这个稳定后的状态。

如果你是一路从上面看完本文的话，相信你肯定明白了要求的不就是这个图的稳态么？高斯消元吧，不过规模这么小，用法一也是能很快解决的。最后要注意的是该图不连通，那么你必须要把每个连通分量单独计算，然后再左乘一个初始的状态的行向量，就得到答案了。

今天说的只是Markov过程中很简单的一部分。今年南京的区域赛上有一个需要用到Markov来解决的难题，不过这篇文章的内容还无法解决那个题，下次的文章我会叙述怎样去解决那道题，有兴趣的读者可以先去找本讲Stochastic Process的书打打基础。