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

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

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)

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

二进制中的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位整数并输

位运算简介及实用技巧(一):基础篇

    去年年底写的关于位运算的日志是这个Blog里少数大受欢迎的文章之一,很多人都希望我能不断完善那篇文章。后来我看到了不少其它的资料,学习到了更多关于位运算的知识,有了重新整理位运算技巧的想法。从今天起我就开始写这一系列位运算讲解文章,与其说是原来那篇文章的follow-up,不如说是一个remake。当然首先我还是从最基础的东西说起。

什么是位运算?
    程序中的所有数在计算机内存中都是以二进制的形式储存的。位运算说穿了,就是直接对整数在内存中的二进制位进行操作。比如,and运算本来是一个逻辑运算符,但整数与整数之间也可以进行and运算。举个例子,6的二进制是110,11的二进制是1011,那么6 and 11的结果就是2,它是二进制对应位进行逻辑运算的结果(0表示False,1表示True,空位都当0处理):
     110
AND 1011
———-
    0010  –>  2

    由于位运算直接对内存数据进行操作,不需要转成十进制,因此处理速度非常快。当然有人会说,这个快了有什么用,计算6 and 11没有什么实际意义啊。这一系列的文章就将告诉你,位运算到底可以干什么,有些什么经典应用,以及如何用位运算优化你的程序。

Pascal和C中的位运算符号
    下面的a和b都是整数类型,则:
C语言  |  Pascal语言
——-+————-
a & b  |  a and b
a | b  |  a or b
a ^ b  |  a xor b
  ~a   |   not a
a << b |  a shl b
a >> b |  a shr b

    注意C中的逻辑运算和位运算符号是不同的。520|1314=1834,但520||1314=1,因为逻辑运算时520和1314都相当于True。同样的,!a和~a也是有区别的。

各种位运算的使用
    === 1. and运算 ===
    and运算通常用于二进制取位操作,例如一个数 and 1的结果就是取二进制的最末位。这可以用来判断一个整数的奇偶,二进制的最末位为0表示该数为偶数,最末位为1表示该数为奇数.

    === 2. or运算 ===
    or运算通常用于二进制特定位上的无条件赋值,例如一个数or 1的结果就是把二进制最末位强行变成1。如果需要把二进制最末位变成0,对这个数or 1之后再减一就可以了,其实际意义就是把这个数强行变成最接近的偶数。

    === 3. xor运算 ===
    xor运算通常用于对二进制的特定一位进行取反操作,因为异或可以这样定义:0和1异或0都不变,异或1则取反。
    xor运算的逆运算是它本身,也就是说两次异或同一个数最后结果不变,即(a xor b) xor b = a。xor运算可以用于简单的加密,比如我想对我MM说1314520,但怕别人知道,于是双方约定拿我的生日19880516作为密钥。1314520 xor 19880516 = 20665500,我就把20665500告诉MM。MM再次计算20665500 xor 19880516的值,得到1314520,于是她就明白了我的企图。
    下面我们看另外一个东西。定义两个符号#和@(我怎么找不到那个圈里有个叉的字符),这两个符号互为逆运算,也就是说(x # y) @ y = x。现在依次执行下面三条命令,结果是什么?
x <- x # y
y <- x @ y
x <- x @ y

    执行了第一句后x变成了x # y。那么第二句实质就是y <- x # y @ y,由于#和@互为逆运算,那么此时的y变成了原来的x。第三句中x实际上被赋值为(x # y) @ x,如果#运算具有交换律,那么赋值后x就变成最初的y了。这三句话的结果是,x和y的位置互换了。
    加法和减法互为逆运算,并且加法满足交换律。把#换成+,把@换成-,我们可以写出一个不需要临时变量的swap过程(Pascal)。
procedure swap(var a,b:longint);
begin
   a:=a + b;
   b:=a - b;
   a:=a - b;
end;

    好了,刚才不是说xor的逆运算是它本身吗?于是我们就有了一个看起来非常诡异的swap过程:
procedure swap(var a,b:longint);
begin
   a:=a xor b;
   b:=a xor b;
   a:=a xor b;
