Try Again

当是日记好了

当是日记好了

趣的是，我还是没能遵守诺言，再忙也抽点空来写写blog。一晃四个月过去了，也学习了一些有趣的东西，又到了该总结一番的时候。
遗憾的是，本来想做一个纯粹的软工research，最终没能实现。这其中有值得总结的地方，首当其冲的问题便是我先入为主的认为一个问题它的复杂度高，所以我只要能把降下来那就可以发paper。其次，因为我的懒惰，所以惧怕去做一些不太兴趣的事情。所以，当我好不容易把一个很实际的问题提炼成图论模型时，我自然不想松手。
事实证明，其实naïve的方法，才是最好的方法：又好理解又好写。根据Occam剃刀原理，其它方法应该被排除了。可是，刚开始，我并没意识到这一点，因为我没有料到一个非常普通的数据结构： sparse bitmap，在实际中性能表现那么好。我承认我起初并非不知道这个结构，而是出于写起来比较麻烦，而没有去做。后来发现GCC里面实现了，我才赶忙拿出来尝试。记得第一次实验完的时候，我心是完全凉了，知道它会好没想到它会这么好。我心里一横，想反正没做多久，干脆抛弃了赶紧找新topic，离deadline还有一个多月，或许还来得及。
就这样，我悄无声息的开始尝试弄一些新的东西出来。后来我有了些小进展，跟老板说，他也觉得奇怪，怎么我突然就不做以前那个了？我没敢说我的实验，因为当初就是因为老板对我的信任，才让我尝试的，我不能告诉他失败了，我只能说contribution不够大，希望做一些results加强，这样更有希望。
对了，还没介绍我一直在做什么。pointer analysis，一个已经研究了四十年的东西，就是静态分析指针该指向何处，任何一个现代的编译器都会分析，只是算法不同分析结果的精度不同，因为广义上它是个undecidable问题。我选它并非是想challenge the whole world，我自知没这个能力，原因仅仅是这个问题很容易理解。参加过今年TIC决赛的同学估计还记得那道A题，便是exactly我所一直在研究的一个问题，然后出了题去恶心大牛的。虽然9月的ICSE我最终没能做出一些东西，但我肯定还会继续探索，因为我在PL差得还太多太多，做一个researcher不能操之过急，当然更不能临阵退缩。
不过，虽然高水平的研究是一个漫长的过程，但HKUST并不会因此而同情你，你得发paper，失败不得。ICSE过后，我又把我那个已经丢在垃圾桶中的算法捡起来，重新改了改并增加了一些证明，弄得还算fancy，投了出去。没报希望，但也必须尝试。
下一步，该往何处？自然，这有段空隙，是该补补基础的时候了。打算花一个月时间，把compiler中的一些基本问题再弄明白些，比如data flow analysis的lattice theory以及其应用。不能以后出去别人问起这些基本问题我心里还打鼓，讲得个是是而非的，那说我是做PL的不是个千古笑话么？
下一篇blog打算从compiler前端开始复习，第一个topic是DFA minimization和predict parser。看能不能把Hopcroft的O(nlgn)算法讲清楚。

算法学习

相隔近两月才来写这篇报告，估计大家早没有兴趣听这个题目了。最近一直在赶paper的档期，实在是静不下来写一些东西，恐怕这样的状态还要持续好一阵。

先回顾一下题目。一个图，每个点x有两个值f(x)和g(x)。两点x和y之间的距离定义为：d = ( f(y)-f(x) )/( g(y)-g(x) )。现在要往图里面添加边，规则如下：

1. 从所有异于x的点中选取一个顶点y，使得y到x的距离d是所有点中最小的，并且d ≥ 0。如果有多种选择，则选abs(f(y)-f(x))最小的一个(abs(t)表示取t的绝对值)，如果有多个abs(f(y)-f(x))相等，则选择f(y)-f(x)>0的那个点；

2. 添加无向边(x, y)（如果这条边已经存在，则不重复添加）。

现在要回答问题，每个问题询问x和y是否存在一条path。

大家都能很快看出，此题的后半部分就是一个并查集，难处主要在前半部分：这个图中哪两点间有边？因为题目的规模是有50,000个顶点，所以一个直接的O(n^2)算法是不行的。

