首页 >

埃氏筛法搜寻1亿以内素数

作者:fortranboy  日期:10-09
来源:Fcode研讨团队
埃氏筛法搜寻1亿以内素数,标记、统计并输出个数和耗时。网络类似算法不少,但fortran版本未见(至少是中文网站上)。本代码部分减少了重复标记,效率较高。
详见代码注释。


  ! szw_sh@163.com
  ! 2018-10-07
  ! 标记2到1亿所有素数,统计个数
  ! 算法:用埃氏筛法;素数分布特征;部分减少重复标记
  ! 结果:5761455     680 ms
  ! 效率:3.2GHz主频,2G内存,运行时间0.8秒
  ! 编译:FPS4 /Ox ; CVF6 默认选项
Program Main  
  integer(4) :: i,j,k,m,n,d
  integer(4),parameter :: mx=100000000,mx1=mx*1.1
  integer(1) :: a(mx1)      ! mx1拓展数组,防止写溢出
  x=timef()                 ! 计时开始
  a=0                       ! a标记对应下标自然数,1为素数,否则标0
  a(2:3)=1                  ! 2,3 单独标记
  m=2
  n=sqrt(mx*1.0)            ! 搜索到mx平方根。见相关定理                                        
  a(6-1:mx-1:6)=1           ! 以6为步长,初始化标记左右侧为素数
  a(6+1:mx+1:6)=1           ! 6i,6i+2,6i+3,6i+4 等皆为合数  
  do i=6,n,6                ! 以6为步长,搜索素数
    do j=i-1,i+1,2          ! 搜索循环节点的左右2个数
      if(a(j)==0) cycle
      n=j*j
      d=j*2
      do k=n,mx,d           ! 标记搜到的素数的倍数为合数
        if(mod(k,6)==3) exit! 标记从j*j开始,更小的已被前一个素数倍数标记
        a(k)=0              ! 可以j*2为步长,跳过偶数
      end do                ! 遇到余数3时跳出循环,后面只对余数1或5的情形作标记
      n=k+d                 ! 循环标记倍数的起始点
      do k=n,mx,d*3         ! 按照j*2*3步长标记素数的倍数
        if(a(k)/=0) a(k)=0  ! 变化周期为余数5,1,3,省却一次标记动作
        if(a(k+d)/=0) a(k+d)=0! 使用if判断比重复标记更快
      end do
    end do
  end do
  do i=6,mx,6               ! 统计素数个数
    m=m+a(i-1)+a(i+1)
  end do
  n=timef()*1000            ! 计时结束
  write(*,'(1x,2i8,a)') m,n,' ms'! 输出个数,时间
End Program Main
常规|工具|专业|读物|
代码|教学|算法|
首页 >
FortranCoder手机版-导航