取自然数中的素数几种算法

ACM中的经典问题。网上看了大神的的各种算法，整理一下。
主要参考：http://blog.csdn.net/once_hnu/article/details/6302283
http://blog.csdn.net/liukehua123/article/details/5482854（这篇博文里有几处错误，可以看下面的评论）
http://blog.csdn.net/morewindows/article/details/7354571（提供利用位运算压缩素数表）
感谢大神提供的思路。

一下代码均为java语言，且为了简洁，只贴出了核心部分，只要把核心部分放在一个类中，用main函数调用即可。

第一种，根据素数的定义，就是直接除以前面的数，看看是否能整除。

/*
n为自然数的上界，下界默认为1。

array[]为存储素数的数组，方便输出。当数量级大的时候，均不输出。

知识：如果一个数是合数，那么一定至少有一个除了1以外的因子小于等于他的平方根。这点很重要，因为只要判断出一个因子，就可以确定他是合数了。这个结论大大减少了循环次数。

以上结论在后面的代码均适用，就不在重复。
*/
	void tradition(int n,int array[]){
		for (int i=2;i<n;i++){//i为需要判断是否为素数的自然数
			int j=2;//j为除数
			for (;j<=Math.sqrt(i);j++){
				if(i%j==0){//取余，等于0就说明j是i的一个因子
					break;
				}
			}
			int mun=0;//用以循环赋值给array数组时的递增量
			//j大于sqrt(i)说明j是自然循环出来的，而不是break出来的，就说明i没有小于等于sqrt(i)的因子，所以i是素数
			if(j>Math.sqrt(i)){
				array[mun]=i;
				mun++;
			}
		}
	}

传统算法容易理解，但是效率非常低，当n为千万时，耗时已经达到20s左右。

第二种，在第一种的基础上，把偶数给排除，效率快了很多。

	void remove_even(int n,int array[]){
		array[0]=2;//2是偶数，但也是素数，所以特殊对待。
		//和第一种区别就在于i+=2这里，也就是跳过了偶数
		for (int i=3;i<n;i+=2){
			int j=3;
			//只需要除以奇数，因为偶数乘以任何自然数数都是偶数，不可能得到奇数。
			for (;j<=Math.sqrt(i);j+=2){
				if(i%j==0){
					break;
				}
			}
			int mun=1;
			if(j>Math.sqrt(i)){
				array[mun]=i;
				mun++;
			}
		}
	}

较第一种优化了一点点，千万时用时12s左右。

第三种，这里的算法就牛逼了，体现差距的地方就在这里了。

	void shaixuan(int n,int array[]){
		//定义一个布尔类型数组，他的下标就代表自然数
		boolean prime[]=new boolean[n];//默认值为false
		prime[2]=true;//2特殊对待
		//这一步是为了把奇数标为true
		for (int i=3;i<n;i++){
			if ((i&1)==1)//奇数的位运算判断方法
				prime[i]=true;//奇数就标为true
		}
		//这一步的具体原理可以去看前面给出的两篇博文。
		for (int i=3;i<=Math.sqrt(n);i++){//i就是数组的下标，也就是要判断的自然数
			if (prime[i]){//true就说明是奇数，需要判断是否为素数
				for(int j=i*i;j<n;j+=2*i)//后面具体解释
					prime[j]=false;//把i的倍数都标为false
			}
		}
		//同上，为了赋值给array	
		int mun=0;
		for(int i=2;i<n;i++){
			if (prime[i]){
				array[mun]=i;
				mun++;
			}
		}	
	}

这种方法看起来更长了，代码也多了，但是效率比上面的快了不是一点两点。千万的速度是150ms左右。自己体会下效率快了多少。
引用上面博文的解释：
一个简单的筛素数的过程：n=30。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
第 1 步过后4 6 … 28 30这15个单元被标成false,其余为true。
第 2 步开始：
i=3; 由于prime[3]=true, 把prime[9], [15], [21], [27], [30]标为false.
i=4; 由于prime[4]=false,不在继续筛法步骤。
i=5; 由于prime[5]=true, 把prime[15],[25]标为false.
i=6>sqrt(30)算法结束。
第 3 步把prime[]值为true的下标输出来：
for(i=2; i<=30; i++)
if(prime[i]) printf("%d ",i);
结果是 2 3 5 7 11 13 17 19 23 29

自己理解一下，以3为例，就是把因子里包括3的合数都给标记为false了。

for(int j=i*i;j<n;j+=2*i)

