筛选法求质数

#include<stdio.h>
#include<time.h>
#include<stdlib.h>
#define m 10000001
int main()
{
	clock_t start, finish;
	int i,j,p[m],count;
	start=clock();
	for(i=0;i<m;i++)
		p[i]=1;
	for(i=2;i*i<m;i++)
		if(p[i])
		{
			for(j=i+i;j<m;j+=i)
				p[j]=0;
		}
	count=0;
	for(i=2;i<m;i++)
		if(p[i])
			count++;
	finish=clock();
	printf("%f.3",(double)(finish - start) / CLOCKS_PER_SEC);
	return 0;
}

质数算法大全，及C程序实现优化详解  筛选法

三木追风

质数的定义

一个数，如果只有1和它本身两个因数，这样的数叫做质数，又称素数。

在上文 《素数算法大全，及C程序实现优化详解 (一) 试除法》中我们已经探讨了求解素数的一类算法，并且将试除法从最初的低效版本优化的高效的V2。那么，还有没有其它更佳算法呢？这就是下面三藏要和大家探讨的内容

合数过滤筛选法

算法描述：我们知道，素数N不能被2~(N-1)间的任何数整除；反过来看，只要能被2~(N-1)间的任何数整除的N，都不是素数。所以我们可以采用一个简单的排除法：就是对N以内的所有数，只要逐个去除值为2~(N-1)的倍数的数，剩下的就是素数。

C语言实现

// 合数过滤筛选法 Ver1 
// 参数：n 求解n以内(包括n)的素数
// 返回值：n以内素数个数 
int CompositeNumFilterV1(int n)
{
 int i, j;
 // 素数数量统计 
 int count = 0;
 // 分配素数标记空间，结合后文思考为何+1
 char* flag = (char*)malloc( n+1 ); 
 
 // 初始化素数标记 
 for (i=2; i<=n; i++)
 {
  // 为什么*(p+i)要写成flag[i]呢？可读性更佳尔 
  flag[i] = 1;
 }
 
 // 写程序要注意排版和留空，方便阅读，也可减少出错几率
 // 以2~(N-1)为因子过滤合数 
 for (i=2; i < n; i++)
 {
  for (j=2; i*j <= n; j++)
  {
   // i*j是由i,j两整数相乘而得，显然不是素数
   flag[i*j] = 0;
  }
 }
 
 // 统计素数个数
 for (i=2; i<=n; i++)
 {
  // 其实if(flag)就其同样作用了，但这么写是有留言的 
  // 请参阅《C语言程序设计常见错误剖析及解决之道》一文 
  if (1 == flag[i]) count++;
 }
  
 // 因输出费时，且和算法核心相关不大，故略
  
 // 释放内存，别忘了传说中的内存泄漏 
 free(flag);
 
 return count;
} 

在上文给出的main函数中以不同参数调用CompositeNumFilterV1函数，得到执行结果如下：

[100000]以内素数个数：9592, 计算用时：15毫秒
[1000000]以内素数个数：78498, 计算用时：125毫秒
[5000000]以内素数个数：348513, 计算用时：2578毫秒
[10000000]以内素数个数：664579, 计算用时：6281毫秒

注：因程序是非独占性运行的，所以时间不是完全精确的，但基本能反映实情

显然，比上文中的试除法要快，而且谁都可以看到上例是一个未经优化的粗陋版本，好多地方是三藏故意采用比较低效做法，为了与后文的优化版比较，凸显优化之重要，也为了初学者记住别采用类似低效做法，下面我们开始优化之旅

优化分析

上面CompositeNumFilterV1函数存在的问题有：

1. 在外层循环，需要一直执行到n-1吗？不要，因为n/2~n-1间的数显然不能整出n 

2. 在内层循环中重复使用i*j显然是低效的，考虑到计算机中加减运算速度比乘除快，可以考虑变乘法为加法 

