[Shootout-list] Re: nsieve-bits in Fortran

[Janne Blomqvist]

On Tue, Mar 01, 2005 at 09:57:55PM +0200, Janne Blomqvist wrote:
> Hello,
> 
> here is nsieve-bits in Fortran 95 (+ getarg extension). Performance is
> about equal to the C version.

Sorry about following up my own message. Anyway, some small tweaks
made the program about 9 % faster. The array now lies on the stack,
usual caveats about stack limits apply.



-- 
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 num, m
  character(len=8) argv

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

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

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

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


contains

  function nsib (m)
    integer, parameter :: bs = bit_size (1)
    integer, intent(in) :: m
    integer, dimension(m/bs + 1) :: 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