end;

    === 4. not运算 ===
    not运算的定义是把内存中的0和1全部取反。使用not运算时要格外小心,你需要注意整数类型有没有符号。如果not的对象是无符号整数(不能表示负数),那么得到的值就是它与该类型上界的差,因为无符号类型的数是用$0000到$FFFF依次表示的。下面的两个程序(仅语言不同)均返回65435。
var
   a:word;
begin
   a:=100;
   a:=not a;
   writeln(a);
end.

#include <stdio.h>
int main()
{
    unsigned short a=100;
    a = ~a;
    printf( "%dn", a );    
    return 0;
}

    如果not的对象是有符号的整数,情况就不一样了,稍后我们会在“整数类型的储存”小节中提到。

    === 5. shl运算 ===
    a shl b就表示把a转为二进制后左移b位(在后面添b个0)。例如100的二进制为1100100,而110010000转成十进制是400,那么100 shl 2 = 400。可以看出,a shl b的值实际上就是a乘以2的b次方,因为在二进制数后添一个0就相当于该数乘以2。
    通常认为a shl 1比a * 2更快,因为前者是更底层一些的操作。因此程序中乘以2的操作请尽量用左移一位来代替。
    定义一些常量可能会用到shl运算。你可以方便地用1 shl 16 – 1来表示65535。很多算法和数据结构要求数据规模必须是2的幂,此时可以用shl来定义Max_N等常量。

    === 6. shr运算 ===
    和shl相似,a shr b表示二进制右移b位(去掉末b位),相当于a除以2的b次方(取整)。我们也经常用shr 1来代替div 2,比如二分查找、堆的插入操作等等。想办法用shr代替除法运算可以使程序效率大大提高。最大公约数的二进制算法用除以2操作来代替慢得出奇的mod运算,效率可以提高60%。

位运算的简单应用
    有时我们的程序需要一个规模不大的Hash表来记录状态。比如,做数独时我们需要27个Hash表来统计每一行、每一列和每一个小九宫格里已经有哪些数了。此时,我们可以用27个小于2^9的整数进行记录。例如,一个只填了2和5的小九宫格就用数字18表示(二进制为000010010),而某一行的状态为511则表示这一行已经填满。需要改变状态时我们不需要把这个数转成二进制修改后再转回去,而是直接进行位操作。在搜索时,把状态表示成整数可以更好地进行判重等操作。这道题是在搜索中使用位运算加速的经典例子。以后我们会看到更多的例子。
    下面列举了一些常见的二进制位的变换操作。

    功能              |           示例            |    位运算
———————-+—————————+——————–
去掉最后一位          | (101101->10110)           | x shr 1
在最后加一个0         | (101101->1011010)         | x shl 1
在最后加一个1         | (101101->1011011)         | x shl 1+1
把最后一位变成1       | (101100->101101)          | x or 1
把最后一位变成0       | (101101->101100)          | x or 1-1
最后一位取反          | (101101->101100)          | x xor 1
把右数第k位变成1      | (101001->101101,k=3)      | x or (1 shl (k-1))
把右数第k位变成0      | (101101->101001,k=3)      | x and not (1 shl (k-1))
右数第k位取反         | (101001->101101,k=3)      | x xor (1 shl (k-1))
取末三位              | (1101101->101)            | x and 7
取末k位               | (1101101->1101,k=5)       | x and (1 shl k-1)
取右数第k位           | (1101101->1,k=4)          | x shr (k-1) and 1
把末k位变成1          | (101001->101111,k=4)      | x or (1 shl k-1)
末k位取反             | (101001->100110,k=4)      | x xor (1 shl k-1)
把右边连续的1变成0    | (100101111->100100000)    | x and (x+1)
把右起第一个0变成1    | (100101111->100111111)    | x or (x+1)
把右边连续的0变成1    | (11011000->11011111)      | x or (x-1)
取右边连续的1         | (100101111->1111)         | (x xor (x+1)) shr 1
去掉右起第一个1的左边 | (100101000->1000)         | x and (x xor (x-1))

    最后这一个在树状数组中会用到。

Pascal和C中的16进制表示
    Pascal中需要在16进制数前加$符号表示,C中需要在前面加0x来表示。这个以后我们会经常用到。

整数类型的储存
    我们前面所说的位运算都没有涉及负数,都假设这些运算是在unsigned/word类型(只能表示正数的整型)上进行操作。但计算机如何处理有正负符号的整数类型呢?下面两个程序都是考察16位整数的储存方式(只是语言不同)。
