Recursion – The Sieve of Eratosthenes

We have looked at the “sieve” before, and realistically it is actually a good candidate for recursion.

Below is the main (calling) program in Fortran. The program uses a sieve which is stored as a dynamic array, by declaring the array primes as allocatable. This makes it more efficient from a storage perspective. The user can input a value for the upper bound to check for primes, and allocate() is called to create the associated resources for prime. After the recursive subroutine sieve() is called, the remainder of the program deals with printing the primes to a file, primes.txt, in an orderly, column-based manner.

program eratosthenes

    integer :: i, j, N
    integer, dimension(:), allocatable :: primes

    open(unit=9,file='primes.txt',status='replace',action='write')

    write (*,*) "Enter the boundary to check for primes: "
    read (*,*) N

    allocate(primes(N))

    do i = 1, N
        primes(i) = i
    end do

    call sieve(N,primes,2)

    write (*,50) "Prime numbers from 1 to ", N, " written to file primes.txt"
    50 format (A,I8,A)
    j = 1
    do i = 1, N
        if (primes(i) /= 0) then
            write (9,"(I5)",advance="no") primes(i)
            j = j + 1
            if (mod(j,10) == 0) then
              write (9,*)
            end if
        end if
    end do

    close (9,status='keep')

end program eratosthenes

Now, on to the recursive subroutine sieve(). The subroutine has three parameters: N, the size of the sieve, p, the sieve array holding the primes, and x, the starting value, in this case 2.

recursive subroutine sieve(N, p, x)

    integer, intent(in) :: N
    integer, intent(inout), dimension(N) :: p
    integer, intent(in) :: x
    integer :: j, k

    k = sqrt(real(N))

    if (x == 2) then
       do j = 4, N, 2
          p(j) = 0
       end do
       call sieve(N,p,x+1)
    elseif (mod(x,2) == 1 .and. x <= k) then
       do j = x, N, x
          !if (mod(p(j),x) == 0 .and. p(j) /= x) then
          if (p(j) /= x) then
             p(j) = 0
          end if
       end do
       call sieve(N,p,x+2)
    end if

end subroutine sieve

What is happening in this subroutine?

• On Line 10, if the value of x is 2, then all even numbers except 2 are set to 0. After this, sieve() is called with the x+1. This section of code is only actioned in the initial call to sieve().
• On Line 15, the odd values are processed, until x is <= k. All multiples of j are marked in p by setting them to zero. For example if x = 3, then 6, 9, 12, etc. are set to 0. At the end of this section of code (Line 22), sieve() is called with x+2, making sure only odd numbers are processed.

Here is the output from the program:

    1    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  101  103  107
  109  113  127  131  137  139  149  151  157  163
  167  173  179  181  191  193  197  199  211  223
  227  229  233  239  241  251  257  263  269  271
  277  281  283  293  307  311  313  317  331  337
  347  349  353  359  367  373  379  383  389  397
  401  409  419  421  431  433  439  443  449  457
  461  463  467  479  487  491  499  503  509  521
  523  541  547  557  563  569  571  577  587  593
  599  601  607  613  617  619  631  641  643  647
  653  659  661  673  677  683  691  701  709  719
  727  733  739  743  751  757  761  769  773  787
  797  809  811  821  823  827  829  839  853  857
  859  863  877  881  883  887  907  911  919  929
  937  941  947  953  967  971  977  983  991  997