半秒内筛一亿以内的所有素数_筛一亿质数 时间-程序员宅基地

技术标签: C  数据结构  素数  

【问题描述】:

   试编写一个程序,找出2->N之间的所有质数。希望用尽可能快的方法实现。
【问题分析】:
   这个问题可以有两种解法:一种是用“筛子法”,另一种是“除余法”。
   如果要了解“除余法”,请看另一篇文章连接 http://blog.csdn.net/lingling_1/article/details/44459423 《求质数 之 除余法(C语言描述)》。
   这里我们来讨论一下用“筛法”来解决这个问题。
   先来举个简单的例子来介绍一下“筛法”,求2~20的质数,它的做法是先把2~20这些数一字排开:
   2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
   先取出数组中最小的数,是2,则判断2是质数,把后面2的倍数全部删掉。
   2 | 3 5 7 9 11 13 15 17 19
   接下来的最小数是3,取出,再删掉3的倍数
   2 3 | 5 7 11 13 17 19
   一直这样下去,直到结束。剩下的数都是素数。
   筛法的原理是:
   1.数字2是素数。
   2.在数字K前,每找到一个素数,都会删除它的倍数,即以它为因子的整数。如果k未被删除,就表示2->k-1都不是k的因子,那k自然就是素数了。
   (1)除余法那篇文章里也介绍了,要找出一个数的因子,其实不需要检查2->k,只要从2->sqrt(k),就可以了。所有,我们筛法里,其实只要筛到sqrt(n)就已经找出所有的素数了,其中n为要搜索的范围。
   (2)另外,我们不难发现,每找到一个素数k,就一次删除2k, 3k, 4k,..., ik,不免还是有些浪费,因为2k已经在找到素数2的时候删除过了,3k已经在找到素数3的时候删除了。因此,当i<k时,都已经被前面的素数删除过了,只有那些最小的质因子是k的那些数还未被删除过,所有,就可以直接从k*k开始删除。
   (3)再有,所有的素数中,除了2以外,其他的都是奇数,那么,当i时奇数的时候,ik就是奇数,此时k*k+ik就是个偶数,偶数已经被2删除了,所有我们就可以以2k为单位删除步长,依次删除k*k, k*k+2k, k*k+4k, ...。
   (4)我们都清楚,在前面一小段范围内,素数是比较集中的,比如1->100之间就有25个素数。越到后面就越稀疏。
   因为这些素数本身值比较小,所以搜索范围内,大部分数都是它们的倍数,比如搜索1->100,这100个数。光是2的倍数就有50个,3的倍数有33个,5的倍数20个,7的倍数14个。我们只需搜索到7就可以,因此一共做删除操作50+33+20+14=117次,而2和3两个数就占了83次,这未免太浪费时间了。
   所以我们考虑,能不能一开始就排除这些小素数的倍数,这里用2和3来做例子。
   如果仅仅要排除2的倍数,数组里只保存奇数:1、3、5...,那数字k的坐标就是k/2。
   如果我们要同时排除2和3的倍数,因为2和3的最小公倍数是6,把数字按6来分组:6n, 6n+1, 6n+2, 6n+3, 6n+4, 6n+5。其中6n, 6n+2, 6n+4是2的倍数,6n+3是3的倍数。所以数组里将只剩下6n+1和6n+5。n从0开始,数组里的数字就一次是1, 5, 7, 11, 13, 17...。
   现在要解决的问题就是如何把数字k和它的坐标i对应起来。比如,给出数字89,它在数组中的下标是多少呢?不难发现,其实上面的序列,每两个为一组,具有相同的基数n,比如1和5,同是n=0那组数,6*0+1和6*0+5;31和35同是n=5那组,6*5+1和6*5+5。所以数字按6分组,每组2个数字,余数为5的数字在后,所以坐标需要加1。
   所以89在第89/6=14组,坐标为14*2=28,又因为89%6==5,所以在所求的坐标上加1,即28+1=29,最终得到89的坐标i=29。同样,找到一个素数k后,也可以求出k*k的坐标等,就可以做筛法了。
   这里,我们就需要用k做循环变量了,k从5开始,交替与2和4相加,即先是5+2=7,再是7+4=11,然后又是11+2=13...。这里我们可以再设一个变量gab,初始为4,每次做gab = 6 - gab,k += gab。让gab在2和4之间交替变化。另外,2和4都是2的幂,二进制分别为10和100,6的二进制位110,所以可以用k += gab ^= 6来代替。参考代码:
