f(Program,Poet)=Programet

其实这是一份比较无聊的东西……因为民意表决掉了期中考试，于是C语言老师znm便换了法子，要我们写一个实验报告出来，其中重点是结构化程序设计。本来嘛，搜索类题目是可以很好地表现出结构化程序设计的优点的（比如递归……），但是我决定写另外一个题目，而把搜索类题目推荐给其他同学~（比如八皇后啊，骑士游历啊这种网上能找到一沓代码的题目）。我选择了素性测试这个问题。当然不需要太难，为了体现结构化，我觉得将不同的方法分别写一个函数然后比较一下就可以了。在下面贴出来的文字中不免有很多装13的地方，同时估计会有说得不严谨的地方…砖不要拍得太狠了…以及有些东西我很懒，就没有真正地完成所有的函数，同时找了一个理由。另：请尽量不要转载，免得老师以为我抄袭——即便我在报告末尾声明过了……

关于检验素数的算法的探究

素数是数论中最重要的概念，同时在实际中也有着广泛的应用。比如RSA公钥加密系统就利用了大数不好分解成素因数的乘积这一特点——其中就包括了如何判定一个数字是素数这一问题。那么我们就来探究一下如何判定素数。

首先应该想到一个最最简单的方法：即对于整数n，如果我们发现(n%i!=0)对于i属于[2,n-1]时恒成立，那么就可以说n是一个素数了，即用小于n大于1的整数来除n看看能否整除。这个方法实现起来也很简单，效率是O(n)。对应代码中的函数 int method1(int n){}。

但是难道不能优化么？很容易发现，n-1永远都不会整除n——其实，大于[n/2]（此处[x]为高斯取整函数）小于n的整数都不能整除n。那么显然，我们只需要查看2到[n/2]的整数来除n得到的结果就行了。从时间上来说比第一种方法快了一倍，仍然是O(n)。对应代码中的函数int method2(int n){}。

又仔细想了一下，会发现其实第二种方法还是有一点点缺陷。比如我判断25是不是素数的时候，发现2无法除尽25，那么此时所有的偶数都无法除尽25，发现3也无法除尽25，那么所有的3的倍数都无法除尽25。所以，这时发现只要有一个数无法除尽n，则这个数的倍数都无法除尽n。反过来，如果一个数能除尽n，那么它的各个因数也能够除尽。此时不难得到，如果所有小于[n/2]的素数都不能整除n，那么n肯定是一个素数。对这个思路进一步提炼，就可以知道因为素数是递增的，所以用反证法易得当i>sqrt(n)仍然不能整除n时，那么n就是一个素数——此时极大地减少了需要判断的次数。由素数定理可知，小于整数x的素数个数约为x/lnx，则需要列举的素数约等于2*sqrt(n)/ln(n)，算法效率为O(sqrt(n)/ln(n))，效率在此时被极大提高了。不过这时就牵涉到了一个问题：如何找到所有小于n的素数？待会再来讨论这个问题。

检验素数的方法除了确定性的，还有不确定性的——由费马小定理开创的素性检验的概率算法。费马小定理是说假如a是一个整数，p是一个素数，那么。那么，反过来，如果a是一个整数同时有上式成立，那么p是不是素数呢？很遗憾，答案是否定的。但是我们可以通过更换a的值来进行多次检验，减小判断错误的概率。理论上只需要将a取遍所有小于p的素数即可给出答案，但是这样会让效率变得很低。这个方法的缺点有两个：一是计算a^p会比较缓慢，二是即使优化了乘方的写法也会很有可能需要为了非高精度的a和p专门进行高精度处理，减慢了速度。对应代码中的int Fermat(int n,int t){}，其中调用了另外一个求幂次的函数 int power(int a,int p){}。

另外一个著名的素性检验的概率算法就是Miller-Rabin算法。要测试N是否为质数，首先将N − 1分解为d*2^s。在每次测试开始时，先随机选一个介于 [1,N − 1]的整数 a，之后如果对所有的，若且，则N是合数。否则，N有3/4的机率为质数。Miller-Rabin检验的好处在于你可以估计出错的概率是多少，我重复m次那么最后出错的概率就是1/(4^m)。如果将费马小定理和MR检验结合起来，那么准确率将会大大提高。

由于之前说到了使用“小于n的所有素数”之类的方法，那么我们就应该来考虑一下如何求出一张素数表。最容易想到的方法就是先检验出2到n-1的素数。但是这样的话效率就变成了O(n^2)。换一种思路：所有素数的倍数都不是素数，那么好办啦，我们每发现一个素数p，就将2到n-1中的p的倍数给找出来，标记为合数，而这个过程中始终没有被标记过的数自然就是素数了。这个方法其实称为埃拉托斯特尼筛法，是古希腊数学家埃拉托斯特尼所提出的。其实可以知道，我们只需要筛到[n/2]即可。更进一步的，如果使用一个数组prime[i]存储得到的第i个素数，那么和前面类似的理由，我们筛到sqrt(n)就可以了。对应代码中的int Eratosthenes1(int n){}。

接下来大概可以给出两种改进算法。第一种是这样的：考虑到2和3的倍数都不是合数，那么显然，素数只可能是6n±1的形式。我们循环的时候可以直接考虑所有的6n±1形式的数字。其实这样以来还可以再考虑5，7，11，……但这样就没有一个终止的时候了。6n±1比较简洁同时没有增加太多中间运算步骤，可以接受。对应代码中的int Eratosthenes2(int n){}。

