算法

推荐视频：另类素数筛选法

/ 7 条评论

  

    画一个三角形阵列，初始时只有每一行的两头有标记，然后从上到下依次把每一行复制两份，摆放成一个等边三角形。最后你会发现，第i行为空行（除两头外不再有其它标记），当且仅当i为素数。对于其它行，标记的位置也与该行行号的质因子有关。这是为什么呢？
    照惯例，给个YouTube链接：http://youtube.com/watch?v=sbjPwyPT1AI

    设f[i,j]表示第i行左起第j个位置是否有标记。j从0开始计数（即第i行最左边用f[i,0]表示）。对于每个f[i,j]，我们将它的值赋给了f[i+j,j]和f[2*i-j,i]。也就是说，对于每组i和j，我们都进行以下两个操作：
f[i+j,j] <- f[i,j]
f[2*i-j,i] <- f[i,j]
    而这实际上就是辗转相除的变形，不断递归下去后，最终f[i,j]表示的其实就是i和j是否互质。这样一来，上面那些东西就全解释清楚了。
    用这种方法描述数论问题是一件很有趣的事，给人感觉很神奇。如果你感兴趣的话，这里有一个类似的例子，你可以看到Sierpinski三角形、杨辉三角和组合数的奇偶性是如何联系在一起的。

计算阶乘的另一些有趣的算法

/ 15 条评论

    一个正整数n的阶乘就是前n个正整数的乘积，我们通常需要n-1次乘法操作来算出精确的值。不像等差数列求和、a的n次幂之类的东西，目前求阶乘还没有什么巨牛无比的高效算法，我们所能做的仅仅是做一些小的优化。

更少的乘法运算次数？
    在高精度运算中，乘法计算的速度远远慢于加减法，因此我们有必要减少乘法运算的次数。下面我将做一个非常简单的变换，使得计算阶乘只需要n/2次乘法。继续看下去之前，你能自己想到这个算法来吗？

    我们可以把一个数的阶乘转换为若干个平方差的积。例如，假如我想求9!，我可以把前9个正整数的乘积写成这个样子：
   1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9
= (5-4) * (5-3) * (5-2) * (5-1) * 5 * (5+1) * (5+2) * (5+3) * (5+4)
= (5-1) * (5+1) * (5-2) * (5+2) * (5-3) * (5+3) * (5-4) * (5+4) * 5
= (5^2 – 1^2) * (5^2 – 2^2) * (5^2 – 3^2) * (5^2 – 4^2) * 5
    注意到一个有趣的事实：上面的四个平方差算出来分别是24, 21, 16, 9，它们之间的差正好是连续的奇数（因为n^2等于前n个正奇数的和）。因此，我们可以用初始数(n/2)^2不断减去一个个的正奇数，求出所有n/2个平方差，再用n/2次乘法把它们乘起来。这种算法实现起来非常简单，并且（当n不大时）同样只需要单精度乘高精度，但需要的乘法次数大大减少了。假设我们已经有了一个高精度类，求n!只需要下面几句话：
long h=n/2, q=h*h;
long r = (n&1)==1 ? 2*q*n : 2*q;
f = LargeInteger.create(r);
for(int d=1; d<n-2; d+=2)
   f = f.multiply(q-=d);

更少的总运算次数？
    尽量提取阶乘中的因子2，我们可以得到另一种阶乘运算的优化方法。这很可能是不需要分解质因数的阶乘算法中最快的一种。
    假如我们需要计算20!，我们可以把20拆成若干组正奇数的乘积：

  1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 2 * 4 * 6 * 8 * 10 * 12 * 14 * 16 * 18 * 20
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 2^10
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 2 * 4 * 6 * 8 * 10 * 2^10
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 2 * 3 * 4 * 5 * 2^15
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 3 * 5 * 2 * 4 * 2^15
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 3 * 5 * 1 * 2 * 2^17
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 3 * 5 * 1 * 2^18

    只需要一次累乘就可以求到每一组奇数的乘积，最后再花费log(n)次乘法把它们全部乘起来。最后的那个2^18也可以二分计算出来。真正的代码还有很多细节上的优化，另外还借用了递归使得操作变得更加简便。你可以在本文最后附的那个链接里去找Split-Recursive算法。