其实，该题的障眼法我觉得设计得还是不够好：一看那个求d的公式，大家就基本上都猜到这个问题肯定和斜率有关。再加上题目的tips中说每个点的f和g值都不相同，于是基本上就可以确定这个题是一个计算几何了。

我坦白，这个问题其实就是要在平面n个点找形成直线后斜率最大的点对。首先把求d的那个公式倒过来一下，于是就变成了求1/d最大的问题。我们设每个点x的坐标为( f(x), g(x) )，很自然的把这些点表示出来。

然后就是把它们按照X坐标从小到大排序了。因为我们要求的是d ≥ 0（自然1/d ≥ 0）的直线，自然对每个点，只需要考虑如果以它为原点建坐标系，那么在1，3象限中的那些点。由于1和3象限的对称性，我们只说明1象限怎么做。

我们按照X坐标从大到小的顺序访问这些点。现在我们考察点x。我们用T(x)表示那些在X坐标比x大的那些点。如果要从T(x)找一个点y出来和x连线得到的斜率最大，那y肯定要是T(x)形成的凸包的上半部分。这个证明比较容易，因为假设y在凸包内部，那么我们知道x和y的连线一定和凸包上的某条边E相交了。设E的两个顶点是e1和e2，其中e1的Y坐标大于e2的Y坐标。那么我们可以知道，选择e1和x连线一定不会比选择y连线得到的直线斜率小。于是得到y应该从凸包上取。

有了这个分析，下一步就是如何从凸包上去找到y。很巧的是，这里面有个单调性，利用它可以设计巧妙的算法，见下图：

其中e1 – e5是凸包上的点，然后L1 – L5是x到这些点的连线。可能一个很直观的感受是，只要选凸包上Y值最大的点就行，于是上面的图就可以选出e3来。但其实不然，比如下图：

这就像一个人站在X处打探，e1和e2是两座山顶，那他肯定无法看到e2（即X – e1的斜率大于X – e2的斜率）。从数学的角度来描述这个事实非常容易：线段e2 – e1和e1 – X的叉积大于等于0。

于是，我们可以通过枚举凸包上连续的两点e1和e2，然后做一个叉积：如果大于0，我们知道e2和e2右边的所有点都被e1挡住了（因为上凸包的左转性质很容易证明）。于是这个单调性就很明显了：

如果e1挡住了其后面的点的“视线”，那么可知要找的最优点y一定是e1或者e1左边的点。依次，我们可以设计一个简单的二分算法来枚举凸包上的点，这样就可以使得求最优点的每步操作提高到lgn。

当然，计算几何题最烦人的地方就在于细节很多。以上给出了主要的思路，剩下的细节就大家自己想想。其实现在来看，C题明显比B的思路要简单，而且编程也不会更复杂。我本来也以为B应该是整套题最难的。不过从当时比赛的结果来看，大家似乎都只对B感兴趣，而忽略了C，我想主要还是因为B看起来要熟悉一些，而在比赛中一般都比较喜欢偷懒，有做过的题自然就要先照顾了。

夜深，剩下的一篇花絮希望尽快补上，不要再拖两个月-_-~~。

ACM解题报告 TIC, 二分, 凸包

非常荣幸今年能够有机会为这样一个大规模的程序设计竞赛出题，在全国为数众多的OIers和ACMers面前献丑，也实属鼓足了勇气，万幸的是整个比赛比较平静，看各个论坛和群也没有人指责我题目的问题，算是得到了些许安慰。

我计划本篇报告一共分成三次来写，第一和二部分分别是B和C的详细算法，第三部分是我关于这次题目在出题过程中的一个简单总结，或者叫花絮，希望愿意读的朋友来捧场。

好了，先来看看题目。B其实就是IOI2000经典问题Post Office的升级版。很多人可能都看过毛子青的那篇经典的将动态规划优化的论文，但惟独一个遗憾是他那篇文章在分析用四边形不等式优化该问题时最后得到的复杂度结论有错，那个算法应该是O(n^2)的复杂度而非O(np)，这也是为什么很多人吃了TLE的原因。

其实，更早有一篇IOI集训队论文就已经讲到了一个非常优美的算法，不过那篇文章非常简短，所以很多人可能没有在意。我标程的算法其实就是那个，在这篇报告中我稍微详细的讲讲。