var
   a,b:integer;
begin
   a:=$0000;
   b:=$0001;
   write(a,' ',b,' ');
   a:=$FFFE;
   b:=$FFFF;
   write(a,' ',b,' ');
   a:=$7FFF;
   b:=$8000;
   writeln(a,' ',b);
end.

#include <stdio.h>
int main()
{
    short int a, b;
    a = 0x0000;
    b = 0x0001;
    printf( "%d %d ", a, b );
    a = 0xFFFE;
    b = 0xFFFF;
    printf( "%d %d ", a, b );
    a = 0x7FFF;
    b = 0x8000;
    printf( "%d %dn", a, b );
    return 0;
}

    两个程序的输出均为0 1 -2 -1 32767 -32768。其中前两个数是内存值最小的时候,中间两个数则是内存值最大的时候,最后输出的两个数是正数与负数的分界处。由此你可以清楚地看到计算机是如何储存一个整数的:计算机用$0000到$7FFF依次表示0到32767的数,剩下的$8000到$FFFF依次表示-32768到-1的数。32位有符号整数的储存方式也是类似的。稍加注意你会发现,二进制的第一位是用来表示正负号的,0表示正,1表示负。这里有一个问题:0本来既不是正数,也不是负数,但它占用了$0000的位置,因此有符号的整数类型范围中正数个数比负数少一个。对一个有符号的数进行not运算后,最高位的变化将导致正负颠倒,并且数的绝对值会差1。也就是说,not a实际上等于-a-1。这种整数储存方式叫做“补码”。

最后还有两句话
    Matrix67原创
    转贴请注明出处

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

    一个数是素数(也叫质数),当且仅当它的约数只有两个——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原创
转贴请注明出处

Matrix67生日邀请赛 标程公布

之所以没公布标程,是因为个人觉得标程写得比较丑。
既然有人需要就发布一下吧,标程丑总比没有标程好。

Problem 1
program whyleast;

procedure solve(t,a,b:integer);
begin
   if t=0 then exit else
   begin
      solve(t-1,a,b);
      writeln(a,' ',2);
      solve(t-1,b,a);
      writeln(2,' ',b);
      solve(t-1,a,b);
   end;
end;

{====main====}
var
   n,i:integer;
   ans:longint=1;
begin
  assign(input,'whyleast.in');
  reset(input);
  assign(output,'whyleast.out');
  rewrite(output);
  
  readln(n);
  for i:=1 to n do ans:=ans*3;
  writeln(ans-1);
  solve(n,1,3);
  
  close(input);
  close(output);
end.

Problem 2
program height;

const
   OutputString:array[boolean]of string=('YES','NO');

type
   rec1=record
          h,p:longint;
        end;

   pointer=^rec2;
   rec2=record
          x,y:longint;
          dir:boolean;
          next:pointer;
        end;
var
   orderHeight : array[1..100000]of longint;
   SortHeight  : array[1..100000]of rec1;
   Deg,DegHash : array[0..100000]of longint;
   EdgeHash    : array[1..100000]of pointer;
   n,m,DegCount:longint;

procedure SwapRec(var t1,t2:rec1);
var
   t3:rec1;
begin
   t3:=t1;
   t1:=t2;
   t2:=t3;
end;

procedure SwapInt(var t1,t2:longint);
var
   t3:longint;
begin
   t3:=t1;
   t1:=t2;
   t2:=t3;
end;

function InsertEdgeHash(x,y:longint):integer;
var
   newp:pointer;
begin
   new(newp);
   newp^.x:=x;
   newp^.y:=y;

   if orderHeight[x] > orderHeight[y] then
      newp^.dir:=( 1.25*OrderHeight[y] <= orderHeight[x] )
   else newp^.dir:=( 1.25*OrderHeight[x] > orderHeight[y] );
   newp^.dir:=not newp^.dir;

   newp^.next:=EdgeHash[x];
   EdgeHash[x]:=newp;
   exit( ord( newp^.dir ) );
end;

function FindEdgeHash(x,y:longint):integer;
         { -1: Not Found;  0: x-->y  1: x<--y
            x Always Smaller than y }
var
   now:pointer;