gab = 4;
for (k = 5; k * k <= N; k += gab ^= 6)
{
...
}
   但我们一般都采用下标i从0->x的策略,如果用i而不用k,那应该怎么写呢?
   由优化策略(1)可知,我们只要从k2开始筛选。n=i/2,我们知道了i对应的数字k是素数后,根据(2),那如何求得k2的坐标j呢?这里假设i为偶数,即k=6n+1。
   k2 = (6n+1)*(6n+1) = 36n2 + 12n + 1,其中36n2+12n = 6(6n2+2n)是6的倍数,所以k2除6余1。
   所以k2的坐标j = k2/6*2 = 12n2+4n。
   由优化策略(2)可知,我们只要依次删除k2+2l×k, l = 0, 1, 2...。即(6n+1)×(6n+1+2l)。
   我们发现,但l=1, 4, 7...时,(6n+1+2l)是3的倍数,不在序列中。所以我们只要依次删除k2, k2+4l, k2+4l+2l...,又是依次替换2和4。
   为了简便,我们可以一次就删除k2和k2+4l两项,然后步长增加6l。所以我们需要求len=4l和stp=6l。不过这里要注意一点,k2+4k=(6n+1)*(6n+5),除以6的余数是5,坐标要加1。
   len = k*(k+4)/6*2 - k2/6*2 = (6n+1)*(6n+1+4)/6*2+1 - (6n+1)*(6n+1)/6*2 = (12n2+12n+1) - (12n2+4n) = 8n+1;
   stp = k*(k+6)/6*2 - k2/6*2 = 12n+2;
   最终,我们得到:
   len = 8n+1;
   stp = 12n+2;
    j = 12n2+4n;
   同理可以求出k=6n+5时的情况:
   len = 4n+3;
   stp = 12n+10;
    j = 12n2+20n+8;
   下面的代码在实现上用了位运算,可能有点晦涩。
★注:第5种优化方法还是理论阶段,下面的代码中并未采用这种优化算法,仅供大家参考。
   (5)由(2)可知,如果每找到一个素数k,能依次只删除以k为最小素数因子的数,那么每个数字就都只被删除一次,那这个筛法就能达到线性的O(n)效率了。比如数字600 = 2*2*3*5*11,其中2是它的最小素数因子。那这个数就被2删除了。3、5、11虽然都是它的因子,但不做删除它的操作。要实现这种策略,那每找到一个素数k,那从k开始,一次后面未被删除的数字来与k相乘,删除它们的积。比如要筛出2~60之间的素数:
   1.先列出所有的数。
   2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... 50 51 52 53 54 55 56 57 58 59 60
   2.选出序列中的第一个数,即2,判断它是素数,然后从2开始,依次与剩下的未被删除的数相乘,删除它们的积。即2*2=4, 2*3=6,2*4=8...。
   2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... 50 51 52 53 54 55 56 57 58 59 60
   02 | 03 05 07 09 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59
   3.去掉2后,再选出序列中第一个数,即3,判断它是素数,然后从3开始,依次与剩下的数相乘,即3*3=9,3*5=15,3*7=21...
   02 | 03 05 07 09 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59
   02 03 | 05 07 11 13 15 17 19 23 25 29 31 35 37 41 43 47 49 53 55 59
   4.去掉3后,选出最小的数5,为素数,依次删除5*5=25,5*7=35,5*11=55,...
   02 03 | 05 07 11 13 15 17 19 23 25 29 31 35 37 41 43 47 49 53 55 59
   02 03 05 | 07 11 13 15 17 19 23 29 31 37 41 43 47 49 53 59
   5.去掉5后,选出最小的数7,为素数,删除7*7=49,...
   02 03 05 | 07 11 13 15 17 19 23 29 31 37 41 43 47 49 53 59
   02 03 05 | 07 11 13 15 17 19 23 29 31 37 41 43 47 53 59
   6.去掉7后,第一个数11的平方121大于60,所以结束。剩下的数字全为素数。
   02 03 05 07 11 13 15 17 19 23 29 31 37 41 43 47 53 59 |
   上面的操作效率很高,但在计算机中模拟的时候却又很大的障碍:
   首先,计算机内存是一维的空间,很多时候我们不能随心所欲,要实现上面的算法,要求这个数据结构既能很高效地查找某个特定的值,又能不费太大代价对序列中的元素进行删除。高效地查找,用数组是最合适的了,能在O(1)的时间内对内存进行读写,但要删除序列中一个元素却要O(n);单链表可以用O(1)的时间做删除操作,当然要查找就只能是O(n)了。所以这个数据结构很难找。
   其次,筛法的一个缺点就是空间浪费太大,典型的以空间换时间。如果我们对数组进行压缩,比如初始时就排除了所有偶数,数组0对应数字1,1对应3,...。这样又会因为多了一道计算数字下标的工序而浪费时间。这又是一个矛盾的问题。
   也许我们可以试试折中的办法:数据结构综合数组和链表2种,数组用来做映射记录,链表来记录剩下的还未被删除的数据,而且开始也不必急着把链表里的节点释放掉,只要在数组里做个标记就可以了。下次遍历到这个数字时才删除。这样为了删除,可以算只遍历了一次链表,不过频繁地使用free()函数,也许又会减低效率。总之,我们所做的,依然是用空间来换时间,记录更多的信息,方便下次使用,减少再次生成信息所消耗的时间。
