『转』求N内的所有素数

转自梦醒潇湘       http://blog.chinaunix.net/uid-26548237-id-3364131.html

先收藏起来,往后慢慢研究。。

*****************************************以下为正文******************************************************

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

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

  

1、素数及相关

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

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

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

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

2、试除法求素数

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

 

   C语言实现如下所示。

  1. #include

  2. #include

  3. #include

  4. #include

  5. //试除法

  6. #define NUM 10000

  7. int Test_prime(int n)

  8. {

  9.     int count = 0;

  10.     int i, j;

  11.     int *num= (int*)malloc(sizeof(int)* n);

  12.     num[count++]= 2;

  13.     

  14.     for(i = 3; i <= n; i++)

  15.     {

  16.         for(j= 2; j <= sqrt(i); j++)

  17.         {

  18.             if(i% j == 0)

  19.             {

  20.                 break;

  21.             }

  22.         }

  23.         if(j> sqrt(i))

  24.         {

  25.             num[count++]= i;

  26.         }

  27.     }

  28.     free(num);

  29.     return count;

  30. }

  31. int main()

  32. {

  33.     int count;

  34.     clock_t start,end;

  35.     start = clock();

  36.     count = Test_prime(NUM);

  37.     end = clock();

  38.     printf("%d 内的素数个数为:%d, 总共耗时为:%d 毫秒\n", NUM, count,end - start);

  39.     return 0;

  40. }

   测试结果如下所示(测试在VC6.0下进行)。

26548237_1349254349Rjz7.jpg

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

3、试除法的优化方案

  

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

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

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

   C语言实现如下所示。

点击(此处)折叠或打开

  1. //求N内的所有素数
  2. #include

  3. #include

  4. #include

  5. #include

  6. //试除法

  7. #define NUM 1000000

  8. int Test_prime(int n)

  9. {

  10.     int count = 0;

  11.     int i, j, k, stop;

  12.     //分配空间

  13.     int *num= (int*)malloc(sizeof(int)* n);

  14.     //2肯定是素数

  15.     num[count++]= 2;

  16.     stop = 0;

  17.     for(i = 3; i <= n; i++)

  18.     {

  19.         //在循环中重复调用sqrt是低效做法,故引入k

  20.         k = (int)sqrt(i);

  21.         //stop的作用是:统计小于当前k值的质数的数目

  22.         while(num[stop]<= k && stop < count)

  23.         {

  24.             stop++;

  25.         }

  26.         for(j= 0; j < stop; j++)

  27.         {

  28.             if( i% num[j]== 0)

  29.             {

  30.                 //i不能被2-sqrt(i)间的素数整除,自然不能被其他整数整除,所以为素数

  31.                 break;

  32.             }

  33.         }

  34.         if(j== stop)

  35.         {

  36.             num[count++]= i;

  37.         }

  38.        

  39.     }

  40.     free(num);

  41.     return count;

  42. }

  43. int main()

  44. {

  45.     int count;

  46.     clock_t start,end;

  47.     start = clock();

  48.     count = Test_prime(NUM);

  49.     end = clock();

  50.     printf("%d 内的素数个数为:%d, 总共耗时为:%d 毫秒\n", NUM, count,end - start);

  51.    

  52.     return 0;

  53. }

  测试结果如下所示。

26548237_1349259935jEnd.jpg

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