首先对所有数据进行排序，这步证明很基本，反证得之，略过。然后列出基本的DP方程：f(k, j) = f(k-1, i-1) + W(i, j)。其中，f(k,j)表示从前j个城市中选择k个作为服务中心的最小开销，W(i,j)表示从数据排行第i的城市至第j的城市中选取一个作为服务中心的最小代价。假设排序后的数据存放在X中。

下面首先证明一个性质：使得W(i,j)最小的城市一定是第(i+j)/2个。假设P(k)表示选取第k个城市作为中心，如果k < (i+j)/2，那么有：

P(k+1) = P(k) + (k+1-i)*( X(k+1) – X(k) ) – (j-k)*( X(k+1) – X(k) )

= P(k) + (k*2+1-i-j)*( X(k+1) – X(k) )

因为 k < (i+j)/2，所以有k ≤ (i+j)/2 – 1，就有(k*2+1-i-j) < 0。因此得到P(k+1) < P(k)。同理，当k > (i+j)/2时也可以如法炮制的证明。于是原命题得证。

接着我们来回顾一下所谓的Convex Monge Condition，也就是常说的四边形不等式。任然用上面的W来描述，就是：

W(a,c) + W(b,d) ≤ W(b,c) + W(a,d)，对所有的a < b 且 c < d成立。

接着介绍的另外一个很类似的概念：Convex Totally Monotone，中文叫凸完全单调性。定义如下：

W(a,c) ≥ W(b,c) => W(a,d) ≥ W(b,d)，对所有的a < b且c < d都成立，其中=>是逻辑蕴含的意思。

这个东西有什么神奇的性质呢？来看两个：

1. 假设W是个n行m列的矩阵（以后都这么说），设Ri表示第i列的所有元素的最小值所在的行。那么有R1<=R2<=R3…<=Rm；

2. 四边形不等式能推导出凸完全单调性，但是反之则不然。

我现在要说是，这个题目中的W满足这个凸完全单调性吗？或者更甚，满足四边形不等式吗？这两个问题的答案都是肯定的。那么如何证明呢？因为后一个命题的更强，所以我们只需要证明满足四边形不等式就好，可以通过两步证明：

1. 对所有的i ≤ p ≤ j，有W(i, j+1)-W(i, j) ≥ W(p, j+1)-W(p, j)。

2. 将上式中的i和p替换为a和b，然后将j替换为所有区间[c,d]间的数，于是得到d-c+1个不等式，把它们加起来就可以得到答案。

下面就利用这个性质来设计高效算法。

回顾刚才的基本方程：f(k, j) = f(k-1, i-1) + W(i, j)。当k固定的时候，我们需要枚举i和j，在naïve算法中这是个O(n^2)的过程。假设矩阵B(i, j) = f(k-1, i-1) + W(i, j)，其中k ≤ i ≤ j，那f(k, j)其实就是B的第j列的最小元。因此，我们的目标就是加快求出B的所有列的最小元。

首先可以证明，当W满足四边形不等式时，B也满足。自然，B也同时满足凸完全单调性。我们用一个三元组（start, end, r）表示在[start, end]间的所有列的最小元在第r行上。算法开始时，我们只有一个元组(k, n, k)，即当我们只考虑了第k行的元素时，所有列的最小元都在第k行。接下来我们一次考虑第k+1, k+2……行。设当前考虑第p行，那么我们的工作就是在最快的时间内找到一个t，当所关心的列 ≥ t时，我们有他们第最小元在第p行。

首先，我们用一个队列维持当前的所有元组，并且满足下一个元组的start是上一个的end+1，即列的区间维持一个升序。然后初始化t = n+1，然后从最后一个元组开始考虑，那么只有三种情况：

1. 当B(end, p) ≥ B(end, r)时，我们知道对所有的c < end，有B(c, p) ≥ B(c, r)。这个就是凸单调性的逆否命题。这时我们不需要再向前扫描了；

2. 当B(start, r) ≥ B(start, p)时，根据凸单调性，对所有的c > start，有B(c, r) ≥ B(c, p)。这时我们知道该元组已经作废了，因为[star, end]中的所有列的最小元都不在r，至少有p就比它们小。于是弹出改元组，继续向前扫描；

3. 第三种情况就是，在[start, end]中存在一个t，使得对所有c ≥ t，有B(c, r) ≥ B(c, p)。如何来找这个t呢？通过二分查找就可以轻松的确定。这时也不需要再继续向前扫描。

