位换记号、排列测试与状态图:杂耍中的数学

2016 年 7 月 30 日至 8 月 7 日,第 39 届欧洲杂耍大会(European Juggling Convention)在荷兰的阿尔梅勒举行, 8 月 3 日凌晨的搏击之夜(Fight Night)自然再度成为了众人关注的焦点——它是杂耍斗(combat juggling)这项运动最大的赛事。在杂耍斗的圈子里,有两个响当当的大名你必须要知道:德国选手 Jochen Pfeiffer 目前世界排名第二,之前拿过 6 次搏击之夜的冠军;英国选手 Luke Burrage 目前世界排名第一,之前拿过 8 次搏击之夜的冠军。这一年的比赛中,两位老将均以完胜的成绩轻松进入 32 强,并在淘汰赛阶段过关斩将,最终成功在决赛场上相遇。最终,世界排名第二的 Jochen 以 5 比 4 的成绩击败了世界排名第一的 Luke ,夺得了又一个搏击之夜的冠军。

杂耍斗是一种两人对战类的体育运动。比赛规则非常简单。每局比赛开始时,两名选手各自抛耍 3 个杂耍棒。任何一方都可以故意上前干扰另一方(但只能针对对方手中的或者空中的杂耍棒,不能针对对方的手臂和身体)。谁站到最后,谁就赢得该局。先赢 5 局者获得比赛的胜利。

典型的一局比赛大致就像下面这样。这是 Jochen 和 Luke 的第 6 局比赛。

Read more…

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

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。大家想想这是为什么。

同余运算及其基本性质

    100除以7的余数是2,意思就是说把100个东西七个七个分成一组的话最后还剩2个。余数有一个严格的定义:假如被除数是a,除数是b(假设它们均为正整数),那么我们总能够找到一个小于b的自然数r和一个整数m,使得a=bm+r。这个r就是a除以b的余数,m被称作商。我们经常用mod来表示取余,a除以b余r就写成a mod b = r。
    如果两个数a和b之差能被m整除,那么我们就说a和b对模数m同余(关于m同余)。比如,100-60除以8正好除尽,我们就说100和60对于模数8同余。它的另一层含义就是说,100和60除以8的余数相同。a和b对m同余,我们记作a≡b(mod m)。比如,刚才的例子可以写成100≡60(mod 8)。你会发现这种记号到处都在用,比如和数论相关的书中就经常把a mod 3 = 1写作a≡1(mod 3)。
    之所以把同余当作一种运算,是因为同余满足运算的诸多性质。比如,同余满足等价关系。具体地说,它满足自反性(一个数永远和自己同余)、对称性(a和b同余,b和a也就同余)和传递性(a和b同余,b和c同余可以推出a和c同余)。这三个性质都是显然的。
    同余运算里还有稍微复杂一些的性质。比如,同余运算和整数加减法一样满足“等量加等量,其和不变”。小学我们就知道,等式两边可以同时加上一个相等的数。例如,a=b可以推出a+100=b+100。这样的性质在同余运算中也有:对于同一个模数m,如果a和b同余,x和y同余,那么a+x和b+y也同余。在我看来,这个结论几乎是显然的。当然,我们也可以严格证明这个定理。这个定理对减法同样有效。

    性质:如果a≡b(mod m),x≡y(mod m),则a+x≡b+y(mod m)。
    证明:条件告诉我们,可以找到p和q使得a-mp = b-mq,也存在r和s使得x-mr = y-ms。于是a-mp + x-mr = b-mq + y-ms,即a+x-m(p+r) = b+y-m(q+s),这就告诉我们a+x和b+y除以m的余数相同。

    容易想到,两个同余式对应相乘,同余式两边仍然相等:
    如果a≡b(mod m),x≡y(mod m),则ax≡by(mod m)。
    证明:条件告诉我们,a-mp = b-mq,x-mr = y-ms。于是(a-mp)(x-mr) = (b-mq)(y-ms),等式两边分别展开后必然是ax-m(…) = by-m(…)的形式,这就说明ax≡by(mod m)。

    现在你知道为什么有的题要叫你“输出答案mod xxxxx的结果”了吧,那是为了避免高精度运算,因为这里的结论告诉我们在运算过程中边算边mod和算完后再mod的结果一样。假如a是一个很大的数,令b=a mod m,那么(a * 100) mod m和(b * 100) mod m的结果是完全一样的,这相当于是在a≡b (mod m)的两边同时乘以100。这些结论其实都很显然,因为同余运算只关心余数(不关心“整的部分”),完全可以每一次运算后都只保留余数。因此,整个运算过程中参与运算的数都不超过m,避免了高精度的出现。

    在证明Fermat小定理时,我们用到了这样一个定理:
    如果ac≡bc(mod m),且c和m互质,则a≡b(mod m) (就是说同余式两边可以同时除以一个和模数互质的数)。
    证明:条件告诉我们,ac-mp = bc-mq,移项可得ac-bc = mp-mq,也就是说(a-b)c = m(p-q)。这表明,(a-b)c里需要含有因子m,但c和m互质,因此只有可能是a-b被m整除,也即a≡b(mod m)。

    可能以后还要用到更多的定理,到时候在这里更新。

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原创
转贴请注明出处