算法

魔方问题新进展：26步足以破解任何魔方

/ 18 条评论

    最近，波士顿Northeastern大学的计算机科学家Daniel Kunkle证明了任何一个魔方可以在26步以内解开。这个结果打破了以往所有的记录。在解魔方的处理过程中，他构造了一些非常具有启发性的算法，这篇文章将简单地介绍一下这些算法。
    一个魔方大约有4.3 x 10^19种可能的初始状态，再牛的机器也不可能搜索完所有的可能。因此Kunkle和他的指导员Gene Cooperman想出了一些对魔方状态进行分类筛选的策略。
    Kunkle和Cooperman首先运用了一个小技巧将问题进行简化。如果魔方的每个面全是一种颜色，我们就认为魔方被解开，不管哪一面是哪一种颜色。换句话说，相互之间可以通过颜色置换得到的初始状态都是等价的。这样，“本质不同”的初始状态就减少到10^18种。
    接下来，他们开始观察一类更简单的问题：如果只允许180度转动(half-turn)，有多少状态可以被解决。在10^18种状态中，只有大约15000种状态可以仅用180度旋转来破解。对于普通计算机来说，这个数目也不大，只需要不到一天的时间就可以搜索出解开所有15000多个魔方各自需要的最少步数。他们发现，这类初始状态中任何一个都可以在13步以内解决。
    然后他们需要做的就是找出，需要多少步才能把任意一种状态转化为这15000种特殊状态中的一个。为了完成这一工作，首先他们把所有的初始状态划分为若干个等价类，每个等价类里的状态都可以仅用180度转动互相得到。这样，同一个等价类中如果任一状态可以变换为其中一种特殊状态，同样的转动步骤也可以使该等价类的其它所有状态都变成特殊状态。最后他们找到了1.4 x 10^12个不同的等价类，需要解决的状态数由最初的4.3×10^19减少到1.4×10^12。但无论如何，10^12仍然是一个恐怖的数字。
    现在他们用了一台超级计算机来完成这个工作，并且使用了一些很有技巧性的决策来加速搜索过程。计算机需要耗费大量的时间读取硬盘上的数据，为了加快速度，Kunkle和Cooperman将数据巧妙地进行了处理，使得数据的排列正好与计算机读取的顺序相符，这样就节省了搜索硬盘的时间。
    “这种方法可以应用在任何一个组合问题上”，Kunkle说。他提到了西洋跳棋、国际象棋、航班安排和蛋白质摺叠等一系列问题。一种类似的组合学方法最近被用于寻找西洋跳棋的最优策略中。
    63小时的计算后，超级计算机得到的答案是，任何一种状态都能在16步以内转化为15000种特殊状态。而这些特殊状态又只需要13步达到最终状态，因此这种方法最终得到的结论是：29步以内可以解决任何一个魔方问题。
    但这个数字还不足以创造出新的记录，去年瑞典就曾经得到过27步内解决魔方问题的结论。Kunkle和Cooperman意识到，要想打破这个记录，他们还需要削减3步才行。
    应用他们现有的算法，只有8×10^7个状态集合还不能做到26步以内出解。再次对这些相对较少的状态进行搜索，最终他们找到了26步以内解决所有魔方的方法。
    7月29日他们在ISSAC(International Symposium on Symbolic and Algebraic Computation,国际符号和代数计算会议)上公布了这一结果。
    现在Kunkle和Cooperman希望把最大步骤数减少到25。他们认为他们可以对所有需要26步的状态进行暴力搜索来寻找更优的方案。
    虽然他们已经获得了很大的成功，但这一结果很可能还有改进的空间。许多学者认为20步以内足以解决任何魔方，但现在没有人能够证明。

Matrix67翻译，原文地址
做人要厚道，转贴请注明出处

十个利用矩阵乘法解决的经典题目