然后我们就得到了这么一个t。如果t < n + 1的话，我们就加入一个元组( t, n, p)，否则就什么都不干。当最后一行n都计算完毕后，此时我们在队列中就维护者我们之前想要计算的答案：对于B的每一列，它的最小元在哪一行。

通过均摊分析，我们发现整个计算过程是O(nlgn)的复杂度。当然，B一共需要求P次，所以总的复杂度是O(Pnlgn)。

具体实现还是有好多需要考虑的细节，大家自由发挥。如果愿意的话，可以向我提出挑战我的程序（如果你有更好的算法，希望你能告诉我，因为我知道传说中好像真的有一个O(nP)的算法存在，只是我能力有限未能收获真知），然后咱们一边去晒晒，^_^。

这篇就到此为止，下一篇讲C 题。说实话，其实我更喜欢C一些，因为最终算法所需要的知识非常简单，只是中间的模型变化不容易想到，而B是如果没学过相关知识，从头到尾都不容易想。貌似比赛中C没有被人做掉，这点我还是有一些遗憾，希望大家能再想想。

ACM解题报告

很久没有来写blog了，眼看着有些荒废的迹象。自上篇文章以来，主要是经历了过年的颓废期，玩得有些忘乎所以，所以工作和研究上近乎停滞。

二月份没干几件正事，就想着到处去玩了。殊不知，玩这种运动是很容易上瘾的，于是荒废了大半个月。如果算上一月份的话，整两个月就仅仅把A. Aho的那本自动机和形式语言的前6章温习了，做了些习题。虽然至今又几乎全部还给了课本，但总算是圆了个心愿。

三月份重新拾起以前落下的研究计划。之前由于对Bug Pattern Matching没有什么好的思路，于是才放下去读书。我觉得我没有思路的根本原因不是缺乏研究的能力，而是从没有一头扎入代码中去看足够多的实例，而完成这项工作最好的便是能做一个静态代码的分析工作，于是我便能从代码中很从容的抽出所要的feature。而当我再次拿起这个课题时，我深感必须花时间去做这么个工具，或者是暂时放下，寻求另外一个问题。

那么，做什么好呢？ 我不知哪来的灵感，想到了parallel programming中的data race问题。说实话，我只是在之前的一个个人软件中遇到了这个问题，而且还是最容易解决的那种多线程竞争（我当时没用锁，呵呵，我不太懂啊）。其它时候我就根本没有碰过这个东西。现在要做，完全是来源于某个通宵的结果：偶然的在对着书柜发呆时，突然看到了那本已经落灰的Butenhof经典著作Programming with Posix Threads，然后故作久别重逢状的拾起翻阅了目录。其中有一章title，讲如何避免bug的那章吸引了我，见页数不多，便认真读毕。突然一种豁然开朗的态势，觉得data race是一种量化特征非常明显的bug，如果用control flow graph来表达程序，那其上的数据流方程可以很容易的detect常见的data race错误。

于是，我研究了一个晚上，总算有了一些结果。我不愿意知道是否人们已经做烂了这个题目，因为我只想训练自己的研究问题的能力，所以没有去读任何paper，所有的资料便是那本书。可以想象，大脑已经闲置了近两个月的我，此时看到这么一个可以做的题目，那是何等的激动啊。

第二天匆忙写好draft report，免得才想好的东西马上就消失了。接下来几天便是继续refine和elaborate的过程。这些天不好过，毕竟一晚上的东西有很多结论都是粗糙的。整个过程前后持续了13天，中间细节就不数了，最后完成了一篇长达19页的paper，看得自己发晕，看得老师也发晕。

当然，现在来看，这篇paper中真正有价值的东西不太多，因为是独立研究的结果，所以很多东西其实已经做得很成熟了。当然，还是庆幸自己有那么一次经历，至少开阔了自己的视野。

啰嗦了这么多， 那就列一下后续的一些写作计划吧。前一段时间由于完成一个小项目的原因，前后有进一个月时间认真分析了GCC的部分源代码，并做了部分改动。我打算写一些分析的总结，也能让自己日后有个参考。

