Listing 1.
!================== Seive ===========
program Sieve_Of_Eratosthenes
! Find prime numbers using array processing.
!
! Strike the Twos and strike the Threes
! The Sieve of Eratosthenes!
! When the multiplies sublime
! The numbers that are left, are prime.
! From: Drunkard's Walk, by Frederik Pohl
integer :: last_number
integer, dimension(:), allocatable :: numbers
integer :: i, number_of_primes, ac
do
print *, "What is the last number you "&
"want to check?"
read *, last_number
select case (last_number)
case (0)
exit ! zero ends the testing
case (:-1, 1, 2)
print *, "That's not possible, "&
"try again"
case (3:100000)
allocate (numbers(last_number))
! Initialize numbers array to 0, 2, 3,
! ..., last_number
! Zero instead of 1 because 1 is a
! special case for primes.
numbers = (/ 0, (ac, ac = 2, &
last_number) /)
do i = 2, last_number
! if this number is prime, eliminate
! all multiples
if (numbers(i) /= 0) then
numbers(2*i : last_number : i) = 0
endif
enddo
! Count the primes.
number_of_primes = count (numbers /= 0)
! Gather them into the front of the
! array.
numbers(1:number_of_primes) = &
pack(numbers, numbers /= 0)
! Print them
print *, "There are ",number_of_primes,&
" prime numbers less than or equal to",&
last_number
print "(5i7)",&
numbers(1:number_of_primes)
deallocate (numbers)
case default
print *, "That's too large a value "&
"to try"
end select
end do
! Sample output:
! There are 25 prime numbers less than or equal
! to 100
! 2 3 5 7 11
! 13 17 19 23 29
! 31 37 41 43 47
! 53 59 61 67 71
! 73 79 83 89 97
end program Sieve_Of_Eratosthenes