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