另外就是一些读了些有趣的paper。现在开始走学术路线了，就得多做一些paper的总结。其实现在发现，原来读paper并不是一件很郁闷的事情，特别是那种研究内容很有新意的东西，往往还能给人醍醐灌顶的感觉，让你知道原来这个研究要这么开展。

当是日记好了

不知道大家是不是和我一样有这样的认识：一级指针和一维数组的语义行为是等价的，或者从编译器的角度看，两者是一样的。我想很少有人清楚编译器如何工作，又或许在实际工作中把两者当成同一物没有出现过问题，也就不会去细究两者的异同。

1. 什么是语义行为？

语义行为就是编译器用一系列操作对某个产生式的解释。对于一个指针int *p，它的语义行为有：

1. p所指向的元素类型是int;

2. p变量本身有一个地址，其值为另一个合法的地址，两者并无关系。

另一方面，int p[]有两层含义：

1. p 的类型是int[]，对p不做越界检查;

2. p的值就是p的地址（试试输出p和&p，看看两者是否一样）；

3. p[2]的计算方法是*(p的值+2*sizeof(int))，注意这里说的是p的值。

如果表达式为int p[10]，它相对int p[]要做越界检查。对于高维数组，其语义和一维的完全一样，只是地址计算上稍微繁琐。

2. 不区分一维数组和一维指针会导致问题吗？

一种常见的混用方式是：

void f( char *p ) { printf( “%d\n”, p[2] ); } int main() { int v[] = {1, 2, 3}; f( v ); return 0; }

这段程序不会有问题，因为参数传递是按值复制的，就是说p的值为v的值（就是v数组的地址）。但下面一段程序就要出错了：

A.c: int v[] = {1, 2, 3}; B.c: int main() { extern int* v; printf( “%d\n”, v[2] ); return 0; }

这段代码错误的原因在于：v作为数组时，其v的值语义为v的地址。现在v被告知是一个指针，于是v的值就是v的地址中所包含的值，上例中就是1。v[2]相当于1[2]，即取地址为3中的值，自然会segmentation fault。

3. 二维数组和二级指针的若干问题

看下面的示例代码：

void f( int **v ) { printf( “%d\n”, v[0][1] ); } void g( int v[2][3] ) { printf( "%d\n", v[0][2] ); } int main() { int v[][3] = { {1,2,3}, {3,4,5} }; int **t = v; f( v ); g( t ); return 0; }

很多人可能认为（包括我），二维数组的变量被编译器当做二级指针看待。其实这是错的，上例中f(v)的调用将会导致一个segmentation fault。根据数组的语义，int p[2][3]实际上是p指向一个int[2][3]类型的变量，据此来看看在f中发生的事情：v[0][1]实际上在尝试解析内存地址2，因为参数v的值就是传入数组的起始地址，所以v[0]就是1，然后解析地址1[1]，固然导致错误。

那么，如何改造f，可以得到正确的输出呢？f可以如下实现：

void f( int **v ) { int (*t)[3] = v; printf( "%d\n", *((t+1)[0]) ); }

但如果f在一个已经发布的函数库中（如getopt的第二个参数），那么我们便不能改动f。怎么办呢？此时就需要在调用f前做一件事情：

int **p; p = (int**)malloc( sizeof(int) * 2 );for ( i = 0; i < 2; ++i ) p[i] = (int*)malloc( sizeof(int) * 3 ); for ( i = 0; i < 2; ++i ) for ( j = 0; j < 3; ++j ) p[i][j] = v[i][j]; f( p ); // Release p

为什么这样？p与v不同就在于数据没有连续存放，中间引入了一个转换层。因此，p的地址解析需要通过两步完成，而v只用一步。除此，p还比v多浪费了2个void*的空间用于存放中间转换层数据。

上文中g(t)的调用是对的这个很容易看出，因为t只是起了个保存v的数据的作用，并没有做解析，然后将其传递给了类型匹配的参数。那用到t转换一次的意义在哪里呢？如果将g改造如下：

void g( int v[3][2] ) { printf( "%d\n", v[0][2] ); }

此时v的类型为int[3][2]而并非起初定义的int[2][3]，程序仍然能编译通过。原因就在于用t做转储时，丢失了范围语义信息，所以再转换回来时编译器已经糊涂了，也只好默认它们能正常工作，只是释放一个warning而已。