这一句的判定是这样的：还是以3为例，当i=3时，那以3为因子的合数就是6，9，12，15…，也就是3*2，3*3，3*4，3*5…，也就是3+3，3+2*3，3+3*3，3+4*3…这里再省一步，因为6，12都是偶数，在之前已经标记为false了，所以可以直接跳过，就变成3+2*3，3+4*3…

int j=i*i，以5为例：5*2，5*3，5*4，5*5，5*6…，乘以偶数的就不用考虑了。乘奇数的，小于5的奇数，那肯定也被这个奇数或者这个奇数的因子去掉了。这样就避免重复了。所以起步就从i*i开始了。

j+=2*i就是上面解释的。

第四种，在第三种的基础上再优化一点。就是前面给出的博文里提到的博士独创的方法。

/*
这里的思路是在预设布尔数组的时候，就已经考虑把偶数的排除掉了，
也就是认为里面只有奇数。存储的数据是这样的：
2，3，5，7，9，11，13，15，17，19，21....
且与数组下标对应：
0，1，2，3，4， 5，6， 7， 8， 9， 10....
也就是下意识要认为其中的prime[0]=2;prime[1]=3....
这样的做法是省去了到偶数哪一步的判断。
*/
	void shaixuan(int n, int array[]){
		//因为只有奇数，所以对n的自然数只有一半的奇数，当然这里很可能会丢掉最尾巴的那个数，这里就不做那么细的考虑了。
		boolean prime[]=new boolean[n>>1];//除以2的位运算方式
		
		for(int i=1;i<=Math.sqrt(n>>1);i++){
			//这一步是因为布尔数组默认值为false，为了简便，这里就认为false是要的，true是不要的。
			if(prime[i]==false){
				for(int j=2*i*i+2*i;j<n>>1;j+=2*i+1){//这一步判断比较绕，后面解释。
					prime[j]=true;//思路和上面一样，把倍数赋值为true
				}
			}
		}
		
		int mun=1;
		array[0]=2;
		for(int i=1;i<n>>1;i++){
			if (prime[i]==false){//false才是素数
				array[mun]=2*i+1;//i是下标，2*i+1才是实际的数值
				mun++;
			}
		}
	}

千万的速度是110ms左右，是不是又快了一点。

其实思路还是第三种的思路，只是这里对偶数的处理更加巧妙了。
for(int j=2*i*i+2*i;j>1;j+=2*i+1)其实判断的本质和原来的判断是一样的，只是这时候不能认为下标就是数字，还需要把下标转化成实际的数值。以3为例，3和3的倍数：
数组的下标：1，4， 7， 10…
实际的数值：3，9，15，21…
思路还是原来的，想把3的倍数赋值为true（这里true和false是反过来的）。那就是把9 ，15，21…赋值为true，但是并不是直接把prime[9],prime[15],prime[21]…赋值为true,还需要找到这些数值对应的下标，也就是prime[4],prime[7],prime[10]…赋值为true。

假设我不知道9对应的下标是多少，设为x。原来的int j =i*i,那这里就是：
(2*i+1)*(2*i+1),但是这只是值，不是下标，还要解(2*i+1)*(2*i+1)=（2*x+1）,这里把x解出来，才是下标，x=2*i*i+2*i。把i=1带进去验证一下，x=4，i=1对应的值是3，x=4对应的值就是9，没错吧。所以int j=2*i*i+2*i就是j的初始值。

对于j+=2*i+1;这个看规律也能看出来。非要计算的话,设要求的式子为j+=x;
原来的是：j+=2*i。对应这里的i则要变化为2*（2*i+1）。意思就是两个紧邻的倍数之间的差值为2*（2*i+1）。
其实对应的式子是：prime[j+x]-prime[j]=2*prime[i]
对应的数值是2*（j+x）+1-(2*j+1)=2*(2*i+1),解得x=2*i+1。

第五种，利用位操作压缩素数表，并且效率也提高了。

	void compress_save(int n,int array[]){
		for(int i=1;i<=Math.sqrt(n>>1);i++){
			if((array[i/32]>>(i%32)&1)==0){
				for(int j=2*i*i+2*i; j<n>>1; j+=2*i+1){
					array[j/32]|=(1<<j%32);
				}
			}
		}
	}

千万的速度是55ms左右，这个效率已经够苛刻了吧。