/ 989 条评论

    好像目前还没有这方面题目的总结。这几天连续看到四个问这类题目的人，今天在这里简单写一下。这里我们不介绍其它有关矩阵的知识，只介绍矩阵乘法和相关性质。
    不要以为数学中的矩阵也是黑色屏幕上不断变化的绿色字符。在数学中，一个矩阵说穿了就是一个二维数组。一个n行m列的矩阵可以乘以一个m行p列的矩阵，得到的结果是一个n行p列的矩阵，其中的第i行第j列位置上的数等于前一个矩阵第i行上的m个数与后一个矩阵第j列上的m个数对应相乘后所有m个乘积的和。比如，下面的算式表示一个2行2列的矩阵乘以2行3列的矩阵，其结果是一个2行3列的矩阵。其中，结果的那个4等于2*2+0*1：
    
    下面的算式则是一个1 x 3的矩阵乘以3 x 2的矩阵，得到一个1 x 2的矩阵：
    

    矩阵乘法的两个重要性质：一，矩阵乘法不满足交换律；二，矩阵乘法满足结合律。为什么矩阵乘法不满足交换律呢？废话，交换过来后两个矩阵有可能根本不能相乘。为什么它又满足结合律呢？仔细想想你会发现这也是废话。假设你有三个矩阵A、B、C，那么(AB)C和A(BC)的结果的第i行第j列上的数都等于所有A(ik)*B(kl)*C(lj)的和（枚举所有的k和l）。

经典题目1 给定n个点，m个操作，构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转
    这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转（两种情况），旋转则以原点为中心。如果对每个点分别进行模拟，那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵，然后每个点与该矩阵相乘即可直接得出最终该点的位置，总共耗时O(m+n)。假设初始时某个点的坐标为x和y，下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来，再乘以(x,y,1)，即可一步得出最终点的位置。
    

经典题目2 给定矩阵A，请快速计算出A^n（n个A相乘）的结果，输出的每个数都mod p。
    由于矩阵乘法具有结合律，因此A^4 = A * A * A * A = (A*A) * (A*A) = A^2 * A^2。我们可以得到这样的结论：当n为偶数时，A^n = A^(n/2) * A^(n/2)；当n为奇数时，A^n = A^(n/2) * A^(n/2) * A （其中n/2取整）。这就告诉我们，计算A^n也可以使用二分快速求幂的方法。例如，为了算出A^25的值，我们只需要递归地计算出A^12、A^6、A^3的值即可。根据这里的一些结果，我们可以在计算过程中不断取模，避免高精度运算。

经典题目3 POJ3233 (感谢rmq)
    题目大意：给定矩阵A，求A + A^2 + A^3 + … + A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据mod m。k<=10^9。
    这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：
    A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)
    应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A + A^2 + A^3，即可得到原问题的答案。

经典题目4 VOJ1049
    题目大意：顺次给出m个置换，反复使用这m个置换对初始序列进行操作，问k次置换后的序列。m<=10, k<2^31。
    首先将这m个置换“合并”起来（算出这m个置换的乘积），然后接下来我们需要执行这个置换k/m次（取整，若有余数则剩下几步模拟即可）。注意任意一个置换都可以表示成矩阵的形式。例如，将1 2 3 4置换为3 1 2 4，相当于下面的矩阵乘法：
    
    置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方，再乘以初始序列即可。做出来了别忙着高兴，得意之时就是你灭亡之日，别忘了最后可能还有几个置换需要模拟。

经典题目5 《算法艺术与信息学竞赛》207页（2.1代数方法和模型，[例题5]细菌，版次不同可能页码有偏差）
    大家自己去看看吧，书上讲得很详细。解题方法和上一题类似，都是用矩阵来表示操作，然后二分求最终状态。

经典题目6 给定n和p，求第n个Fibonacci数mod p的值，n不超过2^31
    根据前面的一些思路，现在我们需要构造一个2 x 2的矩阵，使得它乘以(a,b)得到的结果是(b,a+b)。每多乘一次这个矩阵，这两个数就会多迭代一次。那么，我们把这个2 x 2的矩阵自乘n次，再乘以(0,1)就可以得到第n个Fibonacci数了。不用多想，这个2 x 2的矩阵很容易构造出来：
    

经典题目7 VOJ1067
    我们可以用上面的方法二分求出任何一个线性递推式的第n项，其对应矩阵的构造方法为：在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n行填对应的系数，其它地方都填0。例如，我们可以用下面的矩阵乘法来二分计算f(n) = 4f(n-1) – 3f(n-2) + 2f(n-4)的第k项：
    
    利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。这里给出的例题是系数全为1的情况。