以上可见，从语义的角度来学习一个语言，会在理解上有本质的提高。对于遇到的“奇怪”问题，多做些这样的分析才可以得到正确的解释。

其它开发

离题，我觉得有必要先定义难和易的概念。

有些题目，它可能要经过好几步的模型转换，而且每一步的处理都需要一些比较高深的算法，比如网络流，自动机等。但是，其中的每一步的转化非常的straight forward，或者说用一些很普通的处理手段就可以看出解决方法。这类的题目非常多，特别是POJ月赛中。但我认为它们并不能算真正意义上的难题，能够做出不值得骄傲。

还有些问题，可能解决方法非常简单，代码量也非常小，但是却需要打破常规思路才能得到答案。我觉得就是要做这类题目，才可能在思维上有质的飞跃。

最近的两场SRM（423和424）就出现了两道我觉得值得一提的题目，下面介绍一下。

第一道题目是300pts。一个无限大的棋盘上放了n (n<=50) 个棋子，初始是任意的放置，所以也可能有位置相同的。问题是，最少移动多少步可以使得至少K个棋子在同一个位置上？其中从A移动到B的步数是两者的曼哈顿距离。

有一个很简单的问题我们或许都知道其答案：将n个棋子放在一起最少需要多少步移动？由于曼哈顿距离的特点，我们可以分开算X坐标和Y坐标。最后所有棋子移动到的点的X坐标一定是它们初始坐标排序后的中位数，Y坐标同理。因此，我们可以在O(n)的时间内解决将所有棋子移动到一起的问题。

现在考虑K<n。显然，普通的做法就是在n个中枚举K个棋子，然后用上面的方法计算移动步数。这样做的复杂度是O(C(n,k)*n)，显然不能接受。那么，如何优化？其实很简单，就是将这个不可接受的算法的计算顺序颠倒一下。仔细想想，由于X和Y坐标都是取中位数，也就是说最终要移动到的位置的X和Y的坐标是可以从已知的这n个点中组合出来的！因此，n个点最多有n^2个候选的位置，而对于选出的目标位置，很显然就是是n个点中离它最近的K个移动到一起。因此，通过简单的调换计算顺序避免了大量重复计算，将复杂度降为O(n^2KlgK)。

其实，这个算法也体现了动态规划思想的精髓：避免重复计算。然而，上面看似简单的算法，却非常不容易想到，原因就是我们先入为主的思想在作祟。

另外一个是600pts的题目。你在玩一个冒险游戏，现在有n个关卡。要通过一个关卡，玩家的当前智商和力量的值至少有一个要大于等于该关卡的限制。每通过一关，可以得到一些分数，你可以将这些分数任意的分配到你的智商和力量的属性上。初始你的智商和力量都为1，请问，如果通关的顺序没有要求，请问你最多能通过多少关？其中，关卡的属性值限制<=1000，n<=50。

初看真没什么想法，那就枚举吧。先枚举出通关的顺序，然后再求解这个顺序下最多能通多少关，复杂度是O(n!*n)。想基于这个思路动规，好啊，那总得状态压缩一下吧，答案是O(1000*1000*n*2^n)，没觉得这个方法有多好，怎么办？

其实，第一题的优化已经提示我们，盲目的枚举必然隐藏了大量的重复计算。回顾一下，第一题的重复计算在哪里？其实就在于，虽然一共有C(n,K)种点的组合，但是它们却一共只有O(n^2)中不同的最终集合位置。形式化一点，就是这指数阶的组合方案它们的某一个性质一共只有多项式阶种。回到这道题，其实就是所有通关顺序得到的最后的智商和力量属性的组合最多只有1000^2种（去掉>1000的情况，因为1000已经是所有关卡要求的上限）！

那么依葫芦画瓢，先枚举一种属性组合，然后再求出能得到这个属性组合的最优通关序列。很快发现，这样的话，不直接枚举(1000,1000)就可以通过所有关卡了？所以，再这一思路上，我们还应该多做些限制才好。

那么，为什么不能枚举(1000,1000)呢？原因就是这个组合或许根本就不存在。而题目一中不做这个处理的原因是，虽然有些位置不可能作为最终位置，但是计算它们不会影响最后答案。那如何快速判断哪些状态是可能出现的呢？