解释一下原理：
其实实际原理和第四种一模一样，只是存储数据的方式变了，也就是所谓的数据结构吧。第四种采用布尔数组存储奇数，对小于n的自然数，那就是需要n/2长度的布尔数组，布尔数组当做byte数组处理（网上资料，不敢保证一定正确）。也就是n/2个字节。int数据在内存中为32位二进制数，如果我把每一位数都对应为一个自然数，就是一个int数据就可以存储32个整数（实际有一点点变化，后面讲），因为数组在内存中是连续的，所以可以把数组看成是一个整数。例如数组长度为2，那么就是连续的两个int数据在一起，如果第一int数据的第一位标记为0，那最后一位为31，紧接着的下一个int数据的第一位是不是就可以看成是第32位，最后一位看成是63位？这样就是存储了64位数了。一次类推，数组长度变成n，可以存储的数就变成了32*n。这样就大大节省了空间，特别是素数问题里都是要求超大型数据。

我认为的数据存储依然是这样的：
2，3，5，7，9，11，13，15，17，19，21…
且与数组下标对应：
0，1，2，3，4， 5， 6， 7， 8， 9， 10…
实际操作其实就是把数组里的一个数看成32个数，也就是array[0]的第一位为2，第二位为3，仍然满足2*i+1的关系（除了第一位）。然后操作就是之前介绍的。

解释下那几句位运算：
①(array[i/32]>>(i%32)&1)==0
i是从0开始的，当i=31时，对应的刚好是aarray[0]的最后一位，i=32时，对应的是array[1]的第一位，i/32=1,已经到array[1]了，所以这里不是除以31而是32。
这句话的作用是：例如i=5,先将array[5/32]=array[0]这个数右移5%32=5位，再与1进行与运算，判断这一位是否为0（与1进行与运算是位运算里用的最多的一种方式之一）。因为默认的所有位都为0，判断为合数的为就会置1，这一步判断为0的数就是素数。就相当于第四种中的prime[i]==false，功能一模一样看懂第四种，这里就好理解了。

②int j=2*i*i+2*i; j<n>>1; j+=2*i+1和第四种一样，没变。

③array[j/32] |= (1<<j%32);功能与第四种的prime[j]=true;一模一样，只不过这里对某一位赋值，比较特殊而已。具体原理可以看给的第三个链接。

还得解释一下如何输出素数，这里不像前面，都把它存在一个数组里了，直接输出数组就完事了。需要稍微转一下弯，方式如下：

	System.out.print(2+" ");//2总是特殊对待
	for(int i=1; i<n>>1;i++)
	//判断第i位是不是0，是0说明这一位存储的是素数，而i对应的数值是2*i+1,所以输出2*i+1即可
		if(((array[i/32]>>i%32)&1)==0)
			System.out.print(2*i+1+" "); 

下面是我所有代码采用的统一的main函数写法，修改n的值就可以修改数量级，当然定义的类和函数的名称不一样，用到时记得修改。

public class main {
	public static void main(String[] args) {
		filter s1=new filter();
		int n=10000000;
		int array[] = new int[n>>1];
		long startTime = System.currentTimeMillis();//获取当前时间
		s1.shaixuan(n, array);
//		for(int i=0; i<array.length;i++)
//			System.out.print(array[i]+" ");
		long endTime = System.currentTimeMillis();//获取当前时间
		System.out.println("程序运行时间："+(endTime-startTime)+"ms");
	}
}

最后再扯一下array数组大小的问题。在main函数里，定义是这样的：

int array[] = new int[n];

五种方法里，大小的取值依次为：
①n，②n，③n/2，④n/2，⑤n/40+1
主要分析第5种，理论上n/64是刚好的（n/2再/32位），但是这里不能刚好。

		for(int i=1;i<=Math.sqrt(n>>1);i++){
			if((array[i/32]>>(i%32)&1)==0){
				for(int j=2*i*i+2*i; j<n>>1; j+=2*i+1){
					array[j/32]|=(1<<j%32);

举例：n=100时，sqrt(n>>1)=7,n/64=1。上面代码，当i=6时，(array[i/32]>>(i%32)&1)==0，此时j=2*6*6+2*6=84, array[84/32]=array[2], 然而array长度只为1，所以会报错。
sqrt(n>>1)是不能优化了（不会），那就只能适当的放大array的长度，n/40+1是我随便测试出来的，并不是完全吻合，懒得去追求那么完美了。+1主要是照顾n<40的时候。

ps:以上几种算法均不是自创，借鉴前面提到的三篇博文的思路，自己一点点写的。难免有错误的地方，还请不吝指教。如果有其他更好的方法，更请赐教。

（为什么有部分代码的注释不会变成红色了？不懂什么情况）