# 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: "

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
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```

This site uses Akismet to reduce spam. Learn how your comment data is processed.