【程序清单】:
  1. #include <time.h>  
  2. #include <stdio.h>  
  3. #define N 100000000  
  4. #define size (N/6*2 + (N%6 == 5? 2: (N%6>0)))  
  5. int p[size / 32 + 1] = {1};  
  6. int creat_prime(void)  
  7. {  
  8.     int i, j;  
  9.     int len, stp;  
  10.     int c = size + 1;  
  11.     for (i = 1; ((i&~1)<<1) * ((i&~1) + (i>>1) + 1) < size; i++)  
  12.     {  
  13.         if (p[i >> 5] >> (i & 31) & 1) continue;  
  14.         len = (i & 1)? ((i&~1)<<1) + 3: ((i&~1)<<2) + 1;  
  15.         stp = ((i&~1)<<1) + ((i&~1)<<2) + ((i & 1)? 10: 2);  
  16.         j = ((i&~1)<<1) * (((i&~1)>>1) + (i&~1) + 1) + ((i & 1)? ((i&~1)<<3) + 8 + len: len);  
  17.         for (; j < size; j += stp)  
  18.         {  
  19.             if (p[j >> 5] >> (j & 31) & 1 ^ 1)  
  20.                 p[j >> 5] |= 1L << (j & 31), --c;  
  21.             if (p[(j-len) >> 5] >> ((j-len) & 31) & 1 ^ 1)  
  22.                 p[(j-len) >> 5] |= 1L << ((j-len) & 31), --c;  
  23.         }  
  24.         if (j - len < size && (p[(j-len) >> 5] >> ((j-len) & 31) & 1 ^ 1))  
  25.             p[(j-len) >> 5] |= 1L << ((j-len) & 31), --c;  
  26.     }  
  27.     return c;  
  28. }  
  29. int main( )  
  30. {  
  31.     clock_t t = clock();  
  32.     printf("%d ", creat_prime());  
  33.     printf("Time: %f ", 1.0 * (clock() - t) / CLOCKS_PER_SEC);  
  34. }  
#include <time.h>
#include <stdio.h>
#define N 100000000
#define size (N/6*2 + (N%6 == 5? 2: (N%6>0)))
int p[size / 32 + 1] = {1};
int creat_prime(void)
{
    int i, j;
    int len, stp;
    int c = size + 1;
    for (i = 1; ((i&~1)<<1) * ((i&~1) + (i>>1) + 1) < size; i++)
    {
        if (p[i >> 5] >> (i & 31) & 1) continue;
        len = (i & 1)? ((i&~1)<<1) + 3: ((i&~1)<<2) + 1;
        stp = ((i&~1)<<1) + ((i&~1)<<2) + ((i & 1)? 10: 2);
        j = ((i&~1)<<1) * (((i&~1)>>1) + (i&~1) + 1) + ((i & 1)? ((i&~1)<<3) + 8 + len: len);
        for (; j < size; j += stp)
        {
            if (p[j >> 5] >> (j & 31) & 1 ^ 1)
                p[j >> 5] |= 1L << (j & 31), --c;
            if (p[(j-len) >> 5] >> ((j-len) & 31) & 1 ^ 1)
                p[(j-len) >> 5] |= 1L << ((j-len) & 31), --c;
        }
        if (j - len < size && (p[(j-len) >> 5] >> ((j-len) & 31) & 1 ^ 1))
            p[(j-len) >> 5] |= 1L << ((j-len) & 31), --c;
    }
    return c;
}
int main( )
{
    clock_t t = clock();
    printf("%d ", creat_prime());
    printf("Time: %f ", 1.0 * (clock() - t) / CLOCKS_PER_SEC);
}