3. 在循环修改flag过程中，其实有很多数会被重复计算若干次，比如6=2*3=3*2，会被重复置0，类似操作很多，所以我们得设法避免或减少flag重复置0

据上述分析，我们可将程序优化如下：

// 合数过滤筛选法 Ver2 
// 参数：n 求解n以内(包括n)的素数
// 返回值：n以内素数个数 
int CompositeNumFilterV2(int n)
{
 int i, j;
 // 素数数量统计 
 int count = 0;
 // 分配素数标记空间，明白+1原因了吧，因为浪费了一个flag[0]
 char* flag = (char*)malloc( n+1 ); 
 
 // 初始化素数标记，要高效点咯
 flag[2] = 1;
 // 注意是i<n不是上例中的i<=n了，理由自思 
 for (i=3; i<n; i++)
 {
  flag[i++] = 1;
  // 偶数自然不是素数，直接置0好了 
  flag[i] = 0;
 }
 // n为奇数 
 if (n%2 != 0)
 {
  flag[n] = 1;
 }
 
 // 从3开始filter，因为2的倍数早在初始化时代就干掉了
 // 到n/2止的理由还要说吗 
 for (i=3; i <= n/2; i++)
 {
  // i是合数，请歇着吧，因为您的工作早有您的质因子代劳了 
  if (0 == flag[i]) continue;
  
  // 从i的2倍开始过滤，变乘法为加法  
  for (j=i+i; j <= n; j+=i)
  {
   flag[j] = 0;
  }
 }
 
 // 统计素数个数
 for (i=2; i<=n; i++)
 {
  if (flag[i]) count++;
 }
  
 // 因输出费时，且和算法核心相关不大，故略
  
 // 释放内存，别忘了传说中的内存泄漏 
 free(flag);
 
 return count;
} 

再来调用CompositeNumFilterV2得到执行结果：

[100000]以内素数个数：9592, 计算用时：n太小，时间精度不够
[1000000]以内素数个数：78498, 计算用时：31毫秒
[5000000]以内素数个数：348513, 计算用时：453毫秒
[10000000]以内素数个数：664579, 计算用时：1062毫秒
[100000000]以内素数个数：5761455, 计算用时：12973毫秒

哇哇，比昨天的试除发快了好多倍，可见算法的威力，值得好好学习，别说学算法没用咯。

上例着那个计算一亿以内的素数只要约13秒，应该算不错了，今天是否可以休息了呢？No，我们要追求极限！

int CompositeNumFilterV3(int n)
{
 int i, j;
 // 素数数量统计 
 int count = 0;
 // 分配素数标记空间，明白+1原因了吧，因为浪费了一个flag[0]
 char* flag = (char*)malloc( n+1 );
 // 干嘛用的，请仔细研究下文
 int mpLen = 2*3*5*7*11*13;
 char magicPattern[mpLen];
 // 奇怪的代码，why，思考无法代劳，想！ 
 for (i=0; i<mpLen; i++)
 {
  magicPattern[i++] = 1;
  magicPattern[i++] = 0;
  magicPattern[i++] = 0;
  magicPattern[i++] = 0;
  magicPattern[i++] = 1;
  magicPattern[i] = 0;
 }
 for (i=4; i<=mpLen; i+=5)
  magicPattern[i] = 0;
 for (i=6; i<=mpLen; i+=7)
  magicPattern[i] = 0;
 for (i=10; i<=mpLen; i+=11)
  magicPattern[i] = 0;
 for (i=12; i<=mpLen; i+=13)
  magicPattern[i] = 0;
 
 // 新的初始化方法,将2,3,5,7,11,13的倍数全干掉
 // 而且采用memcpy以mpLen长的magicPattern来批量处理 
 int remainder = n%mpLen;
 char* p = flag+1;
 char* pstop = p+n-remainder;
 while (p < pstop)
 {
  memcpy(p, magicPattern, mpLen);
  p += mpLen;
 }
 if (remainder > 0)
 {
  memcpy(p, magicPattern, remainder);
 }
 flag[2] = 1;
 flag[3] = 1;
 flag[5] = 1;
 flag[7] = 1;
 flag[11] = 1;
 flag[13] = 1;
 
 // 从17开始filter，因为2,3,5,7,11,13的倍数早被kill了 
 // 到n/13止的，哈哈，少了好多吧
 int stop = n/13;
 for (i=17; i <= stop; i++)
 {
  // i是合数，请歇着吧，因为您的工作早有您的质因子代劳了 
  if (0 == flag[i]) continue;
  
  // 从i的17倍开始过滤
  int step = i*2;
  for (j=i*17; j <= n; j+=step)
  {
   flag[j] = 0;
  }
 }
 
 // 统计素数个数
 for (i=2; i<=n; i++)
 {
  if (flag[i]) count++;
 }
  
 // 因输出费时，且和算法核心相关不大，故略
  
 // 释放内存，别忘了传说中的内存泄漏 
 free(flag);
 
 return count;
}