经典题目8 给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数mod p的值
    把给定的图转为邻接矩阵，即A(i,j)=1当且仅当存在一条边i->j。令C=A*A，那么C(i,j)=ΣA(i,k)*A(k,j)，实际上就等于从点i到点j恰好经过2条边的路径数（枚举k为中转点）。类似地，C*A的第i行第j列就表示从i到j经过3条边的路径数。同理，如果要求经过k步的路径数，我们只需要二分求出A^k即可。

经典题目9 用1 x 2的多米诺骨牌填满M x N的矩形有多少种方案，M<=5，N<2^31，输出答案mod p的结果
    
    我们以M=3为例进行讲解。假设我们把这个矩形横着放在电脑屏幕上，从右往左一列一列地进行填充。其中前n-2列已经填满了，第n-1列参差不齐。现在我们要做的事情是把第n-1列也填满，将状态转

位运算简介及实用技巧（三）：进阶篇(2)

/ 68 条评论

今天我们来看两个稍微复杂一点的例子。

n皇后问题位运算版
    n皇后问题是啥我就不说了吧，学编程的肯定都见过。下面的十多行代码是n皇后问题的一个高效位运算程序，看到过的人都夸它牛。初始时，upperlim:=(1 shl n)-1。主程序调用test(0,0,0)后sum的值就是n皇后总的解数。拿这个去交USACO，0.3s，暴爽。
procedure test(row,ld,rd:longint);
var
      pos,p:longint;
begin

{ 1}  if row<>upperlim then
{ 2}  begin
{ 3}     pos:=upperlim and not (row or ld or rd);
{ 4}     while pos<>0 do
{ 5}     begin
{ 6}        p:=pos and -pos;
{ 7}        pos:=pos-p;
{ 8}        test(row+p,(ld+p)shl 1,(rd+p)shr 1);
{ 9}     end;
{10}  end
{11}  else inc(sum);

end;
    乍一看似乎完全摸不着头脑，实际上整个程序是非常容易理解的。这里还是建议大家自己单步运行一探究竟，实在没研究出来再看下面的解说。

  
    和普通算法一样，这是一个递归过程，程序一行一行地寻找可以放皇后的地方。过程带三个参数，row、ld和rd，分别表示在纵列和两个对角线方向的限制条件下这一行的哪些地方不能放。我们以6×6的棋盘为例，看看程序是怎么工作的。假设现在已经递归到第四层，前三层放的子已经标在左图上了。红色、蓝色和绿色的线分别表示三个方向上有冲突的位置，位于该行上的冲突位置就用row、ld和rd中的1来表示。把它们三个并起来，得到该行所有的禁位，取反后就得到所有可以放的位置（用pos来表示）。前面说过-a相当于not a + 1，这里的代码第6行就相当于pos and (not pos + 1)，其结果是取出最右边的那个1。这样，p就表示该行的某个可以放子的位置，把它从pos中移除并递归调用test过程。注意递归调用时三个参数的变化，每个参数都加上了一个禁位，但两个对角线方向的禁位对下一行的影响需要平移一位。最后，如果递归到某个时候发现row=111111了，说明六个皇后全放进去了，此时程序从第1行跳到第11行，找到的解的个数加一。

    ~~~~====~~~~=====   华丽的分割线   =====~~~~====~~~~

Gray码
    假如我有4个潜在的GF，我需要决定最终到底和谁在一起。一个简单的办法就是，依次和每个MM交往一段时间，最后选择给我带来的“满意度”最大的MM。但看了dd牛的理论后，事情开始变得复杂了：我可以选择和多个MM在一起。这样，需要考核的状态变成了2^4=16种（当然包括0000这一状态，因为我有可能是玻璃）。现在的问题就是，我应该用什么顺序来遍历这16种状态呢？
    传统的做法是，用二进制数的顺序来遍历所有可能的组合。也就是说，我需要以0000->0001->0010->0011->0100->…->1111这样的顺序对每种状态进行测试。这个顺序很不科学，很多时候状态的转移都很耗时。比如从0111到1000时我需要暂时甩掉当前所有的3个MM，然后去把第4个MM。同时改变所有MM与我的关系是一件何等巨大的工程啊。因此，我希望知道，是否有一种方法可以使得，从没有MM这一状态出发，每次只改变我和一个MM的关系（追或者甩），15次操作后恰好遍历完所有可能的组合（最终状态不一定是1111）。大家自己先试一试看行不行。
    解决这个问题的方法很巧妙。我们来说明，假如我们已经知道了n=2时的合法遍历顺序，我们如何得到n=3的遍历顺序。显然，n=2的遍历顺序如下：