还能再快一点么？
    继续扩展上面的算法，我们可以想到，如果把每个数的质因数都分解出来，并且统计每种质因子有多少个，我们就可以多次使用二分求幂，再把它们的结果乘起来。注意这里并不是真的要老老实实地去分解每个数的质因子。对于每个质数x，我们可以很快算出前n个正整数一共包含有多少个质因子x（记得如何求n!末尾有多少个0么）。这种算法的效率相当高，已经能够满足大多数人的需要了。

另一种诡异的阶乘算法：
    这个算法可能是所有有名字的阶乘算法中最慢的一个了（Additive Moessner算法），它对一个数列进行重复的累加操作，一次次地计算前缀和，总共将花费O(n^3)次加法操作。但是，令人费解的是，这个简单的程序为什么可以输出前n个正整数的阶乘呢？
a[0]:=1;
for i:=1 to n do
begin
   a[i]:=0;
   for j:=n downto 1 do
   begin
      for k:=1 to j do
         a[k]:=a[k]+a[k-1]
      write(a[i],' ');
   end;
end;
    我在网上搜索相关的东西时找到了另一个有趣的东西。对一个初始时全为1的数列反复进行这两个操作：累加求前缀和，然后以1,2,3,…的间隔划掉其中一部分数（即划去所有位置编号为三角形数的数）形成新的序列。类似的数列操作方法最先由Alfred Moessner提出的，我们这里不妨把它叫做Moessner数列。你会发现，第n轮操作开始前，数列的第一个数恰好是n! 。看看下面的例子吧：

1 1 1 1 1 1 1 1 1  1  1  1  1  1  1 …
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 …
x 2 x 4 5 x 7 8 9  x 11 12 13 14  x …

2 4  5  7  8  9 11 12 13 14 …
2 6 11 18 26 35 46 58 71 85 …
x 6  x 18 26  x 46 58 71  x …

6 18 26 46  58  71 …
6 24 50 96 154 225 …
x 24  x 96 154   x …

24  96 154 …
24 120 274 …
x 120  x  …

120 …
…..

    当然，发现前面O(n^3)的程序和这个Moessner数列的关联时我很是吃了一惊：在前面的程序里，如果你输出每一次i循环末所得到的数列，你会发现输出的这些数正好就是后面这个问题里被我们划掉的数，而它们其实就是第一类Stirling数！
    这到底是为什么呢？是什么东西把阶乘、第一类Stirling数、Moessner数列和那个O(n^3)的程序联系在一起的呢？昨天，我想这个问题想了一天，最后终于想通了。如果把Moessner数列排列成这个样子，一切就恍然大悟了：

  
    仔细观察上图，我们会发现：
    1. 按照Moessner数列的定义，每个数都应该等于它左边的数和左上角的数的和（这个“左边”可以跳过若干空格）。例如，35 = 9 + 26，46 = 11 + 35。排成一系列三角形后，每个三角形最右边一列的数就是被划去的数，它永远不能参与它下面的那些行的运算。
    2. 设a[n,i,j]表示左起第n个三角形阵列中的第i行右起第j列上的数，则a[n,i,j]=a[n-1,i-1,j]*n + a[n-1,i,j]，例如274=50*5+24。如果递推时遇到空白位置而它左边隔若干空格的地方还有数的话，则需要用左边的数来补，例如18=4*4+2。对于每个三角形的最后一列来说，这个性质实际上就是第一类Stirling数的递推关系，因此Moessner数列中才会出现第一类Stirling数。
    3. 在第一类Stirling数中，s(n,1)=n! ，也即左起第n个三角形最底端的那个数等于n!。从上面的第二个性质来看，这也是显然的。
    4. O(n^3)的算法实际上就是在绘制上面这个图。每一次j循环末，我们得到
的序列是第i个三角形中每一行左起第j个数组成的序列。例如，计算第5个三角形内的数时，程序首先累加出1, 11, 46, 96, 120, 120，这样便算出了a[5]=120，数列的前5个数再次累加即得到1, 12, 58, 154, 274，由此算出a[4]=274。
    第二个性质可以利用第一个性质进行数学归纳法证明，证明很简单，我就不多说了。现在我尽可能少写一些繁琐的细节，节约一些时间用来复习古代汉语。