再看CompositeNumFilterV3执行结果：

[1000000]以内素数个数：78498, 计算用时：15毫秒
[5000000]以内素数个数：348513, 计算用时：203毫秒
[10000000]以内素数个数：664579, 计算用时：515毫秒
[100000000]以内素数个数：5761455, 计算用时：6421毫秒

再次优化后速度提升了又一倍左右，三藏不禁有点满足了，睡觉也！

素数筛选法优化

邓达生

原始的素数筛选法是在{2,……,n}里进行筛选。代码：

(注：本文章的代码都不进行for循环优化，例如上述代码改为

若要优化请自行优化)

本文章的优化仅仅是加快筛选的速度(因为基础都是筛选)，如果原算法时间和空间复杂度是f(n)和g(n)，则优化后为k1f(n)和k2g(n)，k1、k2＜1。

优化：

因为2是已知的素数，所以{2,……,N}里的所有偶数(2的倍数)都可以在筛选前排除，也就是说待筛选列表里不应该存在2的倍数，并且在筛选中不应该存在检测2的倍数是否素数的行为或者将2的倍数标记为合数的行为，例如当检测到3是素数后，要把3的倍数都标记为合数时，不应该把6标记为合数，因为在待筛选列表里没有6。

我的基本思想就是在原始的待筛选列表里删除某(几)个已知素数的倍数。设删除后的待筛选列表为D，很容易证明若x、kx∈D，则k∈D。证明：若k∉D，则k是某个已知素数的倍数，则kx是某个已知素数的倍数，所以kx∉D，与已知kx∈D矛盾，所以k∈D。

现在原始待筛选列表删除2的倍数(即D={3,5,7,9……})，对于标记数组mark[]，mark[i]代表的就已经不能是i了，mark[i]代表的是2i+1是否合数。

筛选过程：假设现在mark[i]是false，代表2i+1是素数。然后就要把2i+1的倍数都标记为合数。应该被标记为合数的数分别是

应标记的数/2i+1的倍数	对应mark的下标
3(2i+1)=2(3i+1)+1	3i+1
5(2i+1)=2(5i+2)+1	5i+2
7(2i+1)=2(7i+3)+1	7i+3
……	……
(2k+1)(2i+1)=2((2k+1)i+k)+1	(2k+1)i+k

所以待标记的数的下标初始值为3i+1，增量为2i+1.

算法代码如下(m = n / 2)：

对于求n以内的素数，原始算法m=n，标记合数for循环增量是i，优化后的算法m=n/2，标记合数for循环增量是2i+1，效率提高大约一倍。

继续优化：原始待筛选列表删除2的倍数后m=n/2，再删除3的倍数后m=n*(2/6)=n/3，删除5的倍数后m=n*(8/30)=4n/15，……可以看出来，到后面即使再优化，效率也不会提高多少，而其后可以看到，效率在以很小的幅度提高时，代码的复杂度却以很高的速度提高。(以上为个人简单估计)