4、合数过滤筛选法

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

   C语言实现如下所示。

  1. //合并筛选法
  2. #include

  3. #include

  4. #include

  5. #include

  6. //试除法

  7. #define NUM 10000

  8. int Test_prime(int n)

  9. {

  10.     int count = 0;

  11.     int i, j;

  12.     //分配空间,之所以是n+1,是因为浪费了一个num[0]

  13.     char *num =(char *)malloc(sizeof(char)* (n + 1));

  14.     //初始化素数标记

  15.     for(i = 2; i<= n; i++)

  16.     {

  17.         num[i]= 1;

  18.     }

  19.     //以2-(N-1)为因子过滤合数

  20.     for(i = 2; i <= n-1; i++)

  21.     {

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

  23.         {

  24.             //i*j是由两整数相乘而得,显然不是素数

  25.             num[i*j]= 0;

  26.         }

  27.     }

  28.     //统计素数个数

  29.     for( i = 2; i<= n; i++)

  30.     {

  31.         if( 1== num[i])

  32.         {

  33.             count++;

  34.         }

  35.     }

  36.     free(num);

  37.     return count;

  38. }

  39. int main()

  40. {

  41.     int count;

  42.     clock_t start,end;

  43.     start = clock();

  44.     count = Test_prime(NUM);

  45.     end = clock();

  46.     printf("%d 内的素数个数为:%d, 总共耗时为:%d 毫秒\n", NUM, count,end - start);

  47.    

  48.     return 0;

  49. }

   测试结果如下所示。

26548237_1349259935jEnd.jpg

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

5、合并筛选法优化方案

   上述算法存在的问题是:

   1> 在外层循环,需要一直执行到n-1嘛?不要,因为n/2-(n-1)之间的数显然不能整除出n;

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

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

  

   C语言实现如下所示。

点击(此处)折叠或打开

  1. //合并筛选法的优化方案

  2. #include

  3. #include

  4. #include

  5. #include

  6. #define NUM 300000

  7. int Test_prime(int n)

  8. {

  9.     int count = 0;

  10.     int i, j;

  11.     //分配空间

  12.     char *num =(char *)malloc(sizeof(char)* (n + 1));

  13.    

  14.     //初始化素数标记

  15.     num[2] = 1;

  16.     //注意此处是i<n,上例中的i<=n

  17.     for(i = 3; i< n; i++)

  18.     {

  19.         num[i++]= 1;

  20.         num[i]= 0;//偶数自然不是素数

  21.     }

  22.     //如果n为奇数

  23.     if(n % 2 != 0)

  24.     {

  25.         num[n]= 1;

  26.     }

  27.     //从3开始过滤,因为,2的倍数在初始化中去掉了

  28.    for(i = 3; i <= n/2; i++)

  29.     {

  30.         if(0 == num[i] )

  31.         {

  32.             continue;

  33.         }

  34.         //从i的2倍开始过滤

  35.         for(j = i + i; j <= n;j+=i)

  36.         {

  37.             num[j] = 0;

  38.         }

  39.     }

  40.     //统计素数个数

  41.     for( i = 2; i<= n; i++)

  42.     {

  43.         if( 1== num[i])

  44.         {

  45.             count++;

  46.         }

  47.     }

  48.     free(num);

  49.     return count;

  50. }

  51. int main()

  52. {

  53.     int count;

  54.     clock_t start,end;

  55.     start = clock();

  56.     count = Test_prime(NUM);

  57.     end = clock();

  58.     printf("%d 内的素数个数为:%d, 总共耗时为:%d 毫秒\n", NUM, count,end - start);

  59.    

  60.     return 0;

  61. }

   测试如下所示。

26548237_1349262322scAj.jpg

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

  

   后来想到进行添加补充:

 

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

点击(此处)折叠或打开

  1. //从3开始过滤,因为,2的倍数在初始化中去掉了

  2.     for(i = 3; i <= n/2;i = i + 2)

  3.     {

  4.         //在这里进行判断,就已经具有剔除了偶数的功能

  5.         if(0== num[i])

  6.         {

  7.             continue;

  8.         }

  9.         //从i的2倍开始过滤

  10.         for(j= i + i; j<= n;j+=i)

  11.         {

  12.             //是直接进行赋值快呢?还是在此处加上判断快呢??不晓得啊?求解。。

  13.            if( j % 2 == 0)

  14.             {

  15.                 continue;

  16.             }

  17.             else

  18.             {

  19.                 num[j]= 0;

  20.             }

  21.         }

  22.     }

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

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

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

   对NUM = 300000测试如下所示。

26548237_1349275278v22o.jpg

  

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

6、继续优化

  

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