设状态(x,y)是可达的，则(x,y)必然是由某一个可达状态(x',y')通关后得到P分生成的（注意，假设P可以不用全部分配。因为这一假设不影响最后的答案，所以它是可行假设）。因此，很容易想到可以通过广度优先的过程从(1,1)开始，每次检测该状态能通过的所有关卡能得到的总分数，然后把剩余分数，即left=y+x-2分配后求得其它的可行状态。

答案很明显了，不再多说。

再总结一下这两道题目的特点。它们将一个具有指数阶以上定义域的问题转化到只有多项式阶的状态空间中，再在此状态空间中设计算法以得到原问题的最优解。我们可以看到，上面的算法的都不难证明，而且代码也容易写就，难就难在如何转化问题的状态空间。我觉得，只有多做这样的题目，才能迅速成长。

最后，即使你实现上面的算法很容易，也请抽空看看楼教主的代码，非常漂亮。

ACM解题报告

目的链接放在文末，大家可以试试再来看本文。

简述一下题意，就是有n个人，站成一个圈。每个人开枪击中别人的概率是固定的（设数组P[i]表示第i个人的命中率），从第一个人开始轮流开枪，求每个人生还的概率是多少？（假设每人都以最有利自己的方式选择射击对象，并且当有多个目标都一样好时，随机选择其中一个。并且，每轮都必须射击，虽然故意放弃射击机会可能增大存活的概率）。

一个引导我们想到方法的条件是，n <= 13且时限是4s，这提示我们这道题不像是个推数学公式的题目，并且很有可能是状态空间递推相关的方法。说到状态空间，必然就要设计状态的表示。直觉告诉我，状态表示是带集合的（因为n太小了），因此设计状态F[a][b][S]，表示当S集合的人一起玩这个游戏，并且下一个开枪的人是b时，a最终胜出的概率。

再来研究状态转移。如果a == b，那么显然b不会射击a（不允许自杀=.=），于是选择一个射杀后最大概率使a存活的人，设为c，于是有T = S ^ ( 1 << c )。再计算到c死后下一个射击的player为d，于是转移到F[a][d][T]，由此可得hit(a,b,S) = P[b] * F[a][d][T]；

另外，如果a != b，情况就不容乐观了。因为b为了使自己生存的几率更大，是可能选择射杀a的。假设现在有m个b中意的（何为中意？即射杀后b生还的概率最大）射杀对象（a是否是一员不用特设处理），对于对象k，设射杀后下一个射击的人是kk，集合变为T(k) = S ^ ( 1 << k )，于是转移到F[a][kk][T(k)]。最终有公式hit(a,b,S) = sum( F[a][kk][T[k]], for all k ) / m。

更有趣的是，如果b没有击中他想要击毙的人，怎么办呢？此时，状态集合S没有变，改轮到下一个人c射击了，于是转移到F[a][c][S]，有miss(a,b,S) = ( 1 – P[b] ) * F[a][c][S]。再进一步，如果S中的人射击完这一轮后，都没有命中，是不是又该轮到b射击了，也就是说，这里状态图出现了循环依赖，因此不能直接DP。

要解决循环依赖，其实也不难，假设现在集合中S有三个人a,b,c，从a开始射击，把方程列出来看看：

F[a][a][S] = hit(a,a,S) + ( 1 - P[a] ) * F[a][b][S]

F[a][b][S] = hit(a,b,S) + (1 - P[b] ) * F[a][c][S]

F[a][c][S] = hit(a,c,S) + ( 1 – P[c] ) * F[a][a][S]

移项再整理一下，得：

hit(a,a,S) = F[a][a][S] - ( 1 - P[a] ) * F[a][b][S]

hit(a,b,S) = F[a][b][S] - (1 - P[b] ) * F[a][c][S]

hit(a,c,S) = F[a][c][S] - ( 1 – P[c] ) * F[a][a][S]

 

这便是一个标准的常系数非齐次线性方程组，用高斯消元解决之。可以证明这个系数矩阵是可逆的，因此它有唯一解。还有一个小优化（其实狠重要的优化）是这个系数矩阵是比较特殊的，据此可以把GE跑到O(n)。

对于每个人x，最后的答案就是F[x][0][(1 << n) – 1 ]。

其实，写此文，我还有一个目的，就是我自知这个方法不是很优秀，特别是上述对于循环依赖的处理，疑有数学方法，或者，还有其它更好的状态表示就能自身解决这个问题，所以希望大家回复赐教，不胜感激。

 

