求素数

今天，关于素数问题纠结了好久好久，倍感知识缺乏啊。因此，通过自己的了解和网上查阅资料，加上自己的啰嗦，在这里整理一下，日后可以翻阅。

   首先，感谢网上的前辈，如果没有您们，我不会获得关于素数的比较全面的知识。非常感谢。

   

1、素数及相关

   素数，又称质数，在一个大于1的自然数中，除了1和此整数自身之外，不能被其他自然数整除的数。

   比1大但不是素数的数称为合数。

   1和0既不是素数，也不是合数。

   算术基本定理证明每个大于1的正整数都可以写成素数的乘积，并且这种乘积的形式是唯一的。

2、试除法求素数

   算法描述：根据素数的定义可知，不能被1和自身外的整数整除的数为素数。所以，我们可以得知，判断一个素数是否为素数只要看它是否能被2~sqrt(i)间的数整数即可。而求N内所有素数则是循环重复上述过程。

  

   C语言实现如下所示。

• #include <stdio.h>
• #include <time.h>
• #include <math.h>
• #include <malloc.h>
• //试除法
• #define NUM 10000
• int Test_prime(int n)
• {
•     int count= 0;
•     int i, j;
•     int *num =(int*)malloc(sizeof(int)* n);
•     num[count++]= 2;
•      
•     for(i= 3; i<= n; i++)
•     {
•         for(j= 2; j<= sqrt(i); j++)
•         {
•             if(i% j == 0)
•             {
•                 break;
•             }
•         }
•         if(j> sqrt(i))
•         {
•             num[count++]= i;
•         }
•     }
•     free(num);
•     return count;
• }
• int main()
• {
•     int count;
•     clock_t start,end;
•     start = clock();
•     count = Test_prime(NUM);
•     end = clock();
•     printf("%d内的素数个数为：%d, 总共耗时为：%d毫秒\n", NUM, count,end - start);
•     return 0;
• }

   测试结果如下所示（测试在VC6.0下进行）。

   从上面可以看出，当数据很大时，时间消耗增长的比较快。

3、试除法的优化方案

   

   仔细研究试除法，可以发现以下几个问题：

   1> 在循环条件中重复调用了sqrt(i)，这显然是比较浪费时间的；

   2> 判断素数，真的需要拿2-sqrt(i)间的所有整数去除吗？我们知道，合数都可以分解成若干质数，所以，只要2-sqrt(i)间的质数不能整除i即可；

   C语言实现如下所示。

点击(此处)折叠或打开

• //求N内的所有素数
• #include <stdio.h>
• #include <time.h>
• #include <math.h>
• #include <malloc.h>
• //试除法
• #define NUM 1000000
• int Test_prime(int n)
• {
•     int count= 0;
•     int i, j, k, stop;
•     //分配空间
•     int *num =(int*)malloc(sizeof(int)* n);
•     //2肯定是素数
•     num[count++]= 2;
•     stop = 0;
•     for(i= 3; i<= n; i++)
•     {
•         //在循环中重复调用sqrt是低效做法，故引入k
•         k =(int)sqrt(i);
•         //stop的作用是：统计小于当前k值的质数的数目
•         while(num[stop]<= k && stop< count)
•         {
•             stop++;
•         }
•         for(j= 0; j< stop; j++)
•         {
•             if( i% num[j]== 0)
•             {
•                 //i不能被2-sqrt(i)间的素数整除，自然不能被其他整数整除，所以为素数
•                 break;
•             }
•         }
•         if(j== stop)
•         {
•             num[count++]= i;
•         }
•         
•     }
•     free(num);
•     return count;
• }
• int main()
• {
•     int count;
•     clock_t start,end;
•     start = clock();
•     count = Test_prime(NUM);
•     end = clock();
•     printf("%d内的素数个数为：%d, 总共耗时为：%d毫秒\n", NUM, count,end - start);
•     
•     return 0;
• }

   测试结果如下所示。

   相对于优化前的算法，时间提供了很多。特别是在时间增长曲线的幅度变小了，N值越大，优化后的算法比优化后的算法效率更高。

4、合数过滤筛选法

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

   C语言实现如下所示。

• //合并筛选法
• #include <stdio.h>
• #include <time.h>
• #include <math.h>
• #include <malloc.h>
• //试除法
• #define NUM 10000
• int Test_prime(int n)
• {
•     int count= 0;
•     int i, j;
•     //分配空间，之所以是n+1，是因为浪费了一个num[0]
•     char *num= (char*)malloc(sizeof(char)* (n+ 1));
•     //初始化素数标记
•     for(i= 2; i<= n; i++)
•     {
•         num[i]= 1;
•     }
•     //以2-(N-1)为因子过滤合数
•     for(i= 2; i<= n-1; i++)
•     {
•         for(j= 2; j* i <= n; j++)
•         {
•             //i*j是由两整数相乘而得，显然不是素数
•             num[i*j]= 0;
•         }
•     }
•     //统计素数个数
•     for( i= 2; i<= n; i++)
•     {
•         if( 1== num[i])
•         {
•             count++;
•         }
•     }
•     free(num);
•     return count;
• }
• int main()
• {
•     int count;
•     clock_t start,end;
•     start = clock();
•     count = Test_prime(NUM);
•     end = clock();
•     printf("%d内的素数个数为：%d, 总共耗时为：%d毫秒\n", NUM, count,end - start);
•     
•     return 0;
• }
•    测试结果如下所示。

   上述程序好多地方采用了比较低效的做法，为了与后文的优化作比较，这也是像我一样的初学者通常采用的版本，因此，要学会优化。