begin
   now:=EdgeHash[x];
   while now<>nil do
   begin
      if now^.y=y then
      begin
         now^.dir:=not now^.dir;
         exit( ord( now^.dir ) );
      end;
      now:=now^.next;
   end;
   exit(-1);
end;

procedure UpdateDeg(t,c:longint);
begin
   if deg[t]<>c then
   begin
     dec(DegHash[Deg[t]]);
     if DegHash[Deg[t]]=0 then dec(DegCount);
     inc(DegHash[c]);
     if DegHash[c]=1 then inc(DegCount);
     Deg[t]:=c;
   end;
end;

procedure ReadHeight;
var
   i:longint;
begin
   readln(n,m);
   for i:=1 to n do
   begin
      readln(OrderHeight[i]);
      SortHeight[i].h:=OrderHeight[i];
      SortHeight[i].p:=i;
   end;
end;

procedure Sort(l,r:longint);
var
   i,j,mid:longint;
begin
   i:=l;
   j:=r;
   mid:=SortHeight[(i+j)div 2].h;
   repeat
      while SortHeight[i].h<mid do inc(i);
      while SortHeight[j].h>mid do dec(j);
      if i<=j then
      begin
         SwapRec(SortHeight[i],SortHeight[j]);
         inc(i);
         dec(j);
      end;
   until i>j;
   if l<j then Sort(l,j);
   if i<r then Sort(i,r);
end;

procedure Init;
var
   low:longint=1;
   high:longint=1;
   i:longint;
begin
   DegHash[0]:=n;
   DegCount:=1;
   for i:=1 to n do
   begin
      while SortHeight[low].h*1.25 < SortHeight[i].h do inc(low);
      while ( high<n ) and ( SortHeight[i].h*1.25 >= SortHeight[high+1].h ) do inc(high);
      UpdateDeg( SortHeight[i].p, high+low-i );
   end;
end;

procedure Solve;
var
   i,x,y:longint;
   dir:integer;
   newp:pointer;
begin
   for i:=1 to m do
   begin
      readln(x,y);
      if x>y then SwapInt(x,y);
      dir:=FindEdgeHash(x,y);
      if dir=-1 then dir:=InsertEdgeHash(x,y);
      if dir=0 then SwapInt(x,y);
      UpdateDeg(x,Deg[x]+1);
      UpdateDeg(y,Deg[y]-1);
      writeln( OutputString[DegCount=n] );
   end;
end;

{====main====}
begin
   assign(input,'height.in');
   reset(input);
   assign(output,'height.out');
   rewrite(output);

   ReadHeight;
   Sort(1,n);
   Init;
   Solve;

   close(input);
   close(output);
end.

Problem 3
program wolf;

type
   rec=record
          left,right:integer;
       end;

const
   infinite=maxlongint div 3-100000;
   //Make sure no overflows occur

var
   k,n,m  : integer;
   map    : array[1..1000,1..1000]of longint;
   dist   : array[1..1000]of longint;
   hash   : array[1..1000]of boolean;
   father : array[1..1000]of longint;

 
;  tree : array[1..1000]of rec;
   attk : array[1..1000]of longint;
   cost : array[1..1000]of integer;
   minf : array[1..1000,0..100]of longint;

procedure readp;
var
   i,x,y,d:longint;
begin
   readln(k,n,m);
   for i:=2 to n do
      readln(attk[i],cost[i]);
   for i:=1 to m do
   begin
      readln(x,y,d);
      map[x,y]:=d;
      map[y,x]:=d;
   end;
end;

procedure init;
var
   i,j:longint;
begin
   for i:=2 to n do dist[i]:=infinite;
   for i:=2 to n do hash[i]:=false;
   dist[1]:=0;
   hash[1]:=true;

   for i:=1 to n do
   for j:=1 to n do
      if map[i,j]=0 then map[i,j]:=infinite;

   for i:=1 to n do
   for j:=1 to k do
      minf[i,j]:=-infinite;
end;

procedure sssp;
var
   i,j:longint;
   min:longint=0;
   minj:longint=1;