【运行结果】:
5761455
Time: 0.300000
运行环境:Linux debian 2.6.26-1-686、GCC (Debian 4.3.2-1.1) 4.3.2
【算法比较】:
   现在,我们已经拥有初步改进的“筛法”和“除余法”的函数了,把它们加到自己的函数库里。方便下次调用。
   这里,我想说一下个人对这两种算法的使用经验:
   就时间效率上讲,筛法绝对比除余法高。比如上面的代码,可以在半秒内筛一亿以内的所有素数。如果用除余法来解决这样的问题,绝对可以考验一个人的耐性。因此,在搜索空间比较大的时候,“筛法”无疑会是首选。
   但筛法是以空间换时间,用除余法,我们只要开一个可以容纳结果的数组就可以了,而筛法开的数组要求可以容纳整个搜索范围;另外,我们用“除余法”得到的结果,是一个已经排好序的素数序列,如果要解决的问题需要用到这些连续的素数,而且搜索范围也不大,那显然除余法很适合。而“筛法”得到的结果,是一个布尔型的表格,通过它,你可以很轻松的判断某个数是不是素数,但如果你想知道这个素数的下一个素数是多大,可能要费点劲了。

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/lingling_1/article/details/44462185

智能推荐

【ICPC济南区域赛】签到题题解_icpc2021济南题解-程序员宅基地

文章浏览阅读406次。M Cook Pancakes!对于N个饼,显然需要fry 2N次;考虑最优的情况,显然在每次都同时fry K次即答案为 【2N/K】#include<bits/stdc++.h>using namespace std; int main(){ int N,K; cin>>N>>K; if(N<=K)cout<<2<<endl; else{ cout<<ceil(2_icpc2021济南题解

大数据分析模型构建步骤_构建大数据分析模型-程序员宅基地

文章浏览阅读4.1k次。  我们知道做认识事情都有个流程顺序,正确的流程可以事半功倍,错误的流程往往会导致事情重新来做,越来越多的企业都实现了大数据营销推广。今天我们就来了解一下,大数据分析中的模型构建步骤。  大数据环境下的数据分析模型构建步骤  常用的数据挖掘方法主要是基于客户画像体系与结果,选取相关性较大的特征变量,通过分类模型、聚类模型、回归模型、神经网络和关联规则等机器算法进行深度挖掘。常用算法的基本内容如下:  1、分类和聚类  分类算法是极其常用的数据挖掘方法之一,其核心思..._构建大数据分析模型

Xdebug安装与使用_xdebug 验证脚本的目标目录不是 public。-程序员宅基地

文章浏览阅读1.1k次。Xdebug安装与使用为什么需要Debugger?很多PHP程序员调试使用echo、print_r()、var_dump()、printf()等,其实对 于有较丰富开发经验的程序员来说这些也已经足够了,他们往往可以在程序执行的过程中,通过输出特定变量的值可以判断程序执行是否正确,甚至效率高低也可以 看出来(当然可能还需要使用一些时间函数)。那么我们为什么还需要一个专门的_xdebug 验证脚本的目标目录不是 public。

Java配合Tabula框架实现上传并解析PDF表格_java tabula-程序员宅基地

文章浏览阅读1.5k次,点赞3次,收藏6次。功能:解析上传的pdf表格,并存入数据库最近有个需求,功能如上。百度了一下主要推荐的框架有两个。一个是Itext,听说很厉害,但是商业使用需要花钱就没有太多了解。另一个是PdfBox,简单的写了个demo,可以获取到pdf内的所有文字并返回String。返回的数据位置会错乱,且api没有中文版,例子也不多,使用起来很麻烦(是我太菜)。在之后发现了Tabula,功能实现的很强大,就搜了搜看有没有给Java调用的方法,就找到了下边的贴子:https://blog.csdn.net/qq_3695600_java tabula

使用VsCode打造C#开发IDE_vscode 做c#开发-程序员宅基地

