20世纪10最好的算法

一、算法一词的来源

　　Algos是希腊字，意思是“疼”，A1gor是拉丁字，意思是“冷却”。这两个字都不是Algorithm(算法)一词的词根，a1gorithm一词却与9世纪的阿拉伯学者al-Khwarizmi有关，他写的书《al-jabr w’al muqabalah》(代数学)演变成为现在中学的代数教科书。Ad-Khwarizmi强调求解问题的有条理的步骤。如果他能活到今天的话，他一定会被以他的名字而得名的方法的进展所感动。

二、20世纪10最好的算法

　　20世纪最好的算法，计算机时代的挑选标准是对科学和工程的研究和实践影响最大。下面就是按年代次序排列的20世纪最好的10个算法。

1. Monte Carlo方法

1946年，在洛斯阿拉莫斯科学实验室工作的John von Neumann，Stan Ulam和Nick Metropolis编制了Metropolis算法，也称为Monte Carlo方法。Metropolis算法旨在通过模仿随机过程，来得到具有难以控制的大量的自由度的数值问题和具有阶乘规模的组合问题的近似解法。数字计算机是确定性问题的计算的强有力工具，但是对于随机性（不确定性）问题如何当时并不知晓，Metropolis算法可以说是最早的用来生成随机数，解决不确定性问题的算法之一。

2. 线性规划的单纯形方法

1947年，兰德公司的Grorge Dantzig创造了线性规划的单纯形方法。就其广泛的应用而言，Dantzig算法一直是最成功的算法之一。线性规划对于那些要想在经济上站住脚，同时又有赖于是否具有在预算和其他约束条件下达到最优化的能力的工业界，有着决定性的影响(当然，工业中的“实际”问题往往是非线性的；使用线性规划有时候是由于估计的预算，从而简化了模型而促成的)。单纯形法是一种能达到最优解的精细的方法。尽管理论上讲其效果是指数衰减的，但在实践中该算法是高度有效的——它本身说明了有关计算的本质的一些有趣的事情。

3. Krylov子空间叠代法

1950年，来自美国国家标准局的数值分析研究所的Magnus Hestenes, Eduard Stiefel和Cornelius Lanczos开创了Krylov子空间叠代法的研制。这些算法处理看似简单的求解形为

Ax=b

的方程的问题。当然隐藏的困难在于A是一个巨型的n*n 矩阵，致使代数解

x=b/A

是不容易计算的(确实，矩阵的“相除”不是一个实际上有用的概念)。叠代法——诸如求解形为

Kx(k+1)=Kx(k)+b-Ax(k)

的方程，其中K 是一个理想地“接近”A 的较为简单的矩阵——导致了Krylov子空间的研究。以俄罗斯数学家Nikolai Krylov命名的Krylov子空间由作用在初始“余量”向量

r(0)=b-Ax(0)