点击(此处)折叠或打开

  1. //合并筛选法的优化方案

  2. #include

  3. #include

  4. #include

  5. #include

  6. #include

  7. #define NUM 10000

  8. int Test_prime(int n)

  9. {

  10.     int i, j;

  11.     // 素数数量统计

  12.     int count = 0;

  13.     // 分配素数标记空间,明白+1原因了吧,因为浪费了一个num[0]

  14.     char *num =(char*)malloc( n+1);

  15.     // 干嘛用的,请仔细研究下文

  16.     int mpLen = 2*3*5*7*11*13;

  17.     char magicPattern[2*3*5*7*11*13];

  18.     // 奇怪的代码,想!

  19.     for (i=0; i<mpLen; i++)

  20.     {

  21.         magicPattern[i++]= 1;

  22.         magicPattern[i++]= 0;

  23.         magicPattern[i++]= 0;

  24.         magicPattern[i++]= 0;

  25.         magicPattern[i++]= 1;

  26.         magicPattern[i]= 0;

  27.     }

  28.     for (i=4; i<=mpLen; i+=5)

  29.     {

  30.         magicPattern[i]= 0;

  31.     }

  32.     for (i=6; i<=mpLen; i+=7)

  33.     {

  34.         magicPattern[i]= 0;

  35.     }

  36.     for (i=10; i<=mpLen; i+=11)

  37.     {

  38.         magicPattern[i]= 0;

  39.     }

  40.     for (i=12; i<=mpLen; i+=13)

  41.     {

  42.         magicPattern[i]= 0;

  43.     }

  44.    

  45.     // 新的初始化方法,将2,3,5,7,11,13的倍数全干掉

  46.     // 而且采用memcpy以mpLen长的magicPattern来批量处理

  47.     int remainder = n%mpLen;

  48.     char* p = num+1;

  49.     char* pstop = p+n-remainder;

  50.     while (p< pstop)

  51.     {

  52.         memcpy(p, magicPattern, mpLen);

  53.         p += mpLen;

  54.     }

  55.     if (remainder> 0)

  56.     {

  57.         memcpy(p, magicPattern, remainder);

  58.     }

  59.     num[2] = 1;

  60.     num[3] = 1;

  61.     num[5] = 1;

  62.     num[7] = 1;

  63.     num[11]= 1;

  64.     num[13]= 1;

  65.    

  66.     // 从17开始过滤,因为2,3,5,7,11,13的倍数早被去掉了

  67.     // 到n/13止的

  68.     int stop = n/13;

  69.     for (i=17; i<= stop; i++)

  70.     {

  71.         // i是合数

  72.         if (0== num[i])

  73.         {

  74.             continue;

  75.         }

  76.        

  77.         // 从i的17倍开始过滤

  78.         int step= i*2;

  79.         for (j=i*17; j<= n; j+=step)

  80.         {

  81.             num[j]= 0;

  82.         }

  83.     }

  84.    

  85.     // 统计素数个数

  86.     for (i=2; i<=n; i++)

  87.     {

  88.         if (num[i])

  89.         {

  90.             count++;

  91.         }

  92.     }

  93.    

  94.     // 释放内存

  95.     free(num);

  96.    

  97.     return count;

  98. }

  99. int main()

  100. {

  101.     int count;

  102.     clock_t start,end;

  103.     start = clock();

  104.     count = Test_prime(NUM);

  105.     end = clock();

  106.     printf("%d 内的素数个数为:%d, 总共耗时为:%d 毫秒\n", NUM, count,end - start);

  107.    

  108.     return 0;

  109. }

   测试结果如下所示。

26548237_13492661023IWp.jpg

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

7、其他

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

  

   通过今天的纠结,对于求素数有了更加深刻的了解和认识,感觉自己还差很多,需要更加的努力。

   另外,感谢来自百度空间的作者 doforfun_net,给我了很大的启发,学到了很多。

                                                      梦醒潇湘

                                                  2012-10-3 20:15

 

本文章迁移自http://blog.csdn.net/timberwolf_2012/article/details/8066111