再删除3的倍数：D={5,7,11,13,17,19……}

令mark[i]代表6i±1是否合数。char类型是1字节8位，用低两位来标记(0位代表6i-1,1位代表6i+1)。

6i-1的倍数	mark下标	mark下标增量	位
6i-1	i		0
5(6i-1)=30i-5=6(5i-1)+1	5i-1	4i-1	1
7(6i-1)=42i-7=6(7i-1)-1	7i-1	2i	0
11(6i-1)=66i-11=6(11i-2)+1	11i-2	4i-1	1
13(6i-1)=78i-13=6(13i-2)-1	13i-2	2i	0
……	……	……	……
(6j-1)(6i-1)=6((6j-1)i-j)+1	(6j-1)i-j	4i-1	1
(6j+1)(6i-1)=6((6j+1)i-j)-1	(6j+1)i-j	2i	0

6i+1的倍数	mark下标	mark下标增量	位
6i+1	i		1
5(6i+1)=30i+5=6(5i+1)-1	5i+1	4i+1	0
7(6i+1)=42i+7=6(7i+1)+1	7i+1	2i	1
11(6i+1)=66i+11=6(11i+2)-1	11i+2	4i+1	0
13(6i+1)=78i+13=6(13i+2)+1	13i+2	2i	1
……	……	……	……
(6j-1)(6i+1)=6((6j-1)i+j)-1	(6j-1)i+j	4i+1	0
(6j+1)(6i+1)=6((6j+1)i+j)+1	(6j+1)i+j	2i	1

最终算法代码(m = n / 6)：

再删除5的倍数：D={7,11,13,17,19,23,29,31……}，而mark[i]代表30i+1、30i+7、30i+11、30i+13、30i+17、30i+19、30i+23、30i+29 共8个数，刚刚一个字节可以存完，没有浪费任何空间。

注：除法 "/" 是整除，"%"运算是模运算(即求余)

30i+k的30j+c倍数	mark下标增量
(30j+1)(30i+k)=30((30j+1)i+jk)+k
(30j+7)(30i+k)=30((30j+7)i+jk+7k/30)+7k%30	6i+7k/30
(30j+11)(30i+k)=30((30j+11)i+jk+11k/30)+11k%30	4i+11k/30-7k/30
(30j+13)(30i+k)=30((30j+13)i+jk+13k/30)+13k%30	2i+13k/30-11k/30
(30j+17)(30i+k)=30((30j+17)i+jk+17k/30)+17k%30	4i+17k/30-13k/30
(30j+19)(30i+k)=30((30j+19)i+jk+19k/30)+19k%30	2i+19k/30-17k/30
(30j+23)(30i+k)=30((30j+23)i+jk+23k/30)+23k%30	4i+23k/30-19k/30
(30j+29)(30i+k)=30((30j+29)i+jk+29k/30)+29k%30	6i+29k/30-23k/30
(30(j+1)+1)(30i+k)=30((30j+31)i+jk+k)+k	2i+k-29k/30

这样的表一共有8个，共有8*8=64项，所以我之前说过代码复杂度会提高很多。如果就想前面那样写代码的话这个函数就会很长，所以就制表来解决。

代码很长，后面再列出来。为了测试上述3种算法的效率，我在中山大学的在线测评系统sicily的第一题"A-B"测试了运行效率，各自计算33554432(225)内的素数(共2063689个)，只删除2的倍数的用时大约0.9s，再删除3的倍数的用时大约0.6s，再删除5的倍数的用时只有0.2s！（全部代码在本机测试结果没错后再去sicily测试时间）难道是n不够大？还是我的估计误差太大呢？不过由于代码复杂度都已经比较高了，而且内存空间得到充分利用(mark的每一项代表8个数，刚好一个字节8位二进制用完)，我就不深入探讨了，有兴趣的读者可以自己算一下。下面附上最后一份代码(m = n / 30)：