C语言求素数问题算法.doc

如何求素数1 自然数是0，1，22 素数是2，3，5（不包括1的只能背1和它本身整除的自然数）#include#include void main()int i ,j, flag=1;for(i=101; i200; i+)flag = 1;for(j=2; jNNN。而这是不可能的，所以，d1和d2中必有一个小于或等于N。基于上述分析，设计算法如下：(1)用2，3，5，7逐个试除N的方法求出100以内的所有素数。(2)用100以内的所有素数逐个试除的方法求出10000以内的素数。首先，将2，3，5，7分别存放在a1、a2、a3、a4中，以后每求出一个素数，只要不大于100，就依次存放在A数组中的一个单元 中。当我们求10010000之间的素数时，可依次用a1a2的素数去试除N，这个范围内的素数可以不保存，直接打印。【2】用筛法求素数。简单介绍一下厄拉多塞筛法。厄拉多塞是一位古希腊数学家，他在寻找素数时，采用了一种与众不同的方法：先将2N的各数写在纸上：在2的上面画一个圆圈，然后划去2的其他倍数；第一个既未画圈又没有被划去的数是3，将它画圈，再划去3的其他倍数；现在既未画圈又没有被划去的第一个数 是5，将它画圈，并划去5的其他倍数依次类推，一直到所有小于或等于N的各数都画了圈或划去为止。这时，表中画了圈的以及未划去的那些数正好就是小于 N的素数。这很像一面筛子，把满足条件的数留下来，把不满足条件的数筛掉。由于这种方法是厄拉多塞首先发明的，所以，后人就把这种方法称作厄拉多塞筛法。在计算机中，筛法可以用给数组单元置零的方法来实现。具体来说就是：首先开一个数组：ai，i=1，2，3，同时，令所有的数组元素都等于下标 值，即ai=i，当i不是素数时，令ai=0 。当输出结果时，只要判断ai是否等于零即可，如果ai=0，则令i=i+1，检查下一个ai。筛法是计算机程序设计中常用的算法之一。【3】用6N1法求素数。任何一个自然数，总可以表示成为如下的形式之一：6N，6N+1，6N+2，6N+3，6N+4，6N+5 (N=0，1，2，)显然，当N1时，6N，6N+2，6N+3，6N+4都不是素数，只有形如6N+1和6N+5的自然数有可能是素数。所以，除了2和3之外，所有的素数都可以表示成6N1的形式(N为自然数)。根据上述分析，我们可以构造另一面筛子，只对形如6 N1的自然数进行筛选，这样就可以大大减少筛选的次数，从而进一步提高程序的运行效率和速度。在程序上，我们可以用一个二重循环实现这一点，外循环i按3的倍数递增，内循环j为01的循环，则2(i+j)-1恰好就是形如6N1的自然数。浅析求素数算法 注意: 如果没有特殊说明, 以下讨论的都是针对n为素数时的时间复杂度1. 根据概念判断:如果一个正整数只有两个因子, 1和p，则称p为素数.代码: 时间复杂度O(n).bool isPrime(int n) if(n 2) return false; for(int i = 2; i n; +i) if(n%i = 0) return false; return true;2. 改进, 去掉偶数的判断代码: 时间复杂度O(n/2), 速度提高一倍.bool isPrime(int n) if(n 2) return false; if(n = 2) return true; for(int i = 3; i n; i += 2) if(n%i = 0) return false; return true;3. 进一步减少判断的范围定理: 如果n不是素数, 则n有满足1d=sqrt(n)的一个因子d.证明: 如果n不是素数, 则由定义n有一个因子d满足1dn.如果d大于sqrt(n), 则n/d是满足1n/d=sqrt(n)的一个因子.代码: 时间复杂度O(sqrt(n)/2), 速度提高O(n-sqrt(n)/2).bool isPrime(int n) if(n 2) return false; if(n = 2) return true; for(int i = 3; i*i = n; i += 2)/使用 isqrt(n) if(n%i = 0) return false; return true;4. 剔除因子中的重复判断.?例如: 11%3 != 0 可以确定 11%(3*i) != 0.定理: 如果n不是素数, 则n有满足1d=sqrt(n)的一个素数因子d.证明: I1. 如果n不是素数, 则n有满足1 d sqrt(n)范围内的所有素数bool isPrime(int primes, int n) if(n 2) return false; for(int i = 0; primesi*primesi = n; +i) if(n%primesi = 0) return false; return true;假设n范围内的素数个数为PI(n), 则时间复杂度O(PI(sqrt(n).函数PI(x)满足素数定理: ln(x)-3/2 x/PI(x) = 67时.因此O(PI(sqrt(n)可以表示为O(sqrt(x)/(ln(sqrt(x)-3/2),O(sqrt(x)/(ln(sqrt(x)-3/2)也是这个算法的空间复杂度.5. 构造素数序列primesi: 2, 3, 5, 7, .由4的算法我们知道, 在素数序列已经被构造的情况下, 判断n是否为素数效率很高;但是, 在构造素数序列本身的时候, 是否也可是达到最好的效率呢?事实上这是可以的! - 我们在构造的时候完全可以利用已经被构造的素数序列!假设我们已经我素数序列: p1, p2, . pn现在要判断pn+1是否是素数, 则需要(1, sqrt(pn+1)范围内的所有素数序列,而这个素数序列显然已经作为p1, p2, . pn的一个子集被包含了!代码:/ 构造素数序列primesvoid makePrimes(int primes, int num) int i, j, temp; primes0 = 2; primes1 = 3; for(i = 5, temp = 2; temp num; i += 2) int flag = true; for(j = 1; primesj*primesj = i; +j) if(i%primesj = 0) flag = false; break; if(flag) primestemp+ = i; makePrimes的时间复杂度比较复杂, 而且它只有在初始化的时候才被调用一次.在一定的应用范围内, 我们可以把近似认为makePrimes需要常数时间.在后面的讨论中, 我们将探讨一种对计算机而言更好的makePrimes方法.6. 更好地利用计算机资源.当前的主流PC中, 一个整数的大小为232. 如果需要判断232大小的数是否为素数,则可能需要测试2, 216范围内的所有素数(216 = sqrt(232).由4中提到的素数定理我们可以大概确定2, 216范围内的素数个数.由于216/(ln(216)-1/2) = 6138, 216/(ln(216)-3/2) = 6834,我们可以大概估计出2, 216范围内的素数个数6138 PI(216) 6834.在对2, 216范围内的素数进行统计, 发现只有6542个素数:p_6542: 65521, 655212 = 4293001441 232, (232 = 4294967296)在实际运算时unsigned long x = 4295098369;将发生溢出, 为131073.在程序中, 我是采用double类型计算得到的结果.分析到这里我们可以看到, 我们只需要缓冲6543个素数, 我们就可以采用4中的算法高效率地判断2, 232如此庞大范围内的素数!(原本的232大小的问题规模现在已经被减小到6543规模了!)虽然用现在的计算机处理2, 216范围内的6542个素数已经没有一点问题,虽然makePrimes只要被运行一次就可以, 但是我们还是考虑一下是否被改进的可能?!我想学过java的人肯定想把makePrimes作为一个静态的初始化实现, 在C+中也可以模拟java中静态的初始化的类似实现:#define NELEMS(x) (sizeof(x) / (sizeof(x)0)static int primes6542+1;static struct _Init _Init()makePrimes(primes, NELEMS(primes); _init;如此, 就可以在程序启动的时候自动掉用makePrimes初始化素数序列.但, 我现在的想法是: 为什么我们不能在编译的时候调用makePrimes函数呢?完全可以! 代码如下:/ 这段代码可以由程序直接生成const static int primes = 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,.65521, 65537;有点不可思议吧, 原本makePrimes需要花费的时间复杂度现在真的变成O(1)了!(我觉得叫O(0)可能更合适!)7. 二分法查找现在我们缓存了前大约sqrt(232)/(ln(sqrt(232)-3/2)个素数列表, 在判断232级别的素数时最多也只需要PI(sqrt(232)次判断(准确值是6543次), 但是否还有其他的方式判断呢?当素数比较小的时候(不大于216), 是否可以直接从缓存的素数列表中直接查询得到呢?答案是肯定的! 由于primes是一个有序的数列, 因此我们当素数小于216时, 我们可以直接采用二分法从primes中查询得到(如果查询失败则不是素数).代码:/ 缺少的代码请参考前边#include static bool cmp(const int *p, const int *q) return (*p) - (*q);bool isPrime(int n) if(n = 67 & n = primesNELEMS(primes)-1) return NULL != bsearch(&n, primes, NELEMS(primes), sizeof(n), cmp); else for(int i = 1; primesi*primesi = n; +i) if(n%primesi = 0) return false; return true; 时间复杂度:if(n = 67): O(log2(NELEMS(primes) primesNELEMS(primes)-1): O(PI(sqrt(n) = NELEMS(primes).8. 素数定理+2分法查找在7中, 我们对小等于primesNELEMS(primes)-1的数采用2分法查找进行判断.我们之前针对232缓冲的6453个素数需要判断的次数为13次(log2(1024*8) = 13). 对于小的素数而言(其实就是216范围只内的数), 13次的比较已经完全可以接受了.不过根据素数定理: ln(x)-3/2 x/PI(x) = 67时, 我们依然可以进一步缩小小于232情况的查找范围(现在是0到NELEMS(primes)-1范围查找).我们需要解决问题是(n = primesNELEMS(primes)-1):如果n为素数, 那么它在素数序列可能出现的范围在哪?- (n/(ln(n)-1/2), n/(ln(n)-3/2), 即素数定理!上面的代码修改如下:代码:bool isPrime(int n) if(n = 67 & hi NELEMS(primes) int lo = (int)floor(n/(ln(n)-1/2); return NULL != bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp); else for(int i = 1; primesi*primesi = n; +i) if(n%primesi = 0) return false; return true; 时间复杂度:if(n = 67): O(log2(hi-lo) primesNELEMS(primes)-1): O(PI(sqrt(n) = NELEMS(primes).9. 打包成素数库(给出全部的代码)到目前为止, 我已经给出了我所知道所有改进的方法(如果有人有更好的算法感谢告诉我).这里需要强调的一点是, 这里讨论的素数求法是针对0-232范围的数而言, 至于像寻找成百上千位大小的数不在此讨论范围, 那应该算是纯数学的内容了.代码保存在2个文件: prime.h, prime.cpp.代码:/ file: prime.h#ifndef PRIME_H_2006_10_27_#define PRIME_H_2006_10_27_extern int Prime_max(void); / 素数序列的大小extern int Prime_get (int i); / 返回第i个素数, 0 = i Prime_maxextern bool Prime_test(int n); / 测试是否是素数, 1 = n INT_MAX#endif/ file: prime.cpp#include #include #include #include #include prime.h/ 计算数组的元素个数#define NELEMS(x) (sizeof(x) / (sizeof(x)0)/ 素数序列, 至少保存前6543个素数!static const int primes = 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,.65521, 65537;/ bsearch的比较函数static int cmp(const void *p, const void *q) return (*(int*)p) - (*(int*)q);/ 缓冲的素数个数int Prime_max() return NELEMS(primes);/ 返回第i个素数int Prime_get(int i) assert(i = 0 & i 0); if(n 2) return false; if(n = 2) return true; if(!(n&1) return false; / 如果n为素数, 则在序列hi位置之前 int lo, hi = (int)ceil(n/(log(n)-3/2.0); if(hi = 67是才满足素数定理 if(n = 67) lo = (int)floor(n/(log(n)-1/2.0); else lo =