00
01
11
10

    你可能已经想到了如何把上面的遍历顺序扩展到n=3的情况。n=3时一共有8种状态，其中前面4个把n=2的遍历顺序照搬下来，然后把它们对称翻折下去并在最前面加上1作为后面4个状态：

000
001
011
010  ↑
——–
110  ↓
111
101
100

    用这种方法得到的遍历顺序显然符合要求。首先，上面8个状态恰好是n=3时的所有8种组合，因为它们是在n=2的全部四种组合的基础上考虑选不选第3个元素所得到的。然后我们看到，后面一半的状态应该和前面一半一样满足“相邻状态间仅一位不同”的限制，而“镜面”处则是最前面那一位数不同。再次翻折三阶遍历顺序，我们就得到了刚才的问题的答案：

0000
0001
0011
0010
0110
0111
0101
0100
1100
1101
1111
1110
1010
1011
1001
1000

    这种遍历顺序作为一种编码方式存在，叫做Gray码（写个中文让蜘蛛来抓：格雷码）。它的应用范围很广。比如，n阶的Gray码相当于在n维立方体上的Hamilton回路，因为沿着立方体上的边走一步，n维坐标中只会有一个值改变。再比如，Gray码和Hanoi塔问题等价。Gray码改变的是第几个数，Hanoi塔就该移动哪个盘子。比如，3阶的Gray码每次改变的元素所在位置依次为1-2-1-3-1-2-1，这正好是3阶Hanoi塔每次移动盘子编号。如果我们可以快速求出Gray码的第n个数是多少，我们就可以输出任意步数后Hanoi塔的移动步骤。现在我告诉你，Gray码的第n个数（从0算起）是n xor (n shr 1)，你能想出来这是为什么吗？先自己想想吧。

    下面我们把二进制数和Gray码都写在下面，可以看到左边的数异或自身右移的结果就等于右边的数。

二进制数   Gray码
   000       000
   001       001
   010       011
   011       010
   100       110
   101       111
   110       101
   111       100

    从二进制数的角度看，“镜像”位置上的数即是对原数进行not运算后的结果。比如，第3个数010和倒数第3个数101的每一位都正好相反。假设这两个数分别为x和y，那么x xor (x shr 1)和y xor (y shr 1)的结果只有一点不同：后者的首位是1，前者的首位是0。而这正好是Gray码的生成方法。这就说明了，Gray码的第n个数确实是n xor (n shr 1)。

&nbsp
;   今年四月份mashuo给我看了这道题，是二维意义上的Gray码。题目大意是说，把0到2^(n+m)-1的数写成2^n * 2^m的矩阵，使得位置相邻两数的二进制表示只有一位之差。答案其实很简单，所有数都是由m位的Gray码和n位Gray码拼接而成，需要用左移操作和or运算完成。完整的代码如下：
var
   x,y,m,n,u:longint;
begin
   readln(m,n);
   for x:=0 to 1 shl m-1 do begin
      u:=(x xor (x shr 1)) shl n; //输出数的左边是一个m位的Gray码
      for y:=0 to 1 shl n-1 do
         write(u or (y xor (y shr 1)),' '); //并上一个n位Gray码
      writeln;
   end;
end.

Matrix67原创
转贴请注明出处

位运算简介及实用技巧（二）：进阶篇(1)

/ 66 条评论

=====   真正强的东西来了！   =====

二进制中的1有奇数个还是偶数个
    我们可以用下面的代码来计算一个32位整数的二进制中1的个数的奇偶性，当输入数据的二进制表示里有偶数个数字1时程序输出0，有奇数个则输出1。例如，1314520的二进制101000000111011011000中有9个1，则x=1314520时程序输出1。
var
   i,x,c:longint;
begin
   readln(x);
   c:=0;
   for i:=1 to 32 do
   begin
      c:=c + x and 1;
      x:=x shr 1;
   end;
   writeln( c and 1 );
end.
    但这样的效率并不高，位运算的神奇之处还没有体现出来。
    同样是判断二进制中1的个数的奇偶性，下面这段代码就强了。你能看出这个代码的原理吗？
var
   x:longint;