做人要厚道，
转贴请注明出处。

查看更多：
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
http://www.luschny.de/math/factorial/index.html <—- 巨牛，20多种阶乘算法的代码！

Google面试题中的两道趣题

/ 44 条评论

看看下面两道题，它的解答非常简单，即使没学过信息学的人也可以想到答案。你能在多短的时间内想出问题的算法来？一小时？一分钟？一秒钟？

1. 给你一个长度为N的链表。N很大，但你不知道N有多大。你的任务是从这N个元素中随机取出k个元素。你只能遍历这个链表一次。你的算法必须保证取出的元素恰好有k个，且它们是完全随机的（出现概率均等）。

2. 给你一个数组A[1..n]，请你在O(n)的时间里构造一个新的数组B[1..n]，使得B[i]=A[1]*A[2]*…*A[n]/A[i]。你不能使用除法运算。

Solution:
1. 遍历链表，给每个元素赋一个0到1之间的随机数作为权重（像Treap一样），最后取出权重最大的k个元素。你可以用一个堆来维护当前最大的k个数。
2. 从前往后扫一遍，然后从后往前再扫一遍。也就是说，线性时间构造两个新数组，P[i]=A[1]*A[2]*…*A[i]，Q[i]=A[n]*A[n-1]*…*A[i]。于是，B[i]=P[i-1]*Q[i+1]。

突然想到，给别人(MM?)介绍信息学时，用这两道题当例子挺合适的。

Retargeting：经典动态规划问题的一个出人意料的应用

/ 18 条评论

    计算机的图片处理技术已经越来越牛B了，很多看似不可能的东西现在都已经有了比较成熟的算法，比如无锯齿放大、抠图、智能抹除等等。但你相信吗，现在竟然有这样一种算法，它可以改变图片的长宽比，同时保持画面内容的长宽比不变！
    我们经常遇到这样一个问题：源图片的长宽比与目标长宽不合，把图片剪裁一部分会觉得可惜，而拉伸图片后画面内容就变形了。此时，一种叫做Retargeting的技术或许可以帮助你：当图片长宽比改变后，它能压缩图片中不重要的部分，保持画面主体内容长宽比不变，让人看不出这个图片是被拉伸过的。
    算法的原理来自这篇名为Seam Carving for Content – Aware Image Resizing的论文，它第一次出现在今年8月份举行的SIGGRAPH大会上。下面一段视频简单地介绍了这个算法，看后你会发现其实质非常简单。

YouTube链接：http://www.youtube.com/watch?v=qadw0BRKeMk

    你会在视频中听到一个OIer特别熟悉的词。“从最上面一排的某个像素出发，每次只能向左下、右下和正下方三个方向移动，求出到达最底端的路径中权值和最小的一条”，这是每个OIer学习动态规划的必修课程，它甚至还出现在了前几天的某次NOIp模拟赛中（记得好像是第三题）。
    最近cnBeta的一篇文章提到，这个技术已经用于实践，一个叫做rsizr的Flash网站可以实现上述算法，感兴趣的同学可以去试试。

07年NOIp模拟赛by Matrix67 解题报告

/ 14 条评论