文章浏览阅读6.6k次,点赞3次,收藏25次。用VsCode写了几天Java,还是比较满意的,无论是在智能提示方面,还是在调试跳转文本编辑等方面,个人感觉都不次于IDEA等正牌重型IDE,所以就想顺带用VsCode把VisualStudio也替代了,但是还是发现有点儿小问题,就是必须严格的按照VsCode的新建步骤来,否则运行调试时会报错误。下面就详细分享一下用VsCode开发C#的步骤。_vscode 做c#开发

mysql-数据库字段date datetime timestamp与实体类类型对应关系_数据库的时间与实体类型的关系-程序员宅基地

文章浏览阅读8k次,点赞12次,收藏36次。https://www.cnblogs.com/lrzr/archive/2017/08/07/7299211.htmlhttps://blog.csdn.net/weixin_38336276/article/details/83892408https://www.cnblogs.com/1130136248wlxk/articles/5238538.html_数据库的时间与实体类型的关系

随便推点

存储过程 @与字符连接_存储过程连接符-程序员宅基地

文章浏览阅读289次。&gt; BEGIN -&gt; declare a int; -&gt; declare b varchar(5000); -&gt; set a=1; -&gt; set b=''; -&gt; while a&lt;10 do -&gt; set b = concat(b,',',a); -&gt; set a=a+_存储过程连接符

实验三 顺序图、协作图设计_添加课程顺序图-程序员宅基地

文章浏览阅读6.1k次,点赞6次,收藏36次。实验三 顺序图、协作图设计【实验目的】理解顺序图和协作图的概念及作用; 掌握UML顺序图与协作图的基本图形,了解它们各自的组成元素、特定作用和适用场合; 重点掌握顺序图的画法及其中元素所代表的意义。【实验性质】设计性实验。【实验要求】学习根据指定的用例描述绘制顺序图和协作图的方法; 学习使用Rational Rose绘制顺序图和协作图; 掌握顺序图和协作图的相互转换方法。【实验内容】以网上选课系统中的Select Course(选课)用例为例,设计和实现顺序图、协作图.._添加课程顺序图

Kaggle滑水 - CTR预估(FM_FFM)_ffm在滑坡预测里面是什么-程序员宅基地

文章浏览阅读2.4k次。本文继续以Avazu-CTR赛题为背景,尝试采用FM(Factorization Machine,因子分解机)及FFM(Field-aware Factorization Machine,场感知因子分解机)来进行CTR预估任务。本文的源码托管于我的Github:PnYuan - Kaggle_CTR,欢迎查看交流。1.概念商用推荐场景中的CTR预估工作易面临大规模稀疏数据的挑战..._ffm在滑坡预测里面是什么

安装numpy+mkl报错的处理办法_numpy mkl后报错-程序员宅基地

文章浏览阅读2.6k次。转载自:http://www.fkccp.com/archives/2710.html 非常感谢这个大神,我纠结了好久这个问题!Processing c:\users\rao\downloads\numpy-1.11.2+mkl-cp27-cp27m-win32.whl Installing collected packages: numpy Exception: Traceback (mo_numpy mkl后报错

POJ 2746 约瑟夫问题_poj2746-程序员宅基地

文章浏览阅读525次。题目总时间限制: 1000ms 内存限制: 65536kB 描述 约瑟夫问题:有n只猴子,按顺时针方向围成一圈选大王(编号从1到n),从第1号开始报数,一直数到m,数到m的猴子退出圈外,剩下的猴子再接着从1开始报数。就这样,直到圈内只剩下一只猴子时,这个猴子就是猴王,编程求输入n,m后,输出最后猴王的编号。输入 每行是用空格分开的两个整数,第一个是 n, 第二个是 m ( 0 < m,n <_poj2746

如何为微信小程序添加定位导航和地图标注功能_小程序地图导航-程序员宅基地

文章浏览阅读1.3k次。在JavaScript代码中,我们需要将获取到的用户位置信息设置为map组件的latitude和longitude属性,并创建一个标注点并设置在markers属性中。在获取到用户的地理位置信息后,我们需要将用户的位置在地图上进行标注。在JavaScript代码中,我们需要在数据中添加路线规划相关信息,并在路线规划成功的回调函数中更新数据。同样的,我们在实现地图标注功能时,也需要获取用户的地理位置信息。至此,我们已经成功地实现了获取用户地理位置信息并在地图上标注出用户位置的功能。二、地图标注功能的实现。_小程序地图导航

推荐文章

热门文章

相关标签