begin
   readln(x);
   x:=x xor (x shr 1);
   x:=x xor (x shr 2);
   x:=x xor (x shr 4);
   x:=x xor (x shr 8);
   x:=x xor (x shr 16);
   writeln(x and 1);
end.
    为了说明上面这段代码的原理，我们还是拿1314520出来说事。1314520的二进制为101000000111011011000，第一次异或操作的结果如下：

    00000000000101000000111011011000
XOR  0000000000010100000011101101100
—————————————
    00000000000111100000100110110100

    得到的结果是一个新的二进制数，其中右起第i位上的数表示原数中第i和i+1位上有奇数个1还是偶数个1。比如，最右边那个0表示原数末两位有偶数个1，右起第3位上的1就表示原数的这个位置和前一个位置中有奇数个1。对这个数进行第二次异或的结果如下：

    00000000000111100000100110110100
XOR   000000000001111000001001101101
—————————————
    00000000000110011000101111011001

    结果里的每个1表示原数的该位置及其前面三个位置中共有奇数个1，每个0就表示原数对应的四个位置上共偶数个1。一直做到第五次异或结束后，得到的二进制数的最末位就表示整个32位数里有多少个1，这就是我们最终想要的答案。

计算二进制中的1的个数
    同样假设x是一个32位整数。经过下面五次赋值后，x的值就是原数的二进制表示中数字1的个数。比如，初始时x为1314520（网友抓狂：能不能换一个数啊），那么最后x就变成了9，它表示1314520的二进制中有9个1。
x := (x and $55555555) + ((x shr 1) and $55555555);
x := (x and $33333333) + ((x shr 2) and $33333333);
x := (x and $0F0F0F0F) + ((x shr 4) and $0F0F0F0F);
x := (x and $00FF00FF) + ((x shr 8) and $00FF00FF);
x := (x and $0000FFFF) + ((x shr 16) and $0000FFFF);
    为了便于解说，我们下面仅说明这个程序是如何对一个8位整数进行处理的。我们拿数字211（我们班某MM的生日）来开刀。211的二进制为11010011。

+—+—+—+—+—+—+—+—+
| 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1 |   <—原数
+—+—+—+—+—+—+—+—+
|  1 0  |  0 1  |  0 0  |  1 0  |   <—第一次运算后
+——-+——-+——-+——-+
|    0 0 1 1    |    0 0 1 0    |   <—第二次运算后
+—————+—————+
|        0 0 0 0 0 1 0 1        |   <—第三次运算后，得数为5
+——————————-+

    整个程序是一个分治的思想。第一次我们把每相邻的两位加起来，得到每两位里1的个数，比如前两位10就表示原数的前两位有2个1。第二次我们继续两两相加，10+01=11，00+10=10，得到的结果是00110010，它表示原数前4位有3个1，末4位有2个1。最后一次我们把0011和0010加起来，得到的就是整个二进制中1的个数。程序中巧妙地使用取位和右移，比如第二行中$33333333的二进制为00110011001100….，用它和x做and运算就相当于以2为单位间隔取数。shr的作用就是让加法运算的相同数位对齐。

二分查找32位整数的前导0个数
    这里用的C语言，我直接Copy的Hacker's Delight上的代码。这段代码写成C要好看些，写成Pascal的话会出现很多begin和end，搞得代码很难看。程序思想是二分查找，应该很简单，我就不细说了。
int nlz(unsigned x)
{
   int n;

   if (x == 0) return(32);
   n = 1;
   if ((x >> 16) == 0) {n = n +16; x = x <<16;}
   if ((x >> 24) == 0) {n = n + 8; x = x << 8;}
   if ((x >> 28) == 0) {n = n + 4; x = x << 4;}
   if ((x >> 30) == 0) {n = n + 2; x = x << 2;}
   n = n - (x >> 31);
   return n;
}

只用位运算来取绝对值
    这是一个非常有趣的问题。大家先自己想想吧，Ctrl+A显示答案。
    答案：假设x为32位整数，则x xor (not (x shr 31) + 1) + x shr 31的结果是x的绝对值
    x shr 31是二进制的最高位，它用来表示x的符号。如果它为0（x为正），则not (x shr 31) + 1等于$00000000，异或任何数结果都不变；如果最高位为1（x为负），则not (x shr 31) + 1等于$FFFFFFFF，x异或它相当于所有数位取反，异或完后再加一。