Problem 1: Matrix67的情书（二）
    大家都知道，看一个数是否能被2整除只需要看它的个位能否被2整除即可。可是你想过为什么吗？这是因为10能被2整除，因此一个数10a+b能被2整除当且仅当b能被2整除。大家也知道，看一个数能否被3整除只需要看各位数之和是否能被3整除。这又是为什么呢？答案或多或少有些类似：因为10^n-1总能被3整除。2345可以写成2*(999+1) + 3*(99+1) + 4*(9+1) + 5，展开就是2*999+3*99+4*9 + 2+3+4+5。前面带了数字9的项肯定都能被3整除了，于是要看2345能否被3整除就只需要看2+3+4+5能否被3整除了。当然，这种技巧只能在10进制下使用，不过类似的结论可以推广到任意进制。
    注意到36是4的整数倍，而ZZZ…ZZ除以7总是得555…55。也就是说，判断一个36进制数能否被4整除只需要看它的个位，而一个36进制数能被7整除当且仅当各位数之和能被7整除。如果一个数同时能被4和7整除，那么这个数就一定能被28整除。于是问题转化为，有多少个连续句子满足各位数字和是7的倍数，同时最后一个数是4的倍数。这样，我们得到了一个O(n)的算法：用P[i]表示前若干个句子除以7的余数为i有多少种情况，扫描整篇文章并不断更新P数组。当某句话的最后一个字能被4整除时，假设以这句话结尾的前缀和除以7余x，则将此时P[x]的值累加到最后的输出结果中（两个前缀的数字和除以7余数相同，则较长的前缀多出来的部分一定整除7）。
    上述算法是我出这道题的本意，但比赛后我见到了其它各种各样新奇的算法。比如有人注意到36^n mod 28总是等于8，利用这个性质也可以构造出类似的线性算法来。还有人用动态规划（或者说递推）完美地解决了这个问题。我们用f[i,j]表示以句子i结束，除以28余数为j的文本片段有多少个；处理下一句话时我们需要对每一个不同的j进行一次扫描，把f[i-1,j]加进对应的f[i,j']中。最后输出所有的f[i,0]的总和即可。这个动态规划可以用滚动数组，因此它的空间同前面的算法一样也是常数的。
    如果你完全不知道我在说什么，你可以看看和进位制、同余相关的文章。另外，我之前还曾出过一道很类似的题(VOJ1090)，你可以对比着看一看。
    另外，非常抱歉地告诉大家，这道题的最后一个数据是错的。这个数据的第一个句子是一个空句子（感谢Ai.Freedom的提醒）。很多第一题只得了90分的好同志估计都是由于我的失误造成的，这里再次表示歉意。如果去掉最前面的那个问号，输出应该是19511110。

Problem 2: 送给MM的生日礼物
        
    我们用f[i,j,k]表示以第i行第j列的格子为右下角，边长为k的正方形是否符合要求。要想f[i,j,k]=True，首先必须满足f[i-1,j-1,k-2]为True（灰色部分满足要求）。另外，我们还需要保证图中两块蓝色区域相等，两块绿色区域相等，并且这四个区域自身还必须是对称的。由于动态规划的状态已经是三方的了，因此判断后面的这些条件必须在常数时间里完成。为此，我们可以在动态规划前进行以下6个预处理：

以同列不同行的两个格子(i,j)和(i',j)为右端点，同时向左扩展可以得到多长的相等区域
以同行不同列的两个格子(i,j)和(i,j')为最底端，同时向上扩展可以得到多长的相等区域
以格子(i,j)为中心，向左右扩展可以得到多长的对称区域
以格子(i,j)为中心，向上下扩展可以得到多长的对称区域
以两个横向相邻的格子(i,j-1)和(i,j)为中心，向左右扩展可以得到多长的对称区域
以两个纵向相邻的格子(i-1,j)和(i,j)为中心，向上下扩展可以得到多长的对称区域

    上面的每个预处理都可以在三方的时间里完成，动态规划的决策降到常数级别，因此总的复杂度还是三方。

    同样地，这也只是我出这道题的本意。我事先已经想到，这道题应该还有很多其它的算法，只是没想到会有那么多。一个比较容易想到的算法是，执行与上面相同的预处理后，枚举某个格子（或某个交叉点）为中心，向外一层一层地进行扩展。虽然这样的复杂度仍然是三方的，并且与前面的算法实质相同，但它的效率显然更高，因为你可以在无法再向外扩展时停止最里面的那个循环，继续枚举下一个中心点。
    同样是枚举正方形的中心点，只使用上面的后4个预处理也可以解决这个问题。我们可以在线性的时间内找出，以某个格子（或交叉点）为中心的最大的合法正方形。假如这个最大的正方形边长为k，这相当于我们同时找到了(k+1) div 2个正方形来。至于如何找到这个最大的正方形，还是留给大家思考吧。
    Cockhorse想到了一个非常巧妙的算法，预处理结束后它可以在常数时间内判断任一个给定的小正方形是否满足题目要求。用f[i,j,k]表示，以第i行中(i,j)及其右边相邻的k-1个格子（共k个格子）来作为底边，所能得到的左右对称的矩形其最大高度是多少。则当f[i,j+1,k-2]为True且(i,j)格与(i,j+k-1)格花色相等时f[i,j,k]=f[i-1,j,k]+1（否则f[i,j,k]=0）。同样地，再用g[i,j,k]表示以(j,i)..(j+k-1,i)作为右边界，使矩形上下对称的最大宽度。这样，判断任意一个正方形是否满足题目要求，只需要检查它的底边和右边能够产生的最大对称区域是否可以覆盖该正方形即可。
    当然，这道题目还有很多其它的算法，这里就不一一列举了。

Problem 3: 流言的传播
    给定一个图，把图中所有点构成的点集叫做U，另指定一个不等于全集的子集S，那么所有一个顶点在S里另一个顶点在U-S里的边就叫做S点集的“割边”，因为去掉所有这样的边后S集将孤立出来。这道题就是需要你找出一个边集E，使得不管S集是什么，对应的割边中权值最小的那一条一定在边集E中。这样的边集肯定是存在的，比如所有边构成的集合就是一个满足条件的边集。这道题希望你找到的边集E里所含的边尽可能的少。
    一个错误的贪心法是，去掉每个点的邻接边中权值最小的那一条。这种算法显然不对，比如下面就是一个鲜活的反例：

o—–o—–o—–o
   2     3     1

    在上图中，每个点的邻接边中权值最小的不是1就是2，这样的话中间那条边就被保留了下来，于是令S集为左边两个点（或右边两个点），割边仍然一条不少。很容易想到，要想让任意S集的割边中权值最小的被去掉，首先得保证边集E连通了所有的点，否则令S等于任一连通分量，边集E里都不会包含任何一条割边。这样，边集E至少要有n-1条边。
&n
bsp;   考虑到最优解至少是n-1条边，这n-1条边必须连通所有的点，而所选边的权值都应该很小，于是想到最后的答案很可能就是最小生成树。现在假设我们已经选择了某些边，这些边形成了若干个连通分量。考虑所有连接两个不同的连通分量的边（两顶点不在同一连通分量的边）中权值最小的一条，那这条边必须要选，否则令S集为其中一个连通分量，题目条件就不能达到了。这不就是Kruskal算法吗？
    且慢，我们还没有证明，对任意一个S集，最小生成树都符合条件。证明其实很简单，假设权值最小的割边e不在最小生成树E中，添加割边e将形成一个回路。这条回路将从S集的某个点出发，经过e跑到U-S里，最后必须还要回到S集。回到S集必然要经过另一条割边e'，而显然e'>e（因为e是权值最小的割边，且题目说了没有权值相同的边），于是边集E+{e}-{e'}就成了更小的生成树。

Problem 4: 表白机器人
    第四题是整个模拟赛中最简单的题，它不需要你构造任何算法，你只需要按照题目意思进行模拟即可。和去年的NOIp一样，纯粹考编程能力的题目并没有放在第一道题的位置上。这提示大家拿到题目后要先看完所有的题，特别是看到最后一道题时千万别慌，很可能把它的衣服扒光了一看，发现它居然是一道赤裸裸的送分题。
    为了加快速度，先预处理出一个数组wall[i][j][k]，表示第i行第j列往k(1<=k<=4)方向上走是否走得通。然后枚举所有可能的命令序列，模拟机器人的行走过程，看它是否能到达终点。在模拟过程中，你需要用一个数组hash[i][j][k]表示机器人曾在序列的第k个命令后到达第i行第j列的位置，并在模拟过程中不断更新hash数组；当某一次命令后机器人还没走到右上角，而对应的hash值已经为True了，则机器人的行动将从这里开始循环，永远也到不了右上角了。
    虽然这道题不需要任何剪枝，但我还想再说一句。这道题有一个非常有趣的剪枝：命令序列的第一个指令只能是U或R，最后一个指令也只能是U或R。大家想想这是为什么。