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

埃氏筛法搜寻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