另一种优化方法则巧妙很多。考虑到最原始的筛法中6会被2和3筛两次，60就会被2，3，5筛3次。那么我们就想办法让每一个合数都被筛一次。结合代码加以说明。

void makeprime() { for (i=0;i<=n;i++) isprime[i]=1; len=0; isprime[1]=0; for (i=2;i<=n;i++) { if (isprime[i]) prime[len++]=i; for (j=0;j<=n;j++)//标记一 { isprime[prime[j]*i]=0; if (i%prime[j]==0)//标记二 break; } } }

以上便是生成素数的函数的主体代码。

利用了每个合数必有一个最小素因子。每个合数仅被它的最小素因子筛去正好一次。所以为线性时间。代码中体现在：

if(i%prime[j]==0) break;

prime[]数组中的素数是递增的,当i能整除prime[j]，那么i*prime[j+1]这个合数肯定被prime[j]乘以某个数筛掉。因为i中含有prime[j],prime[j]比prime[j+1]小，即i=k*prime[j]，那么i*prime[j+1]=(k*prime[j])*prime[j+1]=k’*prime[j]，接下去的素数同理。所以不用筛下去了。因此，在满足i%prime[j]==0这个条件之前以及第一次满足改条件时,prime[j]必定是prime[j]*i的最小因子。这个方法改良后的时间效率是O(n)！int Eratosthenes3(int n){}，其中因为是一旦判定了n是合数那么就退出，如果设n的最小质因数是p，那么时间效率应该是O(p+n/p)。

但是到这里还没完，我在维基百科上找到了一个更厉害的算法，它的时间效率是。这个算法的名字叫做Sieve of Atkin，大概能够翻译为“阿特金筛法”，是埃氏筛法的极大改良版。我还没能证明这个算法，只能将步骤写下了：

首先定义r(n)=n%60，再定义反转筛中某数的状态就是让它从是素数到不是素数或相反。

预处理：素数表中有2，3，5三个元素，其余的元素一直在筛子中。

• 如果r(n) 是 1, 13, 17, 29, 37, 41, 49 或53之一，反转方程4x² + y² = n所有的正整数解的状态；
• 如果r(n) 是 7, 19, 31 或 43 之一，反转方程3x² + y² = n. 所有的正整数解的状态；
• 如果 r(n) 是 11, 23, 47 或 59 之一，反转方程{3x² − y² = n;x > y.} 所有的正整数解的状态；
• 对于其他的r(n)可以不予理睬。

将筛子中的最小的标记为“是素数”的数字添加到素数表中，同时将其平方以及平方的倍数都标记为非素数。

经过以上的步骤，就可以得到n以内的素数表。但是由于与O(n)实际上相差无几，所以这个方法的理论意义大于实际意义。维基这个词条的英文链接：http://en.wikipedia.org/wiki/Sieve_of_atkin 。我上面的文字来自那些文字的我的翻译。

在实际的测试中，由于数据范围较小，所以说每一种算法的时间都是一瞬间。不过在实际的应用中，还需要结合高精度运算，滚动数组等方法扩展数据的计算与存储范围。如果前文所讨论的问题要想运用到实际的RSA等应用方面，算法的效率以及复杂而高性能的编码就更为重要了。

通过上文的分析可以发现，想要优化一个算法，那么最重要的就是减少冗余运算（比如第一到第三种方法的优化），或者是做足够多的预处理(生成素数表)，甚至是放弃准确度以求得一个概率上的解。这两种思想可以非常广泛地应用到搜索剪枝、动态规划、记忆化搜索等更高级的方法中去。因此，这个关于检验素数的算法的讨论更重要的是优化算法的思想以及对时间复杂度、实际代码复杂度的权衡。

以上是算法分析，下面是对代码的简要说明：

1. 所有的求解函数中都以返回1表示所求的n是素数，返回0则表示不是素数。
2. 每一个函数都记录了这个函数运行始末使用的时间。
3. 代码中为了让各个函数之间的停顿不至于太短，因此使用了#include<windows.h>，调用了其中的Sleep函数，让程序暂停运行。
4. 费马小定理无法在简单的编码下发挥实际的检验作用，因此只是作为一个概率算法的展示，Miller-Rabin因为同样的原因，便不再进行编码。
5. SieveOfAtkin在编码上较为困难，且我暂时不明白其机理，因此不准备进行编码，以免程序运行结果不良好的时候难以查找错误。
6. 代码基本上由我独立完成，能够成功编译并成功运行。
7. 由于数据范围所限，各个算法之间的比较还不能凸现出来。因此其实还可以对代码进行如下变换进行比较：即使用判定某一个数是否为素数的程序来生成素数表，即可大致看出其与专门生成素数表的程序的效率高低。或者还有另外一个更简单的思路：将每一个程序运行10000遍，也可以得到比较明显的时间差别。

本文可能会在实验报告提交之后放到我的博客上去：http://blog.programet.cn/，所以如果在网上的其他地方出现，日期不可能早于12月11日——而且应该都是转载。

这里说一下，O(n)的筛法来自于以前我在这里看到的代码。

与此报告配套的代码请点击下载，因为老师要求使用VS编译，所以请大家自行更改一些细节后再使用gcc编译。

本文由 严酷的魔王 创作，转载或引用前请联系我们。

相关文章：

1. 谜之系列(0):线性规划
2. 统计与推理(2)：置信区间
3. 无穷中的二分（一）
4. Euler筛法
5. 从(0,1]×(0,1]到(0,1]的双射

2009年12月11日 星期五