! This file:
! http://ftp.aset.psu.edu/pub/ger/fortran/hdk/primesbench.f90
!
!From: "David Frank" <dave_frank@hotmail.com>
!Newsgroups: comp.lang.fortran
!Subject: Using Primes calculation  as benchmark
!Date: Sun, 26 Jan 2003 11:17:02 GMT
!
!My revised James VanBuskirk source (he posted his original in a recent
!topic) has been posted below and in comp.lang.pl1 newsgroup in
!response to a challenge and may be ofinterest here as well.
!
!
!--- CVF6.6b results running above on my 833Mhz PC ---
!
!#Primes    = 664579
!1st Prime  = 2
!Last Prime = 9999991
!time(sec)  = 2.35
!-------------------------------------
module primes_1
implicit none
contains
subroutine compute_primes(n, numbers)
implicit none
integer, intent(in) :: n
integer count
integer, allocatable :: numbers(:)
integer, parameter :: NB = bit_size(1)  ! usually = 32
integer :: sieve(0:(n-1)/(2*NB)+1)    ! bit table for prime candidates
integer :: i, j, upper

sieve = 0
sieve(0) = 1  ! 1 isn't prime
upper = sqrt(n+0.5d0)
do i = 1, (upper-1)/2
   if (.not.btest(sieve(i/NB),i-i/NB*NB)) then
      do j = 2*i*(i+1), (n-1)/2, 2*i+1
!         sieve(j/NB) = ibset(sieve(j/NB),j-j/NB*NB)
         sieve(j/NB) = ibset(sieve(j/NB),modulo(j,NB))
!         sieve(j/NB) = ior(sieve(j/NB),ishft(1,modulo(j,NB)))
      end do
   end if
end do
count = 1        ! 1st prime =2 count
do i = 1,(n-1)/2     ! All entries are odd
   if (.not.btest(sieve(i/NB),i-i/NB*NB)) count = count+1
end do
allocate( numbers(count) )
count = 1 ; numbers(1) = 2   ! preset 1st prime
do i = 1,(n-1)/2     ! All entries are odd
   if (.not.btest(sieve(i/NB),i-i/NB*NB)) then
      count = count+1
      numbers(count) = 2*i+1
   end if
end do
end subroutine compute_primes
end module primes_1
! ---------------------------------------
program primes
use primes_1
implicit none
integer,allocatable :: numbers(:)
integer :: clock1,clock2, rate, n = 100000000

call system_clock(clock1,rate)
call compute_primes(n,numbers)
call system_clock(clock2,rate)

write(*,91) ' #Primes    = ', size(numbers), &
            ' 1st Prime  = ', numbers(1), &
            ' Last Prime = ', numbers(size(numbers)), &
            ' Last Test# = ', n, &
            ' time(sec)  = ',(clock2-clock1)/dble(rate)
91 format (4(a,i0/),a,f0.2)
end program primes