高低位交换
    这个题实际上是我出的，做为学校内部NOIp模拟赛的第一题。题目是这样：

    给出一个小于2^32的正整数。这个数可以用一个32位的二进制数表示（不足32位用0补足）。我们称这个二进制数的前16位为“高位”，后16位为“低位”。将它的高低位交换，我们可以得到一个新的数。试问这个新的数是多少（用十进制表示）。
　　例如，数1314520用二进制表示为0000 0000 0001 0100 0000 1110 1101 1000（添加了11个前导0补足为32位），其中前16位为高位，即0000 0000 0001 0100；后16位为低位，即0000 1110 1101 1000。将它的高低位进行交换，我们得到了一个新的二进制数0000 1110 1101 1000 0000 0000 0001 0100。它即是十进制的249036820。

    当时几乎没有人想到用一句位操作来代替冗长的程序。使用位运算的话两句话就完了。
var
   n:dword;
begin
   readln( n );
   writeln( (n shr 16) or (n  shl 16) );
end.
    而事实上，Pascal有一个系统函数swap直接就可以用。

二进制逆序
    下面的程序读入一个32位整数并输

数论部分第一节：素数与素性测试

/ 75 条评论

    一个数是素数（也叫质数），当且仅当它的约数只有两个——1和它本身。规定这两个约数不能相同，因此1不是素数。对素数的研究属于数论范畴，你可以看到许多数学家没事就想出一些符合某种性质的素数并称它为某某某素数。整个数论几乎就围绕着整除和素数之类的词转过去转过来。对于写代码的人来说，素数比想像中的更重要，Google一下BigPrime或者big_prime你总会发现大堆大堆用到了素数常量的程序代码。平时没事时可以记一些素数下来以备急用。我会选一些好记的素数，比如4567, 124567, 3214567, 23456789, 55566677, 1234567894987654321, 11111111111111111111111 (23个1)。我的手机号前10位是个素数。我的网站域名的ASCII码连起来(77 97 116 114 105 120 54 55 46 99 111 109)也是个素数。还有，我的某个MM的八位生日也是一个素数。每次写Hash函数之类的东西需要一个BigPrime常量时我就取她的生日，希望她能给我带来好运。偶尔我叫她素MM，没人知道是啥意思，她自己也不知道。
    素数有很多神奇的性质。我写5个在下面供大家欣赏。

1. 素数的个数无限多（不存在最大的素数）
  证明：反证法，假设存在最大的素数P，那么我们可以构造一个新的数2 * 3 * 5 * 7 * … * P + 1（所有的素数乘起来加1）。显然这个数不能被任一素数整除（所有素数除它都余1），这说明我们找到了一个更大的素数。

2. 存在任意长的一段连续数，其中的所有数都是合数（相邻素数之间的间隔任意大）
  证明：当0<a<=n时，n!+a能被a整除。长度为n-1的数列n!+2, n!+3, n!+4, …, n!+n中，所有的数都是合数。这个结论对所有大于1的整数n都成立，而n可以取到任意大。

3. 所有大于2的素数都可以唯一地表示成两个平方数之差。
  证明：大于2的素数都是奇数。假设这个数是2n+1。由于(n+1)^2=n^2+2n+1，(n+1)^2和n^2就是我们要找的两个平方数。下面证明这个方案是唯一的。如果素数p能表示成a^2-b^2，则p=a^2-b^2=(a+b)(a-b)。由于p是素数，那么只可能a+b=p且a-b=1，这给出了a和b的唯一解。

4. 当n为大于2的整数时，2^n+1和2^n-1两个数中，如果其中一个数是素数，那么另一个数一定是合数。
  证明：2^n不能被3整除。如果它被3除余1，那么2^n-1就能被3整除；如果被3除余2，那么2^n+1就能被3整除。总之，2^n+1和2^n-1中至少有一个是合数。

5. 如果p是素数，a是小于p的正整数，那么a^(p-1) mod p = 1。
  这个证明就有点麻烦了。
    首先我们证明这样一个结论：如果p是一个素数的话，那么对任意一个小于p的正整数a，a, 2a, 3a, …, (p-1)a除以p的余数正好是一个1到p-1的排列。例如，5是素数，3, 6, 9, 12除以5的余数分别为3, 1, 4, 2，正好就是1到4这四个数。
    反证法，假如结论不成立的话，那么就是说有两个小于p的正整数m和n使得na和ma除以p的余数相同。不妨假设n>m，则p可以整除a(n-m)。但p是素数，那么a和n-m中至少有一个含有因子p。这显然是不可能的，因为a和n-m都比p小。
    用同余式表述，我们证明了：
