! This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/primes.f90
!
!  by James Van Buskirk
program sieve0
   implicit none
   integer, parameter :: dp = selected_real_kind(15,300)
   integer, allocatable :: table(:)
   integer N
   integer, parameter :: nbits = bit_size(1)
   integer i
   integer j
   integer :: np = 0
   character(80) filename

   write(*,'(a)',advance='no') ' Enter max N:> '
   read(*,*) N
   allocate(table(0:N/nbits))
   table = not(0)
   table(0) = not(3)
   do i = 0, floor(sqrt(size(table)*nbits-0.5_dp))
      if(iand(table(i/nbits), ishft(1, mod(i,nbits))) == 0) cycle
      do j = i**2, nbits*size(table)-1, i
        table(j/nbits) = iand(table(j/nbits),not(ishft(1,mod(j,nbits))))
      end do
   end do
   write(*,'(a)',advance='no') ' Enter the output file name:> '
   read(*,'(a)') filename
   open(10,file=filename,status='replace')
   do j = 1, N
      if(iand(table(j/nbits), ishft(1,mod(j,nbits))) == 0) cycle
      np = np+1
      write(10,'(1x,i0)') j
   end do
   write(*,'(a,i0,a,i0)') ' There are ', np,' primes from 2 to ', N
   write(*,'(a)') ' Compare at '// &
      'http://www.utm.edu/research/primes/howmany.shtml'
end program sieve0