5、合并筛选法优化方案

   上述算法存在的问题是：

   1> 在外层循环，需要一直执行到n-1嘛？不要，因为n/2-（n-1）之间的数显然不能整除出n；

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

   3> 在循环修改flag的过程中，其实有很多数被重复计算若干次，比如6 = 2*3 = 3*2,被重复置零，所以，可以进行避免；

   

   C语言实现如下所示。

点击(此处)折叠或打开

• //合并筛选法的优化方案
• #include <stdio.h>
• #include <time.h>
• #include <math.h>
• #include <malloc.h>
• #define NUM 300000
• int Test_prime(int n)
• {
•     int count= 0;
•     int i, j;
•     //分配空间
•     char *num= (char*)malloc(sizeof(char)* (n+ 1));
•     
•     //初始化素数标记
•     num[2]= 1;
•     //注意此处是i<n，上例中的i<=n
•     for(i= 3; i< n; i++)
•     {
•         num[i++]= 1;
•         num[i]= 0;//偶数自然不是素数
•     }
•     //如果n为奇数
•     if(n% 2 != 0)
•     {
•         num[n]= 1;
•     }
•     //从3开始过滤，因为，2的倍数在初始化中去掉了
•     for(i = 3; i <= n/2; i++)
•     {
•         if(0 == num[i] )
•         {
•             continue;
•         }
•         //从i的2倍开始过滤
•         for(j = i + i; j <= n;j+=i)
•         {
•             num[j] = 0;
•         }
•     }
•     //统计素数个数
•     for( i= 2; i<= n; i++)
•     {
•         if( 1== num[i])
•         {
•             count++;
•         }
•     }
•     free(num);
•     return count;
• }
• int main()
• {
•     int count;
•     clock_t start,end;
•     start = clock();
•     count = Test_prime(NUM);
•     end = clock();
•     printf("%d内的素数个数为：%d, 总共耗时为：%d毫秒\n", NUM, count,end - start);
•     
•     return 0;
• }

   测试如下所示。

   确实比先前快了很多，优化真的可以带来时间的提高，这样我很是欣喜。

   

   后来想到进行添加补充：

  

   如果我对上述红色部分代码进行优化，如下所示。

点击(此处)折叠或打开

• //从3开始过滤，因为，2的倍数在初始化中去掉了
•     for(i= 3; i<= n/2;i = i + 2)
•     {
•         //在这里进行判断，就已经具有剔除了偶数的功能
•         if(0== num[i])
•         {
•             continue;
•         }
•         //从i的2倍开始过滤
•         for(j= i + i; j<= n;j+=i)
•         {
•             //是直接进行赋值快呢？还是在此处加上判断快呢？？不晓得啊？求解。。
•             if( j % 2 == 0)
•             {
•                 continue;
•             }
•             else
•             {
•                 num[j]= 0;
•             }
•         }
•     }

   第一部分红色，我将原来的奇数和偶数进行判断，变为了只对奇数进行判断；

   第二部分红色，我将奇数的倍数为偶数的直接剔除，变成只对倍数为奇数的进行赋值；

   以上两者改变，都基于开始时，已经将偶数剔除。

   对NUM = 300000测试如下所示。

   

   时间仅为7毫秒，比 优化前NUM = 300000时，时间更快。

 

6、继续优化

   

   C语言实现代码如下所示。

点击(此处)折叠或打开

• //合并筛选法的优化方案
• #include <stdio.h>
• #include <time.h>
• #include <math.h>
• #include <malloc.h>
• #include <string.h>
• #define NUM 10000
• int Test_prime(int n)
• {
•     int i, j;
•     // 素数数量统计 
•     int count= 0;
•     //分配素数标记空间，明白+1原因了吧，因为浪费了一个num[0]
•     char *num= (char*)malloc( n+1);
•     //干嘛用的，请仔细研究下文
•     int mpLen= 2*3*5*7*11*13;
•     char magicPattern[2*3*5*7*11*13];
•     // 奇怪的代码，想！ 
•     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= num+1;
•     char* pstop= p+n-remainder;
•     while(p < pstop)
•     {
•         memcpy(p, magicPattern, mpLen);
•         p += mpLen;
•     }
•     if (remainder > 0)
•     {
•         memcpy(p, magicPattern, remainder);
•     }
•     num[2]= 1;
•     num[3]= 1;
•     num[5]= 1;
•     num[7]= 1;
•     num[11]= 1;
•     num[13]= 1;
•     
•     // 从17开始过滤，因为2,3,5,7,11,13的倍数早被去掉了 
•     // 到n/13止的
•     int stop= n/13;
•     for (i=17; i<= stop; i++)
•     {
•         // i是合数
•         if(0 == num[i]) 
•         {
•             continue;
•         }
•         
•         //从i的17倍开始过滤
•         intstep = i*2;
•         for(j=i*17; j<= n; j+=step)
•         {
•             num[j]= 0;
•         }
•     }
•     
•     // 统计素数个数
•     for (i=2; i<=n; i++)
•     {
•         if(num[i]) 
•         {
•             count++;
•         }
•     }
•     
•     // 释放内存
•     free(num);
•     return count;
• }
• int main()
• {
•     int count;
•     clock_t start,end;
•     start = clock();
•     count = Test_prime(NUM);
•     end = clock();
•     printf("%d内的素数个数为：%d, 总共耗时为：%d毫秒\n", NUM, count,end - start);
•     
•     return 0;
• }

   测试结果如下所示。

   说实话，这种思想真的很赞，现在的我是无法想到的，感谢作者，让我有了更广泛的见识。

7、其他

   除了以上几种算法外，如拉宾米勒素数测试算法，感觉这个算法比较难，先好好看看，等弄懂了，然后补上。