begin
   for i:=1 to n-1 do
   begin
      for j:=1 to n do if not hash[j] then
      begin
         if ( min+map[minj,j] = dist[j] ) and ( father[j] > minj ) then
           father[j]:=minj
         else if min+map[minj,j] < dist[j] then
         begin
           dist[j]:=min + map[minj,j];
           father[j]:=minj;
         end;
      end;

      min:=infinite;
      for j:=1 to n do if not hash[j] and (dist[j]<min) then
      begin
         minj:=j;
         min:=dist[j];
      end;

      tree[ minj ].right:=tree[ father[minj] ].left;
      tree[ father[minj] ].left:=minj;
      hash[ minj ]:=true;
   end;
end;

function solve(x,y:longint):longint;  //(node,cost)

   procedure update(var t1:longint;t2:longint);
   begin
      if t1<t2 then t1:=t2;
   end;

var
   ans:longint=-infinite;
   i,t:longint;
begin
   if minf[x,y]<>-infinite then exit(minf[x,y]);
   if y>=cost[x] then ans:=attk[x];

   if tree[x].left>0 then update(ans,solve(tree[x].left,y)+attk[x]);
  
   if tree[x].right>0 then
   begin
      update(ans,solve(tree[x].right,y));
      if y-cost[x]>0 then
         update(ans,solve(tree[x].right,y-cost[x])+attk[x]);
   end;
  
   if (tree[x].left>0) and (tree[x].right>0) then
      for i:=1 to y-1 do
         update(ans,solve(tree[x].left,i)+solve(tree[x].right,y-i)+attk[x]);

   minf[x,y]:=ans;
   exit(minf[x,y]);
end;

procedure writep;
var
   ans:longint=-infinite;
   i,j:integer;
begin
   for i:=0 to k do
     if solve(1,i)>ans then ans:=solve(1,i);
   writeln(ans);
  
   {===For Debug===
   for i:=1 to n do
   begin
      for j:=1 to k do write(minf[i,j]:3);
      writeln;
   end;
   for i:=1 to n do writeln(tree[i].left,' ',tree[i].right);
   }
end;

{====main====}
begin
   assign(input,'wolf.in');
   reset(input);
   assign(output,'wolf.out');
   rewrite(output);

   readp;
   init;
   sssp;
   writep;

   close(input);
   close(output);
end.

Problem 4
program garden;

const
   dir:array[1..4,1..2]of integer=
     ((1,0),(0,1),(-1,0),(0,-1));

type
   arr=array[1..10]of integer;
   rec=record x,y:integer;end;

var
   map:array[0..11,0..11]of boolean;
   ans:array[1..100]of rec;
   n,m,max:integer;
   step:integer=1;
   state:arr;

procedure readp;
var
   i,j:integer;
   ch:char;
begin
   readln(m,n);
   for i:=1 to n do
   begin
      for j:=1 to m do
      begin
         read(ch);
         map[i,j]:=(ch='1');
         inc(max,ord( map[i,j] ))
      end;
   readln;
   end;
end;

procedure writep;
var
   i:integer;
begin
   for i:=1 to step do
      writeln( '(' , ans[i].x , ',' , ans[i].y , ')' );
end;

procedure solve(x,y:integer);
var
   tx,ty,d:integer;
   step_cache:integer;
   state_cache:arr;
begin
   step_cache:=step;
   state_cache:=state;
   if step=max then
   begin
      writep;
      exit;
   end;

   for d:=1 to 4 do
   begin
      tx:=x+dir[d,1];
      ty:=y+dir[d,2];
      while map[tx,ty] and ( not state[tx] and(1 shl (ty-1) )>0) do
      begin
         inc(step);
         ans[step].x:=tx;
         ans[step].y:=ty;
         state[tx]:=state[tx] or ( 1 shl (ty-1) );
         tx:=tx+dir[d,1];
         ty:=ty+dir[d,2];
      end;

      tx:=tx-dir[d,1];
      ty:=ty-dir[d,2];
      if (tx<>x) or (ty<>y) then solve(tx,ty);
      state:=state_cache;
      step:=step_cache;
   end;
end;

{====main====}
var
   i,j:integer;
begin
   assign(input,'garden.in');
   reset(input);
   assign(output,'garden.o
ut');
   rewrite(output);

   readp;
   for i:=1 to n do
   for j:=1 to m do
     if map[i,j] then
     begin
        ans[1].x:=i;
        ans[1].y:=j;
        state[i]:=1 shl (j-1);
        solve(i,j);
        state[i]:=0;
     end;

   close(input);
   close(output);
end.

依然欢迎大家来挑错