上的矩阵幂张成的。当 A是对称矩阵时，Lanczos找到了一种生成这种子空间的正交基的极好的方法。对于对称正定的方程组，Hestenes 和Stiefel提出了称为共轭梯度法的甚至更妙的方法。过去的50年中，许多研究人员改进并扩展了这些算法。当前的一套方法包括非对称方程组的求解技巧，像字首缩拼词为GMRES和Bi-CGSTAB那样的算法。(GMRES和Bi-CGSTAB分别首次出现于1986和1992 SIAM journal on Scientific and Statistical computing(美国工业与应用数学学会的科学和统计计算杂志)。

4. 矩阵计算的分解方法

1951年，橡树岭国家实验室的A1ston Householder系统阐述了矩阵计算的分解方法。研究证明能把矩阵因子分解为三角、对角、正交和其他特殊形式的矩阵是极其有用的。这种分解方法使软件研究人员能生产出灵活有效的矩阵软件包。这也促进了数值线性代数中反复出现的大问题之一的舍入误差分析问题。 (1961年伦敦国家物理实验室的James Wilkinson基于把矩阵分解为下和上三角矩阵因子的积的LU分解，在美国计算机协会(ACM)的杂志上发表了一篇题为“矩阵逆的直接方法的误差分析”的重要文章。)

5. Fortran最优编译程序

1957年，John Backus在IBM领导一个小组研制Fortran最优编译程序。Fortran的创造可能是计算机编程历史上独一无二的最重要的事件：科学家(和其他人)终于可以无需依靠像地狱那样可怕的机器代码，就可告诉计算机他们想要做什么。虽然现代编译程序的标准并不过分――Fortran I只包含23，500条汇编语言指令――早期的编译程序仍然能完成令人吃惊的复杂计算。就像Backus本人在1998年在IEEE annals of the History of computing 发表的有关Fortran I，II, III的近代历史的文章中回忆道：编译程序“所产生的如此有效的代码，使得其输出令研究它的编程人员都感到吓了一跳。”

6. 矩阵本征值计算的QR算法

1959—61年，伦敦Ferranti Ltd.的J.G. F. Francis找到了一种称为QR算法的计算本征值的稳定的方法。本征值大概是和矩阵相连在—起的最重要的数了，而且计算它们可能是最需要技巧的。把—个方阵变换为一个“几乎是”上三角的矩阵――意即在紧挨着矩阵主对角线下面的一斜列上可能有非零元素――是相对容易的，但要想不产生大量的误差就把这些非零元素消去，就不是平凡的事了。QR 算法正好是能达到这一目的的方法，基于QR 分解， A可以写成正交矩阵Q 和一个三角矩阵R 的乘积，这种方法叠代地把 A=Q(k)R(k) 变成 A(k+1)==Q(k)R(k) 就加速收敛到上三角矩阵而言多少有点不能指望。20世纪60年代中期QR 算法把一度难以对付的本征值问题变成了例行程序的计算。

7. 快速分类法

1962：伦敦Elliott Brothers, Ltd.的Tony Hoare提出了快速(按大小)分类法.把n个事物按数或字母的次序排列起来，在心智上是不会有什么触动的单调平凡的事。智力的挑战在于发明一种快速完成排序的方法。Hoare的算法利用了古老的分割开和控制的递归策略来解决问题：挑一个元素作为“主元”、把其余的元素分成“大的”和“小的”两堆(当和主元比较时)、再在每一堆中重复这一过程。尽管可能要做受到严厉责备的做完全部N(N-1)/2 次的比较(特别是，如果你把主元作为早已按大小分类好的表列的第一个元素的话！)，快速分类法运行的平均次数具有O(Nlog(N)) 的有效性，其优美的简洁性使之成为计算复杂性的著名的例子。

8. 快速Fourier变换

1965年，IBM的T. J. Watson研究中心的James Cooley以及普林斯顿大学和AT＆T贝尔实验室的John Tukey向公众透露了快速Fourier变换(方法)(FFT)。应用数学中意义最深远的算法，无疑是使信号处理实现突破性进展的FFT。其基本思想要追溯到Gauss(他需要计算小行星的轨道)，但是Cooley—Tukey的论文弄清楚了Fourier变换计算起来有多容易。就像快速分类法一样，FFT有赖于用分割开和控制的策略，把表面上令人讨厌的O(N*N) 降到令人欢乐的O(Nlog(N)) 。但是不像快速分类法，其执行(初一看)是非直观的而且不那么直接。其本身就给计算机科学一种推动力去研究计算问题和算法的固有复杂性。

9. 整数关系侦查算法

1977年，BrighamYoung大学的Helaman Ferguson 和Rodney Forcade提出了整数关系侦查算法。这是一个古老的问题：给定—组实数，例如说x(1),x(2),...,x(n) ，是否存在整数a(1),a(2),..,a(n) （不全为零），使得

a(1)x(1)+a(2)x(2)+...+a(n)x(n)=0

对于n=2 ，历史悠久的欧几里得算法能做这项工作、计算x(1)/x(2) 的连分数展开中的各项。如果x(1)/x(2) 是有理数，展开会终止，在适当展开后就给出了“最小的”整数a(1)和a(2) 。欧几里得算法不终止——或者如果你只是简单地由于厌倦计算——那么展开的过程至少提供了最小整数关系的大小的下界。Ferguson和Forcade的推广更有威力，尽管这种推广更难于执行(和理解)。例如，他们的侦查算法被用来求得逻辑斯谛(logistic)映射的第三和第四个分歧点，b(3)=3.544090 和 b(4)=3.564407所满足的多项式的精确系数。(后者是120 阶的多项式；它的最大的系数是257^30 。)已证明该算法在简化量子场论中的Feynman图的计算中是有用的。

10. 快速多极算法

1987年，耶鲁大学的Leslie Greengard 和Vladimir Rokhlin发明了快速多极算法。该算法克服了N体模拟中最令人头疼的困难之一：经由引力或静电力相互作用的N个粒子运动的精确计算(想象一下银河系中的星体，或者蛋白质中的原于)看来需要O(N*N) 的计算量——比较每一对质点需要一次计算。该算法利用多极展开(净电荷或质量、偶极矩、四矩，等等)来近似遥远的一组质点对当地一组质点的影响。空间的层次分解用来确定当距离增大时，比以往任何时候都更大的质点组。快速多极算法的一个明显优点是具有严格的误差估计，这是许多算法所缺少的性质。

三、结束语

　　2l世纪将会带来什么样的新的洞察和算法？对于又一个一百年完整的回答显然是不知道的。然而，有一点似乎是肯定的。正如20世纪能够产生最好的l0个算法一样，新世纪对我们来说既不会是很宁静的，也不会是弱智的。

基本的搜索算法

响应李老师的号召，抽空写了一下基本的搜索算法。我的搜索技术很有限，写的东西很难足够充实与深刻，但结合了具体的POJ 的题目，相信对于初学的队员多少还是有一点用处的。我主要说了四个部分：
1． 二分搜索 （Binary Search）
2． 记忆化搜索 (Dynamic Programming)
3． 广度优先搜索(BFS)
4． 深度优先搜索（DFS）

另外，有一点我觉得初学的人要理解，搜索就是穷举，当在做题的时候你发现一个题目需要枚举所有的情况的时候，你在做的就是搜索，当然，搜索也有很多策略和技巧可以优化你的搜索，这就是经常听说的剪枝。惭愧，我的搜索基本没剪枝。根据不同的具体题目剪枝策略很不相同，该怎么剪枝我也就不再多说了，但是在搜索的时候要时刻考虑有哪些剪枝条件可以加上以优化搜索，否则就跟我一样只能平庸的搜索了。

1． 二分搜索 （Binary Search）

二分也是搜索，因为解空间具有单调性，所以可以利用二分的思想去掉一大部分解空间，也即做了一个强有力的剪枝。同时二分也是一种基本的算法思想，当在思考一个问题的时候如果不能灵光一闪可以尝试用一些算法思想来考虑，比如对于一个问题你可以这样思考，这题贪心是否可以呢如果贪心的话该怎样贪呢，找不出合理的贪心方法的时候再想，这题动态规划是否可以呢，有没有第推关系呢，还找不出来再想用图论模型是否可以呢，该用哪种图论模型呢最短路径，最小生成树还是网络流呢，如果还不行还可以想那用二分可不可行呢，该怎样二分，为什么是可行的呢。
什么时候二分是可行的呢？下面简单的探讨一下。
先来看这样的一个小问题，给一个函数 F(x) 在闭区间 [A，B]上单调递增（假设F(x)的定义域均为整数），要求找出在区间[A，B]上最小的一个满足 F(n)>=0 的整数n。对于这个小问题最直观的想法就是从A开始穷举，这样写出来的程序可能会是这个样子。

For (int n=A; n<=B; n++)
{
If (F(n) >= 0)
Break;
}
这就是传说中的纯暴力，但是这个问题是可以用二分优化的，它之所以可以二分是因为 函数 F(n) 满足单调性，这点十分重要，这是一个问题可否二分的先决条件。那为什么满足单调性就可以二分呢，这是因为如果n满足F(n)>=０的话，所有　Ｂ＞＝ｘ＞＝ｎ　的数肯定也满足。所以所要求的目标解肯定在区间　［Ａ，ｎ］中。二分程序我一般都这样写。

Int min = A; // 解空间的最小值
Int max = B;//解空间的最大值
Int ans; //存储目标解
While (min <= max) //当解空间非空
{
Int mid = (min +max) / 2; // 取当前解空间的中值

If（ok(mid) ）// 如果mid 满足要求
{ //继续在 [min, mid-1] 中寻找更优的解
Max = mid-1;
Ans = mid; //记录下当前找到的解
}
Else //否则在区间 [mid+1,max]中寻找最优解
{
Min = mid + 1;
}

}
Ans 即为想要的答案，实际应用时注意在区间里找不到最优值得情况，注意看题目说明这个时候应该输出什么值。

下面结合一道具体的POJ 的题目来说下二分的应用


POJ 3273 Monthly Expense

题意不再多说，说得抽象一点，这个题可以这样思考。
设函数 f(x) 表示 如果要使每个月的钱都没有超过 x ，则划分出来的月的个数，设题目的解空间为[min, max] （相信理解了题意很容易确定）则，题目的要求可以描述为在区间 [min, max]中寻找最小的值ans 使得 f(ans) == M; 当然如果找到最小的ans 使得 f(ans)=m < M 的话肯定也满足要求，很显然，因为用m月来划分尚且没有使得每个月的钱超过 ans 那用 M个月划分就更不会超过ans，于是题目要求更新为在区间[min,max]中寻找最小值 ans 使得 f(ans) <= M； 上面说过要是可以二分必然满足函数 f(x) 在区间上单调， 这个函数 f(x)显然单调，理由，对于解空间的两个值 x1 < x2 ，如果用x1作为每个月钱的上限划分出来的月数为 y1，则用更多的钱 x2作为每个月钱的上限划分出来的月数 y2没有理由比y1会小，所以原函数满足单调性，所以可以用二分来解决。现在只剩下一个问题了，对于一个任意值x判断它的可行性，也即函数f(x)的值该如何求。方法就是尽量往一个月中加入尽量多的连续的天数使得当前月的钱数总合没有超过 x，如果加入一天的钱超出了x则函数值++，置当前月的钱数为该天的钱数，继续进行下去。
具体实现不再多说了。

推荐的题目：

3122：分析方法与上面类似，推荐。
3179：二分+枚举区间，有点难度，有兴趣的队员可以做下。
2922：综合题有一定的难度，BFS+二分+枚举区间，有兴趣的队员可以做下。
2． 记忆化搜索 （Dynamic Programming）

动态规划其实也是搜索，只不过，在搜索过程中它把一些已经搜索到的结果记录下来，当再次需要时如果发现当前状态的结果已经记录下来的时候，没有必要再次搜索，直接取出值来就可以。
来看一道POJ 的例题
POJ 3186


注意到题目的描述，FJ是按照时间顺序，也即天来卖treat的，每一天FJ手中的treat都有一个特定的状态。设为d[i,j]，表示为当前FJ手中的treat 标号为从i到j ，即 i,i+1, i+2,…,j，初始时状态为d[1,n]，每一天FJ都要做一个决策，选择卖出一个treat，对于状态为d[i,j]时他只有两种选择，要么卖标号为i的treat要么卖标号为j的treat，所以如果让dp[i,j]表示当FJ面对的状态为d[i,j]时他所能得到的最大值，则他所能得到的最大值是他所做的两种决策的最大值，则
dp[i,j] = max(dp[i+1,j]+wt(i,j,i), dp[i,j-1]+wt(i,j,j))；由此得到了第推关系，其中wt(i,j,i)表示当状态是d[i,j]时，卖标号为i的treat，标号为i的treat能贡献的值，具体是什么样子的可以看题意。
有两种方式可以实现，一种是第规的方式，一种是非第规，如果用第规的方式的话，更能体会为什么会叫作记忆化搜索，实现的时候，用一个数组dp[i，j]记录下搜索得到的当前状态的最优值，可以赋初值为-1，设函数 f(i,j) 表示寻找dp[i,j]的最优值，可以在进入函数的时候检查dp[i,j]是否为-1，如果不为-1则表示它的值已经计算出来了，直接返回dp[i,j]的值即可，如果为-1则表示它还没有被计算，返回f(i+1,j)+wt(i,j,i)和
f(i,j-1) + wt(i,j,j)的最大值；

这种第规形式的动态规划，我印象中几乎没做过，有我也是改成了非第规的形式。如果谁知道的话可以共享一下。
3． 广度优先搜索 （BFS）

这种题太多了，从易到难，数不过来，所以如果学会了广搜的话你会感觉你可以AC掉很多很多的题目。如果没学会的话，赶紧着学吧，我就推荐几个题目吧。

如果你刚学会BFS，或者不太会的话，可以做一下这几个题目。
POJ 1562，2386，3194，1915，1979
如果你觉得基本的BFS题可以顺利的AC了可以做下这几个题目。
这些题就是考得耐心和细心，其实都不难。
POJ 1101 传说中的连连看
POJ 2935 简单的题目但要求输出路径
POJ 2046 记录状态判重比较麻烦
POJ 2049 不是很容易，数据比较阴险，我做了好长时间。


4． 深度优先搜索 （DFS）

深搜是在万不得已，不得不枚举所有状态的时候才会用到，而且往往要用到比较巧妙的剪枝。我就推荐几个题目吧 （不用剪枝）
BFS推荐题中的1562，2386，1979都可以用DFS做，可以体会一下。
POJ 2676 数独问题，硬搞就行不难，
POJ 1011 很难的题目，要求有精妙的剪枝，推荐早晚把它做了。
POJ 2362 类似POJ 1011 但比那个要简单一些，推荐先做2362再作1011。

基于多项式快速傅立叶变换的大整数乘法

【关键字】

多项式乘法 快速傅立叶变换(FFT) 离散傅立叶变换(DFT) 大整数乘法 分治策略 ACM 高精度算法库
【摘要】

大整数的高精度运算是一个很重要的算法课题。这篇文章讨论了如何实现大整数的快速四则运算问题。目的是为了建立一个大整数算法库。主要讨论了如何实现 乘法。
【正文】

一：引言

大整数的四则运算中，加法和减法没有太多的优化余地；对于除法，则应充分用到乘法的优势，这要通过巧妙的运用牛顿迭代算法。因而乘法运算是大整数四则运算中的核心内容。进行乘法运算时，如果我们模拟平时进行整数乘法的步骤复杂度将达到O(n^2)，应用快速傅立叶变换可以将复杂度降到O(n*lg(n)) ，从而大大提高程序的运行效率并扩大了使用范围。本文先从理论角度介绍了大整数乘法运算的FFT(Fast Fourier Transforms)算法。并从编程角度讨论了一些实现时应注意的边界和精度问题。（针对ACM/ICPC程序设计比赛设计，因而某些地方为了提高效率而牺牲了一部分其它性能）

二：大整数的表示方法

一个N进制的整数可以表示成SUM(ai*(N^i))。也就是说整数运算问题可以转换为多项式的运算问题。例如：对N进制的整数X=P(N)和Y=Q(N)有XY＝R(N)。（其中P(N)和Q(N)是以N为变量的n-1次多项式X和Y的各位是其系数，且P(N)*Q(N)=R(N)）

所以下面的讨论就针对多项式的乘法运算展开。

我们都比较熟悉多项式的这种以系数ai来表示的形式。但我们今天更敢兴趣的是它的另一种表示形式——点值表示法。由n个不同的变量x及其对应的多项式取值y所构成的n点可以唯一确定一个n－1次多项式。这也就形成了在这种联系基础上的另一种多项式表示方法：用这些点(x,y)表示多项式的形式。

三：在多项式系数表示法的基础上做乘法运算

如果我们要把两个多项式乘起来。用系数表示法如右图：我们需要进行O(n^2)次基本乘法运算。这也是我们所手工计算是经常采用的方法。这种方法看上去已经没有什么可以优化的了。因为它看似没有做任何多余的事情。其中的每一组乘法运算都会对最终的结果产生影响。而且我们不能通过调整计算顺序或使用什么技巧降低算法的时间复杂度。最多只能通过一些语法和其它编程技巧使得其常数因子稍微小一些。

四：在点值表示法的基础上做乘法运算

对于两个多项式Y1(x),Y2(x).他们的一组点值表示法形式分别为(x[0],y1[0]),(x[1],y1[1])...(x[n-1], y1[n-1])和(x[0],y2[0]),(x[1],y2[1])...(x[n-1], y2[n-1]) 在这种形式的基础上做乘法运算就可以用O(n*lg(n))的复杂度得到多项式的成绩。因为Y(x)=Y1(x)*Y2(x)在x[i]处的值是y1[i]*y2[i] 。同理可求出其它点的值。

要把n－1次多项式和m－1次多项式做成绩，我们只需要在n＋m－1个不同的的x点分别求处对应的y1，y2值。该点所对应的y值就为y1×y2。只需做 这样的运算就可以得到乘积多项式的点值表示形式。

现在，速度瓶颈已经转移到如何在点值表示和系数表示之间进行快速转换。如果任给n个x值算y值从而得到n个点的话时间复杂度为 。把这n个点值表示法的点对应成为n－1次多项式也得花费相同复杂度的时间。那么我们的努力也就徒劳。怎么找到更好的转换算法就成了解决大整数乘法的关键所在。

五：用快速傅立叶变换和傅氏逆变换对多项式的两种形式进行 的转换

我们利用单位虚根之间的内部联系通过在这些特殊点上进行抽样(取值)以降低采样的时间复杂度。为了后面使用分治策略解决这个问题我们不失一般性的假设多项式的次数为n－1且n为2的幂(低次多项式可以看作高次系数为0的高次多项式)

六：利用FFT的大整数乘法

只要将多项式中的变量看成一个固定的常数N。那么各项系数就对应着N进制各位的数字。所以我们进行大整数相乘的方法如下：

1：把一个乘数(整数)的个,十,百,千。。。位分别放入a数组中的第0,1,2,3…项

2：把另一个乘数按同样的方法对应到b数组中。

3：对于a，b数组分别求FFT得到虚数数组A，B。（按n位处理，n为2的幂，不够的话高位补0）

4：虚数数组C中C[i]=A[i]*B[i]。（该乘法为虚数相乘需要自己实现运算符重载）

5：对C做傅立叶逆变换得到c。

（本来是虚数组，但由于整数乘法这个实际问题的解一定是整数故最终c是整数数组，虚部舍掉，实部四舍五入）

6：这个c即对应着大整数成绩的结果

七：复杂度问题

每次运算时a被分为两个和原来相同大小的数组，经FFT后在以O(n)的速度进行重组得到A，整个傅立叶变换（FFT）和傅立叶逆变换的时间复杂度如为:O(n*lg(n)).

八：实现过程中应该考虑的问题

1：整数溢出问题

由于傅立叶变换法类似于先做不进位乘法，然后处理进位问题。故当两个数位数最大(出于现有计算机计算速度等方面的原因，这里考虑1,000,000以内的数 )且都为N－1时，中间结果最有可能溢出。在C＋＋语言中，N<=43时中间结果不超出int型，N<=61时中间结果不会超过unsigned型。

2：精度问题

由于只进行log(n)次迭代，而double型有18位有效数字，因而误差远远小于0.5。因而不会对运算结果产生影响。

3：空间问题

由于处理位数比较大，并且存储效率比较低（每一位占4个字节）占用内存间比较大。因而当位数接近处理上限（这里考虑1,000,000）时，局部变量已经不能使用了。只能用new去申请空间。否则会在运行时出错。

注：源码可在http://victorcheng.wy8.net下载

【参考文献】
[1]:Thomas H.Cormen, Charles E.Leiserson, Ronald L.Rivest and Clifford Stein. Introduction to Algorithm. The MIT Press, second edition.

[2]:李庆阳，王能超，易大义，《数值分析》