(p-1)! ≡ a * 2a * 3a * … * (p-1)a (mod p)
    也即：
(p-1)! ≡ (p-1)! * a^(p-1) (mod p)
    两边同时除以(p-1)!，就得到了我们的最终结论：
1 ≡ a^(p-1) (mod p)

    可惜最后这个定理最初不是我证明的。这是大数学家Fermat证明的，叫做Fermat小定理(Fermat's Little Theorem)。Euler对这个定理进行了推广，叫做Euler定理。Euler一生的定理太多了，为了和其它的“Euler定理”区别开来，有些地方叫做Fermat小定理的Euler推广。Euler定理中需要用一个函数f(m)，它表示小于m的正整数中有多少个数和m互素（两个数只有公约数1称为互素）。为了方便，我们通常用记号φ(m)来表示这个函数（称作Euler函数）。Euler指出，如果a和m互素，那么a^φ(m) ≡ 1 (mod m)。可以看到，当m为素数时，φ(m)就等于m-1（所有小于m的正整数都与m互素），因此它是Fermat小定理的推广。定理的证明和Fermat小定理几乎相同，只是要考虑的式子变成了所有与m互素的数的乘积：m_1 * m_2 … m_φ(m) ≡ (a * m_1)(a * m_2) … (a * m_φ(m)) (mod m)。我为什么要顺便说一下Euler定理呢？因为下面一句话可以增加我网站的PV：这个定理出现在了The Hundred Greatest Theorems里。

    谈到Fermat小定理，数学历史上有很多误解。很长一段时间里，人们都认为Fermat小定理的逆命题是正确的，并且有人亲自验证了a=2, p<300的所有情况。国外甚至流传着一种说法，认为中国在孔子时代就证明了这样的定理：如果n整除2^(n-1)-1，则n就是素数。后来某个英国学者进行考证后才发现那是因为他们翻译中国古文时出了错。1819年有人发现了Fermat小定理逆命题的第一个反例：虽然2的340次方除以341余1，但341=11*31。后来，人们又发现了561, 645, 1105等数都表明a=2时Fermat小定理的逆命题不成立。虽然这样的数不多，但不能忽视它们的存在。于是，人们把所有能整除2^(n-1)-1的合数n叫做伪素数(pseudoprime)，意思就是告诉人们这个素数是假的。
    不满足2^(n-1) mod n = 1的n一定不是素数；如果满足的话则多半是素数。这样，一个比试除法效率更高的素性判断方法出现了：制作一张伪素数表，记录某个范围内的所有伪素数，那么所有满足2^(n-1) mod n = 1且不在伪素数表中的n就是素数。之所以这种方法更快，是因为我们可以使用二分法快速计算2^(n-1) mod n 的值，这在计算机的帮助下变得非常容易；在计算机中也可以用二分查找有序数列、Hash表开散列、构建Trie树等方法使得查找伪素数表效率更高。
    有人自然会关心这样一个问题：伪素数的个数到底有多少？换句话说，如果我只计算2^(n-1) mod n的值，事先不准备伪素数表，那么素性判断出错的概率有多少？研究这个问题是很有价值的，毕竟我们是OIer，不可能背一个长度上千的常量数组带上考场。统计表明，在前10亿个自然数中共有50847534个素数，而满足2^(n-1) mod n = 1的合数n有5597个。这样算下来，算法出错的可能性约为0.00011。这个概率太高了，如果想免去建立伪素数表的工作，我们需要改进素性判断的算法。

    最简单的想法就是，我们刚才只考虑了a=2的情况。对于式子a^(n-1) mod n，取不同的a可能导致不同的结果。一个合数可能在a=2时通过了测试，但a=3时的计算结果却排除了素数的可能。于是，人们扩展了伪素数的定义，称满足a^(n-1) mod n = 1的合数n叫做以a为底的伪素数(pseudoprime to base a)。前10亿个自然数中同时以2和3为底的伪素数只有1272个，这个数目不到刚才的1/4。这告诉我们如果同时验证a=2和a=3两种情况，算法出错的概率降到了0.000025。容易想到，选择用来测试的a越多，算法越准确。通常我们的做法是，随机选择若干个小于待测数的正整数作为底数a进行若干次测试，只要有一次没有通过测试就立即把
这个数扔回合数的世界。这就是Fermat素性测试。
    人们自然会想，如果考虑了所有小于n的底数a，出错的概率是否就可以降到0呢？没想到的是，居然就有这样的合数，它可以通过所有a的测试（这个说法不准确，详见我在地核楼层的回复）。Carmichael第一个发现这样极端的伪素数，他把它们称作Carmichael数。你一定会以为这样的数一定很大。错。第一个Carmichael数小得惊人，仅仅是一个三位数，561。前10亿个自然数中Carmichael数也有600个之多。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

    Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。新的测试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时x=p-1)。
    我们下面来演示一下上面的定理如何应用在Fermat素性测试上。前面说过341可以通过以2为底的Fermat测试，因为2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，面具后面真实的嘴脸显现了出来，想假扮素数和我的素MM交往的企图暴露了出来。
    这就是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。那么我们需要计算的东西就变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足a^(d * 2^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。这样，Fermat小定理加强为如下形式：
    尽可能提取因子2，把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) （注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）
    Miller-Rabin素性测试同样是不确定算法，我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。
    Miller-Rabin算法的代码也非常简单：计算d和r的值（可以用位运算加速），然后二分计算a^d mod n的值，最后把它平方r次。程序的代码比想像中的更简单，我写一份放在下边。虽然我已经转C了，但我相信还有很多人看不懂C语言。我再写一次Pascal吧。函数IsPrime返回对于特定的底数a，n是否是能通过测试。如果函数返回False，那说明n不是素数；如果函数返回True，那么n极有可能是素数。注意这个代码的数据范围限制在longint，你很可能需要把它们改成int64或高精度计算。
