[Shootout-list] nsieve-bits in Fortran

[Janne Blomqvist]

Hello,

here is nsieve-bits in Fortran 95 (+ getarg extension). Performance is
about equal to the C version.

PS: I couldn't subscribe to the list, pressing the subscribe button on
the webpage only led to a 404 error page..


   
-- 
Janne Blomqvist
-------------- next part --------------
! The Great Computer Language Shootout
! http://shootout.alioth.debian.org/
!
! nsieve-bits 
!
! Janne Blomqvist 2005-03-01, based on nsieve by Simon Geard
!
! Building info.
! ==============
!
! Linux  - using the Intel Fortran compiler:
!
!          ifort nsieve-bits.f90 -O3 -static-libcxa -o nsieve-bits
!
!        - using gfortran:
!
!          gfortran nsieve-bits.f90 -O3 -o nsieve-bits

program nsievebits

  implicit none
  integer, parameter :: bs = bit_size (1)
  integer num, m
  character(len=8) argv
  integer, dimension(:), allocatable :: flags

  call getarg(1,argv)
  read(argv,*) num

  m = (2**num) * 10000
  allocate (flags(m/bs+1))
  write(*,'(A,I8,A,I8)')  'Primes up to ',m,' ',nsib (m,flags)

  m = (2**(num-1)) * 10000
  write(*,'(A,I8,A,I8)')  'Primes up to ',m,' ',nsib (m,flags)

  m = (2**(num-2)) * 10000
  write(*,'(A,I8,A,I8)')  'Primes up to ',m,' ',nsib (m,flags)


contains

  function nsib (m, prime)
    integer, intent(in) :: m
    integer, dimension(:) :: prime
    integer :: nsib, i, j

    forall (i = 1: m/bs + 1)
       prime(i) = not (ishft (prime(i), bs))
    end forall

    nsib = 0

    do i = 2, m
       if (btest (prime(i/bs + 1), mod (i, bs))) then
          do j = i + i, m, i
             prime(j/bs + 1) = ibclr (prime (j/bs + 1), mod (j, bs))
          end do
          nsib = nsib + 1
       end if
    end do

  end function nsib

end program nsievebits