Fortran 埃氏筛法统计1亿以内素数的个数

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'! 输出个数，时间

endprogram main

 

! 2021.10.17

! 标记2到1亿所有素数，统计个数

! 算法：用埃氏筛法；素数分布特征；部分减少重复标记

! 结果：5761455 997 ms

! 硬件：I5-10210U 12G内存

! 编译：Windows10 IVF 2011