题目链接：http://acm.pku.edu.cn/JudgeOnline/problem?id=3028

ACM解题报告

离上一次的发布已是一年有余，本来今年年初就打算要写，但后来不知为何放弃了计划，所以拖到现在，这也说明做事不得拖拉，越往后人越懒散。

好了，废话少说。新版本是完全重写的，因为上次的代码实在是不堪入目，因此结构性上有了很大的提高，有利于开发后续版本。功能上也做了大幅的增加，主要有：

1. 支持三种测试模式（两种数据对拍和一种标准OJ测试方式）；

2. 自动探测输入数据和输出数据的后缀对照关系，不必再向POJ那样一定要命名为xx.in和xx.out；

3. 支持从源代码自动编译，不必再手工编译后传入；

4. 支持Special Judge功能；

5. 统计运行中程序的时间和内存占用。

 

其它的更多细节请参见readme.txt文件。

应该说，基本上功能上还是比较完善了，急于这个程序进一步开发OJ也不是难事。当然不足肯定是有的，比如内存占用的统计目前就不是非常准确，希望有谁能告诉我这个应该怎么做得更好。还有就是可移植性方面考虑很少，只针对Linux的32位平台进行了测试，MacOS和BSD等系统和64位平台没有做，而代码中也没有体现对这些平台特殊的支持，因此后续的改善工作还是有很多的。

上传的代码打的Tag是RC2版本，正式版的发布请继续关注这个帖子。

另外，还有什么功能上的需求，请大家在回复中不吝赐教，是否切合实际不用考虑，毕竟现在功能有限，我还打算继续开发。

 

附使用方法简要说明：

设当前目录下待测试的程序叫test.c，输入数据在input文件夹中，输出在output中：

输入：

bash > tester -m 3 -i test.c -I input -O output

输出：

Compile ./test.c ....... OK

Test 0 = Accepted, time = 40ms, mem = 1228KB

Test 1 = Accepted, time = 36ms, mem = 1212KB

Test 2 = Accepted, time = 32ms, mem = 1212KB

Test 3 = Accepted, time = 36ms, mem = 1216KB

Test 4 = Accepted, time = 36ms, mem = 1216KB

 

如果要对拍程序，设另外一个程序叫std.c，数据生成器为gen.c：

输入：

bash > tester -i quick_sort.c -s std.c -g gen.c

输出：

Compile ./quick_sort.c ....... OK

Compile ./std.c ....... OK

Compile ./gen.c ....... OK

Test 0 = Wrong Answer, time = 316ms, mem = 2280KB

 

注意，一旦某个数据有错误，就不再继续生成新数据，可以使用-D选项转储错误数据。

下载： tester_v02.rar

其它开发

来打算这篇文章就开始讲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 构造算法

算法学习

必每个linux发行版都应该带了这样的小工具，我写这样一个小脚本的原因有3：

1. 字符界面，不必再进入X登录网络；

2. 避免每次都输入冗长的命令核密码；

3. 学习bash编程。

使用前请先保证当前用户有sudo使用网络界面控制命令的权利（iwconfig等），如果不行请先配置/etc/sudoer文件。

欢迎大家使用该工具，也欢迎大家提供自己的代码互相交流（回复留言即可）。

下载请点这里wuxian_ver1.0代码

=================================================

时隔一月，使用当中也发现了很多的bugs。针对这些，现发布过渡版本V1.5，以期解决如下问题：

1. 增加一些错误处理procedure；

2.  支持更多的功能，比如支持Ascii格式的password（这使得现有的密码保存文件的格式和原来的不兼容，建议手工在原有密码的最后空格再加上1）；

3. 修复若干bugs。

当然，还有一些问题没有解决啦，比如DHCP超时不能捕获，导致报道正确链接，这些我在研究当中。

代码下载在wuxian_ver1.5代码。

其它开发

隔半个月之久，我又看了许多关于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的后缀作为子问题排序即可。当然，更详细的数论的结论建议大家看看论文，简单而有启发。
最后，程序的代码在原始论文的后面有，写得很巧妙，基本上没什么可以改的（我试图自己写一个新的实现，结果始终不如论文作者思路清晰），因此我也不贴这里了。

算法学习