function pow( a, d, n:longint ):longint;
begin
   if d=0 then exit(1)
   else if d=1 then exit(a)
   else if d and 1=0 then exit( pow( a*a mod n, d div 2, n) mod n)
   else exit( (pow( a*a mod n, d div 2, n) * a) mod n);
end;

function IsPrime( a,n:longint ):boolean;
var
   d,t:longint;
begin
   if n=2 then exit(true);
   if (n=1) or (n and 1=0) then exit(false);
   d:=n-1;
   while d and 1=0 do d:=d shr 1;
   t:=pow( a, d, n );
   while ( d<>n-1 ) and ( t<>1 ) and ( t<>n-1 ) do
   begin
      t:=(t * t)mod n;
      d:=d shl 1;
   end;
   exit( (t=n-1) or (d and 1=1) );
end;
    对于大数的素性判断，目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取，但当待测数不太大时，选择测试底数就有一些技巧了。比如，如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足够了。当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在OI中非常实用。通常认为，Miller-Rabin素性测试的正确率可以令人接受，随机选取k个底数进行测试算法的失误率大概为4^(-k)。

    Miller-Rabin算法是一个RP算法。RP是时间复杂度的一种，主要针对判定性问题。一个算法是RP算法表明它可以在多项式的时间里完成，对于答案为否定的情形能够准确做出判断，但同时它也有可能把对的判成错的（错误概率不能超过1/2）。RP算法是基于随机化的，因此多次运行该算法可以降低错误率。还有其它的素性测试算法也是概率型的，比如Solovay-Strassen算法。另外一些素性测试算法则需要预先知道一些辅助信息（比如n-1的质因子），或者需要待测数满足一些条件（比如待测数必须是2^n-1的形式）。前几年AKS算法轰动世界，它是第一个多项式的、确定的、无需其它条件的素性判断算法。当时一篇论文发表出来，题目就叫PRIMES is in P，然后整个世界都疯了，我们班有几个MM那天还来了初潮。算法主要基于下面的事实：n是一个素数当且仅当(x-a)^n≡(x^n-a) (mod n)。注意这个x是多项式中的未知数，等式两边各是一个多项式。举个例子来说，当a=1时命题等价于如下结论：当n是素数时，杨辉三角的第n+1行除两头的1以外其它的数都能被n整除。

Matrix67原创
转贴请注明出处